基于Python实现植被覆盖度时空变化规律分析实例【100010990】

# 植被覆盖度时空变化规律分析实例 ## 一、数据集 ISLSCP II GIMMS 月度 NDVI 数据集：[点我](https://daac.ornl.gov/ISLSCP_II/guides/gimms_ndvi_monthly_xdeg.html) * 涵盖了 1981.07 ——2002.12 的数据； * 该数据集包含三个数据文件，分别以 0.25 、0.5 、1.0 度的经纬度空间分辨率提供； * 数据集进行了处理，减少了因校准、视图几何、火山气溶的火山平流层气溶胶矫正，以及使用经验模式分解/重建（ EMD ）改进 NDVI 以最小化轨道漂移的影响； * 数据集空间范围：全局网格 * 经度范围：-180 W —— 180 E ； * 纬度范围：-90 S —— 90 N ； * 数据获取地址：[点我](https://search.earthdata.nasa.gov/search/granules?p=C179002875-ORNL_DAAC&pg[0][v]=f&pg[0][gsk]=-start_date&q=gimms&tl=1637818931.017!3!!) ； * 文件类型：.asc 文件； * 本实例选择了精度为 0.25 度的数据集； * 数据中含有负值，-99 表示水体，-88 表示缺失值，-77 表示永久冰值。 Arcgis 的学习建议参考官方文档：[点我](https://resources.arcgis.com/zh-cn/help/) 。 原始数据展示：  ## 二、数据处理流程 ### 1. 下载文件完整性验证 由于笔者还下载了 MODIS 的数据集，而该数据集是逐条下载的，因此在下载完成后需要判断数据集是否完整。 所以笔者写了一个小程序进行验证，有兴趣的可以看看： ```python # coding=utf-8 import os # check if the file exists f = open(r'D:\Documentation\Project\Grassland ecology\data\MODIS NDVI\download.txt') dir = r'D:\Documentation\Project\Grassland ecology\data\MODIS NDVI\hdf' path = r'D:\Documentation\Project\Grassland ecology\data\MODIS NDVI\check.txt' line = f.readline() save = open(path, mode='w') count = 0 while line: if os.path.exists(dir + os.sep + line.split('/')[-1].replace('\n', '')) == False: # txt file has line breaks, we should ignore it save.writelines(line) count+=1 print count line = f.readline() print 'ok!' ``` 解释一下，上述代码中读取了一个文本类型的文件，该文件中存储了所有数据的下载地址，例如： `https://e4ftl01.cr.usgs.gov//DP131/MOLT/MOD13A3.061/2021.10.01/MOD13A3.A2021274.h24v06.061.2021306183117.hdf https://e4ftl01.cr.usgs.gov//DP131/MOLT/MOD13A3.061/2021.10.01/MOD13A3.A2021274.h28v04.061.2021306183118.hdf https://e4ftl01.cr.usgs.gov//DP131/MOLT/MOD13A3.061/2021.10.01/MOD13A3.A2021274.h26v07.061.2021306182609.hdf` 所以这一步需要做的就是与下载好的文件进行对比，看是否有遗漏，有的话就将未下载的链接存储到一个叫做 check.txt 的文件中，这样我就可以下载它们了。 ### 2. .asc 格式 转 .tif 格式 在后续步骤的处理中，.asc 文件较为不便，因此需要将其转换为栅格类型 TIFF ： 可以使用 Arcgis 工具箱中的 `Conversion Tools/To Raster/ASCII to Raster` 工具转换格式，也可使用代码进行批量处理： ```python #encoding:utf-8 import os import arcpy dir = 'D:/Documentation/Project/Grassland ecology/Grassland_ecology/data/gimms_ndvi_qd_1981-2002_exp/' # 数据集 files = os.listdir(dir) for f in files: if os.path.splitext(f)[1] == '.asc': # Script arguments... input_raster_file = dir + f # 文件名 # output_data_type = "DOUBLE" # 数据类型 # Local variables... raster_format = "TIFF" # 文件类型 output_workspace = 'D:/Documentation/Project/Grassland ecology/Grassland_ecology/data/conversion/' + os.path.splitext(f)[0] # 输出路径 os.mkdir(output_workspace) # file name process output_raster = output_workspace + os.sep + os.path.splitext(f)[0] + ".tif" # 输出文件名 if os.path.exists(output_raster) == False: # print input_raster_file # Process: Raster To Other Format (multiple)... arcpy.RasterToOtherFormat_conversion(input_raster_file, output_workspace, raster_format) # print output_raster print "ALL DONE" ``` 转换后的数据展示：  可以发现没什么区别。 ### 3. 植被覆盖度 植被覆盖度的计算通常都使用像元二分法来计算： 假设每个像元（像素）的 NDVI 值由植被和土壤两部分组成： $$ NDVI=NDVI_vC_i+NDVI_s(1-C_i) $$ 其中，$NDVI_v$ 表示植被覆盖部分的值，$NDVI_s$ 表示土壤部分的值，$C_i$ 表示植被覆盖度： $$ C_i=\frac{NDVI-NDVI_s}{NDVI_v-NDVI_s} $$ 由于植被和土壤部分的 NDVI 值并不能确定，因此通常使用累积概率为 95% 和 5% 来代替，因此公式如下： $$ C_i=\frac{NDVI-NDVI_{0.05}}{NDVI_{0.95}-NDVI_{0.05}} $$ 下面进行实操，首先提取 NDVI 的累积概率，笔者最开始使用 ENVI 软件的统计功能进行提取，但是因为数据中含有负值，因此统计的结果不正确，而笔者又不知道该如何处理这种问题，因此最终还是使用 Python 代码处理了。 笔者一开始写代码处理计算累积概率的时候使用的是原始数据类型 .asc ，其实处理 .tif 更加简便，这里给出 .asc 文件类型的处理代码，后面在对数据进行再分类的时候会提供处理 .tif 文件类型的代码，有需要的朋友自行改写吧。 ```python # coding=utf-8 import pandas as pd import numpy as np def get_confidence(filepath): ASCfile = pd.read_csv(filepath, skiprows=6, engine='python', sep=' ', delimiter=None, index_col=False, header=None, skipinitialspace=True) ndarray = ASCfile.as_matrix().reshape(1, -1) ndarray.sort() ndarray = ndarray[ndarray>=-1] # 取正常的NDVI值 q_5 = np.percentile(ndarray, 5) q_95 = np.percentile(ndarray, 95) return q_5, q_95 ``` 植被覆盖度具体的计算可以使用 Arcgis 软件中的栅格计算器：`Spatial Analyst Tools/Map Algebra/Raster Calculator` 。 由于在计算之前需要处理一下数据中的负值，而且对于超出累积概率范围的值也需要清洗，因此处理上文中给出的计算公式外，还需要增加一些步骤： * 如果 $NDVI<NDVI_{0.05}$ ，则 $NDVI=0$ ； * 如果 $NDVI>NDVI_{0.95}$ ，则 $NDVI=1$ ； * 如果 $NDVI_{0.05}\le NDVI\le NDVI_{0.95}$ ，则 $NDVI=NDVI$ 。 因此使用栅格计算器的时候所需的公式为：`(b1 lt ndvi_min)*0 + (b1 gt ndvi_max)*1 + (b1 ge ndvi_min and b1 le ndvi_max)*((b1 + ndvi_min) / (ndvi_max - ndvi_min))` 。 使用栅格计算器的时候，一次只能处理一个文件，因此需要自行提取累积概率并输入计算器，这样会非常麻烦，因此下面给出批量处理的代码： ```python # coding=utf-8 import os import arcpy from arcpy.sa import * from get_confidence import get_confidence arcpy.CheckOutExtension("spatial") # 检查外部扩展 arcpy.gp.overwriteOutput = 1 # 覆盖之前的文件 dirs = 'D:/Documentation/Project/Grassland ecology/Grassland_ecology/data/conversion/' out_dir = 'D:/Documentation/Project/Grassland ecology/Grassland_ecology/data/Vegetation_coverage/' asc_dir = 'D:/Documentation/Project/Grassland ecology/Grassland_ecology/data/gimms_ndvi_qd_1981-2002_exp/' files = os.listdir(dirs) for file in files: if os.path.exists(out_dir + file) == False: os.mkdir(out_dir + file) out_path = out_dir + file + os.sep dir = dirs + file + os.sep arcpy.env.workspace = dir rasters = arcpy.ListRasters(raster_type='TIF') for raster in rasters: inRaster = Raster(raster) Min, Max = get_confidence(asc_dir + file + ".asc") ans = Con(inRaster < Min, 0, Con((inRaster >= Min) & (inRaster <= Max), (inRaster - Min)/(Max - Min), 1)) ans.save(out_path + file + ".tif") print file + " is done!" print "OK!" ``` 植被覆盖度图像如下：  可以发现，它的颜色变浅了一些。 ### 4. 截取目标区域 实例以中�