| 基于Python的植被覆盖度时空变化规律分析 | 您所在的位置：网站首页 › 植被变化监测 › 基于Python的植被覆盖度时空变化规律分析 |
基于Python的植被覆盖度时空变化规律分析
|
资源下载地址：https://download.csdn.net/download/sheziqiong/85745576 资源下载地址：https://download.csdn.net/download/sheziqiong/85745576 一、数据集ISLSCP II GIMMS 月度 NDVI 数据集：点我 涵盖了 1981.07 ——2002.12 的数据；该数据集包含三个数据文件，分别以 0.25 、0.5 、1.0 度的经纬度空间分辨率提供；数据集进行了处理，减少了因校准、视图几何、火山气溶的火山平流层气溶胶矫正，以及使用经验模式分解/重建（ EMD ）改进 NDVI 以最小化轨道漂移的影响；数据集空间范围：全局网格经度范围：-180 W —— 180 E ；纬度范围：-90 S —— 90 N ；数据获取地址：点我 ；文件类型：.asc 文件；本实例选择了精度为 0.25 度的数据集；数据中含有负值，-99 表示水体，-88 表示缺失值，-77 表示永久冰值。Arcgis 的学习建议参考官方文档：点我 。 原始数据展示：
由于笔者还下载了 MODIS 的数据集，而该数据集是逐条下载的，因此在下载完成后需要判断数据集是否完整。 所以笔者写了一个小程序进行验证，有兴趣的可以看看： # 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 工具转换格式，也可使用代码进行批量处理： #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 值由植被和土壤两部分组成： N D V I = N D V I v C i + N D V I s ( 1 − C i ) NDVI=NDVI_vC_i+NDVI_s(1-C_i) NDVI=NDVIvCi+NDVIs(1−Ci) 其中， N D V I v NDVI_v NDVIv 表示植被覆盖部分的值， N D V I s NDVI_s NDVIs 表示土壤部分的值， C i C_i Ci 表示植被覆盖度： C i = N D V I − N D V I s N D V I v − N D V I s C_i=\frac{NDVI-NDVI_s}{NDVI_v-NDVI_s} Ci=NDVIv−NDVIsNDVI−NDVIs 由于植被和土壤部分的 NDVI 值并不能确定，因此通常使用累积概率为 95% 和 5% 来代替，因此公式如下： C i = N D V I − N D V I 0.05 N D V I 0.95 − N D V I 0.05 C_i=\frac{NDVI-NDVI_{0.05}}{NDVI_{0.95}-NDVI_{0.05}} Ci=NDVI0.95−NDVI0.05NDVI−NDVI0.05 下面进行实操，首先提取 NDVI 的累积概率，笔者最开始使用 ENVI 软件的统计功能进行提取，但是因为数据中含有负值，因此统计的结果不正确，而笔者又不知道该如何处理这种问题，因此最终还是使用 Python 代码处理了。 笔者一开始写代码处理计算累积概率的时候使用的是原始数据类型 .asc ，其实处理 .tif 更加简便，这里给出 .asc 文件类型的处理代码，后面在对数据进行再分类的时候会提供处理 .tif 文件类型的代码，有需要的朋友自行改写吧。www.biyezuopin.vip # 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 。 由于在计算之前需要处理一下数据中的负值，而且对于超出累积概率范围的值也需要清洗，因此处理上文中给出的计算公式外，还需要增加一些步骤： 如果 N D V I < N D V I 0.05 NDVINDVI_{0.95} NDVI>NDVI0.95 ，则 N D V I = 1 NDVI=1 NDVI=1 ；如果 N D V I 0.05 ≤ N D V I ≤ N D V I 0.95 NDVI_{0.05}\le NDVI\le NDVI_{0.95} NDVI0.05≤NDVI≤NDVI0.95 ，则 N D V I = N D V I NDVI=NDVI 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)) 。 使用栅格计算器的时候，一次只能处理一个文件，因此需要自行提取累积概率并输入计算器，这样会非常麻烦，因此下面给出批量处理的代码： # 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) & (inRaster |
【本文地址】