文章目录
简介
本文介绍了利用MOD13A3数据处理加工中国各省级行政区2000年-2021年共22年1km逐月NDVI影像的过程。该处理过程以hdf文件为基础，采用本地处理的方式处理加工各个省的NDVI影像。
在1节中简单介绍了MOD13A3数据集及其下载方式。2.1节介绍了利用易方MODIS处理工具箱在Arcmap软件中基于MOD13A3数据（.hdf）加工得到宁夏2000年-2010年逐月NDVI栅格文件（.tif）的案例。2.2节介绍了固态硬盘、多线程并行在加快MODIS数据处理的应用。3节展示了对中国共34个县级行政区2000-2021年逐月NDVI的各项统计指标。4节提供了2001年各地区逐月NDVI数据和加工完成的NDVI栅格数据的下载地址。
1. 数据来源
1.1. MOD13A3数据简介
本文的植被指数（NDVI）数据来源于美国NASA定期发布的MOD13A3植被指数数据集。MOD13A3提供自2000年2月起逐月1km的NDVI数据。详细的介绍见：https://doi.org/10.5067/MODIS/MOD13A3.006
1.2. 数据下载
数据下载可以参考[MODIS数据处理#0]下个数据能再简单些吗？，还可以在EarthData（https://search.earthdata.nasa.gov/search）搜索MOD13A3进行下载，如下：
1.3. 下载后的文件大小
MOD13A3的文件大小比较小，中国区域单年的hdf文件为22(区块数)*12(月份数)=264个，大小约为4GB。中国区域2000-2021年的MOD13A3原始数据文件（hdf）大小约为87GB。
另附2021年中国区域MOD13A3数据集（.hdf）的百度网盘下载链接，https://pan.baidu.com/s/1gwipPqSsmkSaXpJBYeSY9A (提取码pkog)，可以用于练习处理hdf文件，需要的可以自取。
2. 处理方法
2.1. 处理单个省份
易方MODIS处理工具是本人基于Arcmap的自定义脚本编程制作的用于处理MODIS数据集的工具，可以将下载得到的原始hdf文件完成一系列的预处理操作，得到指定区域的栅格影像。该工具的介绍及下载地址见【MODIS数据处理#15】分享一个自制的MODIS数据处理工具箱
2.2. 加快处理速度的几个方法
用MODIS数据综合处理脚本，可以自动完成对从hdf提取NDVI、拼接区块、投影、裁剪、换算单位一系列操作，速度比MRT工具要快。如果你想进一步提高处理MODIS数据的运行速度，可以参考以下两种途径
2.2.1. 选用固态硬盘读写数据
在使用易方MODIS处理工具时，可以通过将hdf文件存放的位置以及工作空间都设置为固态硬盘来加快处理的速度。
2.2.2. 多线程处理
在pycharm调用arcpy库和multiprocessing库来实现多核并行处理，以下是多核实现投影栅格的代码示例：
# -*- coding: utf-8 -*-
"""
@File : 多核并行-重投影.py
@Author : salierib
@Time : 2022/5/15 11:07
"""
import os
import time
import multiprocessing as mp
import sys
arcpy_path = [r'C:\Program Files (x86)\ArcGIS\Desktop10.7\arcpy',
r'C:\Program Files (x86)\ArcGIS\Desktop10.7\bin',
r'C:\Program Files (x86)\ArcGIS\Desktop10.7\ArcToolbox\Scripts'] # 修改成Arcgis安装对应路径
sys.path.extend(arcpy_path)
import arcpy
def list_split(in_list, num_split):
"""
将列表划分为n份，并返回切分后的嵌套列表
当列表大小不能被num_split整除时，如23个数分为4份，则各组的数字个数为6，6，6，5
Parameters
----------
in_list:List
待切分的列表
num_split:int
需要划分的组数
Returns
-------
list_splitted:List[List[obj]]
分组得到的嵌套列表
Examples
-------
>> list_split(list(range(1,23)), num_split=4)
[[1, 2, 3, 4, 5, 6, 7],
[8, 9, 10, 11, 12, 13, 14],
[15, 16, 17, 18, 19, 20, 21],
[22, 23, 24]]
"""
lenn = len(in_list)
n = lenn // num_split + 1
list_splitted = [in_list[i * n:(i + 1) * n] for i in range((lenn + n - 1) // n)]
return list_splitted
def batch_project_raster(rasters, out_dir, prefix="pr_", out_coor_system="WGS_1984.prj",
resampling_type="NEAREST", cell_size="250 250"):
"""
批量投影栅格工具
将in_dir中的所有栅格文件投影到out_corr_syestem坐标系下
Parameters
----------
rasters:List[str]
由待进行投影栅格操作的栅格文件的绝对路径组成的列表
out_dir:str
批量投影栅格后的输出文件夹
prefix:str
投影后栅格的新文件名的前缀
out_coor_system:
待投影到的目标坐标系文件路径(.prj)
resampling_type:str
要使用的重采样算法。默认设置为 NEAREST。
NEAREST —最邻近分配法
BILINEAR —双线性插值法
CUBIC —三次卷积插值法
MAJORITY —众数重采样法
cell_size:str
新栅格数据集的像元大小。
若输出分辨率为250m，则为”250 250"
Examples
----------
>> in_dir = r"S:\1_merge"
>> tifs = [os.path.join(in_dir,n) for n in os.listdir(in_dir) if n.endswith(".tif")]
>> batch_project_raster(tifs, out_dir=r"S:\test2")
"""
if arcpy.CheckExtension("Spatial") != "Available":
print("Error!!! Spatial Analyst is unavailable")
nums = len(rasters)
num = 1
for raster in rasters:
s = time.time()
raster_name = os.path.split(raster)[1]
out_raster = os.path.join(out_dir, prefix + raster_name) # 投影后栅格的绝对路径
if not os.path.exists(out_raster):
try:
arcpy.ProjectRaster_management(raster, out_raster, out_coor_system, resampling_type, cell_size, "#",
"#", "#")
e = time.time()
print("Step3/5 %d/%d | %s completed, time used %.2fs" % (num, nums, out_raster, e - s))
except Exception as err:
print("Step3/5 %d/%d | %s errored, %s" % (num, nums, out_raster, err))
else:
print("Step3/5 %d/%d | %s already exists" % (num, nums, raster))
num = num + 1
if __name__ == "__main__":
# 线程数,可以通过mp.cpu_count()查看自己电脑的可用总线程数
num_thread = 10
# 待投影的栅格存放文件夹
in_dir = r"S:\MOD13A3_hdf\2_mosaic"
# 投影后栅格的存储文件夹
out_dir = r"S:\MOD13A3_hdf\3_pr"
# 目标投影坐标系的坐标系文件
prj_file = r"S:\utm_shps\utm48ningX.prj"
if not os.path.exists(out_dir):
os.mkdir(out_dir)
pool = mp.Pool(num_thread)
rasters = [os.path.join(in_dir, fname) for fname in os.listdir(in_dir) if fname.endswith(".tif")]
project_param_list = [(raster_set, out_dir, "Pr_", prj_file, "NEAREST", "1000 1000") for raster_set in list_split(rasters, num_thread)]
results = [pool.apply_async(batch_project_raster, args=param) for param in project_param_list]
results = [p.get() for p in results]
3. 2000年至2021年逐月NDVI变化统计结果
各省级行政区的逐月NDVI统计指标（最小值、最大值、平均值、中位数、标准差）分析结果见中国34个省级行政区2000年-2021年逐月NDVI统计分析结果
4. NDVI栅格数据下载地址
4.1. 栅格数据说明
| 属性 | 信息 |
|---|---|
| 发布日期 | 2022.5.15 |
| 地区 | 中国34个省级行政区 |
| 时间跨度 | 2000年2月至2021年12月 |
| 时间分辨率 | 逐月 |
| 空间分辨率 | 1 km |
| 投影坐标系 | WGS 1984 UTM Zone |
| 存储类型 | 栅格图像（.tif） |
| 数据大小 | 12.3GB |
对于给定的省级行政区，投影坐标系选择对应的utm分带，包含已裁剪（行政区边界裁剪得到）和未裁剪（扩大的行政区矩形边界）两部分，可按需选用。
4.2. 样例：各省级行政区2001年逐月NDVI数据下载地址
数据下载地址为https://pan.baidu.com/s/1hinbDqnvlYXMDN_77O7a-A （提取码：uifc ）。解压后的文件目录如下：
4.3. 各省级行政区2000-2021年逐月NDVI栅格数据下载
以下是本人加工的全国34个省级行政区逐月NDVI栅格数据的付费资源地址，也可以参考2.1节处理单个省市的方法自己加工下。
1407
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
