[MODIS数据处理#9]例四：基于MCD12Q2数据集初步分析中国植被物候空间分布特征

最新推荐文章于 2024-04-18 19:00:00 发布
最新推荐文章于 2024-04-18 19:00:00 发布
阅读量5.3k 收藏 108
点赞数 19
分类专栏： MODIS数据处理 文章标签： python MCD12Q2 MODIS
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/qq_37948866/article/details/113819125
版权

文章目录


系列文章目录： MODIS数据处理系列专栏

本文基于Arcmap10.2对MCD12Q2植被物候生长数据集进行处理，并在每个步骤中提供批量处理的python代码。

一、MCD12Q2介绍

MCD12Q2以MODISEVI为数据源，先对 EVI 时间序列数据进行预处理，消减云、气溶胶、雪和观测角度等影响，然后用滑动窗口（连续 5 个时相 EVI 值）来判断持续升高和降低区间，当所在区间的极大值和振幅满足给定阈值条件后，判断该升高和降低区间为一个生长周期过程。对于每一个生长周期，采用分段 Logistic 函数拟合，并利用其曲率值变化的极值点，确定生长始期、EVI 持续增加的中点、成熟期、峰值、末期、 EVI 持续降低中点和生长末期。一年当中最多记录两个生长周期。详情介绍可参考：MCD12Q2数据集说明文档
在这里插入图片描述
MCD12Q2数据集中各子数据集对应的变量如下图所示，其中各物候节点的值为距1970年1月1日的天数：

  • Greenup：生长季始期
  • MidGreenup：EVI持续增加的中点
  • Peak：峰值
  • Maturity：成熟期
  • MidGreendown：EVI持续降低的中点
  • Senescence：衰老
  • Dormancy：生长季末期

在这里插入图片描述

二、问题描述

本文演示使用的是MCD12Q2数据集，空间分辨率为500m，时间分辨率为yearly。为了分析中国区域植被物候空间分布特征，需要提取生长季始期SOS（start of the growing season）、生长季末期EOS(end of the growing season)，这两个变量分别对应MCD12Q2的Greenup、Dormancy，并利用生长季末期SOS减生长季始期EOS得到生长季长度LOS(length of the growing season)

三、预处理

预处理涉及到的操作均可在Arcmap中完成，演示的版本为10.2，每步均提供有批量处理的方法。

1、提取子数据集

利用提取子数据集工具从.hdf文件中提取Greenup、Dormancy子数据集并转为栅格文件。
在这里插入图片描述
批量处理的方法见：【ArcGIS自定义脚本工具】批量提取子数据集
经过此步我们得到了Greenup、Dormancy子数据集的栅格（.tif）文件。
在这里插入图片描述

2、波段提取

MCD12Q2一年当中最多记录两个生长周期，中国部分区域存在两个生长周期，且大都分布在华东地区，这在EVI的全年变化曲线中也有体现，如下：
在这里插入图片描述
这导致提取后的子数据集波段数为2，为了方便分析及之后的操作，我们需要将栅格文件的两个波段分别提取出来。
在这里插入图片描述
利用创建栅格图层工具将Greenup、Dormancy子数据集中的两个波段分别存储为栅格，并以后缀“_b1”与“_b2”区别所在波段，参数设置可参考下图。
在这里插入图片描述
其他提取波段的方法可以参考这篇：利用ArcGIS10.2波段合并及提取。基于ArcGIS的自定义脚本，可以批量提取栅格某个波段，代码如下：

#!/usr/local/bin/python2.7
# -*- coding: utf-8 -*-

"""
@File    : Batch_extractband.py
@Author  : salierib
@Time    : 2020/12/15 22:19
"""
import arcpy
import os
import time

raster_path = arcpy.GetParameterAsText(0)
band_index = arcpy.GetParameterAsText(1)
out_path = arcpy.GetParameterAsText(2)

arcpy.env.workspace = raster_path
rasters = arcpy.ListRasters("*", "tif")

nums = len(rasters)  # number of hdf files
num = 1  # current serial number
if not os.path.exists(out_path):
    os.makedirs(out_path)
for raster in rasters:
    s = time.time()
    new_name = "{0}_b{1}.tif".format(raster, band_index)
    try:
        arcpy.MakeRasterLayer_management(raster, new_name, "#", "#", band_index)
        arcpy.CopyRaster_management(new_name, os.path.join(out_path, new_name), "", "", "", "NONE", "NONE", "")
        e = time.time()
        arcpy.AddMessage("{0}/{1} | {2} Completed, time used {3}s".format(num, nums, new_name, e-s))
    except:
        arcpy.AddMessage("{0}/{1} | {2} Errored".format(num, nums, new_name))
    num += 1

该代码可以将raster_path中的所有栅格（.tif）中波段索引为band_index的波段提取出来，并存储到out_path路径下，提取后的文件名为“原文件名+_波段索引”，工具参数设置及脚本运行截图如下所示：
在这里插入图片描述
还可以利用QGIS右键重排波段工具，以批处理的方式运行，具体步骤如下：
在这里插入图片描述
利用表达式对转换后的文件进行命名，表达式为：
format('%1\\%2_band1.tif',file_path(@INPUT),base_file_name(@INPUT))
在这里插入图片描述
在这里插入图片描述

3、拼接

利用镶嵌至新栅格工具，将同一波段、同时刻不同区块的多个栅格拼接为一个栅格，参数设置可参考：
在这里插入图片描述
批量镶嵌不同MODIS区块栅格的方法可以参考：【ArcGIS自定义脚本工具】按区块批量镶嵌MODIS影像

4、栅格投影

栅格目标坐标系采用Albers等面积割圆锥投影方法（第1、第2纬线、中央经线分别采用27，45，105），地理坐标系选择WGS_1984，可在ArcGIS自定义投影坐标界面新建该坐标系。
在这里插入图片描述
新建适合中国地区的自定义投影坐标系后，利用栅格投影工具将拼接好的包含中国地区的栅格文件投影到目标坐标系，参数设置可参考：
在这里插入图片描述
批量投影栅格到指定坐标系的方法可以参考：【ArcGIS自定义脚本工具】批量重投影栅格脚本

5、裁剪

经过上述1-4步骤，我们得到了包含中国区域的Greenup、Dormancy不同波段的栅格文件，如下所示：
在这里插入图片描述
在这里插入图片描述
利用裁剪工具，仅保留中国区域范围的栅格，参数设置可参考：
在这里插入图片描述
在这里插入图片描述
批量裁剪栅格的方法可以参考：【ArcGIS自定义脚本工具】批量裁剪栅格脚本

四、获取生长季长度

利用数学分析工具箱的减工具用同波段索引的生长季末期Dormancy栅格减去生长季初期Greenup栅格，可以得到反映生长季长度的栅格影像，参数设置如下：
在这里插入图片描述
在这里插入图片描述

批量执行该步骤的代码可以参考：

import sys
arcpy_path = [r'C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy',
              r'C:\Program Files (x86)\ArcGIS\Desktop10.2\bin',
              r'C:\Program Files (x86)\ArcGIS\Desktop10.2\ArcToolbox\Scripts']#修改成Arcgis安装对应路径
sys.path.extend(arcpy_path)

import arcpy
from arcpy import env
from arcpy.sa import *
import os
import datetime
import time

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

# Set environment settings
path1 = r"F:\MCD12Q2v006\start" # 输入，Greenup所在的文件夹
path2 = r"F:\MCD12Q2v006\end" # 输入，Dormancy所在的文件夹
out_path = r"F:\MCD12Q2v006\LOS" # 输出
env.workspace = path1
rasters = arcpy.ListRasters('*', 'tif')

if not os.path.exists(out_path):
    os.makedirs(out_path)
nums = len(rasters)
num = 1

for raster in rasters:
    # start stopwatch
    s = time.time()

    inRaster1 = raster # china_greenup_a2004001_b1.tif
    inRaster2 = os.path.join(path2, raster.replace("greenup", "dormancy")) # china_dormancy_a2004001_b1.tif

    # Execute Minus
    outMinus = Minus(inRaster2, inRaster1)

    # Save the output
    newtif = "LOS" + "_" + raster.split("_")[2]
    newraster = os.path.join(out_path, newtif)
    outMinus.save(newraster)

    # stop stopwatch
    e = time.time()
    arcpy.AddMessage("{0}/{1} | {2} Completed, time used {3}s".format(num, nums, newtif, e-s))
    num += 1

五、获取生长季始、末所在日期

通过上述的预处理操作我们得到了中国区域逐年Greenup（生长期初期）、Dormancy（生长季末期）的栅格，但此时栅格的值代表的是距1970年1月1日的天数。例如值为12518时，对应的日期为2004年4月10日。
在这里插入图片描述
为了方便分析，将栅格的值转换为日期在所在年的日序，2004年4月10日是2004年的第100天，故需将像元值从12518转为100。
在这里插入图片描述
此步可以利用数学分析工具箱的减工具完成
在这里插入图片描述
批量执行该步骤的代码如下：

#!/usr/local/bin/python2.7
# -*- coding: utf-8 -*-

"""
@File    : arcpy_minus.py
@Author  : salierib
@Time    : 2020/12/16 20:02
"""
"""
本脚本可以通过在pycharm中直接调用Arcpy的方式运行
作用：将MCD12Q2中各物候阶段对应的原始值（距1970-1-1的天数），转为在当前年的日序，即生长季开始的日子位于当前年的日序
"""
import sys
arcpy_path = [r'C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy',
              r'C:\Program Files (x86)\ArcGIS\Desktop10.2\bin',
              r'C:\Program Files (x86)\ArcGIS\Desktop10.2\ArcToolbox\Scripts']#修改成Arcgis安装对应路径
sys.path.extend(arcpy_path)

import arcpy
from arcpy import env
from arcpy.sa import *
import os
import datetime
import time

# Set environment settings
in_path = r"F:\MCD12Q2v006\6_China" # 输入，待转换后栅格所在的文件夹
out_path = r"F:\MCD12Q2v006\7_days" # 输入，存储转换后栅格的文件夹，若无，会自动创建
env.workspace = in_path
rasters = arcpy.ListRasters('*', 'tif')

if not os.path.exists(out_path):
    os.makedirs(out_path)
nums = len(rasters)
num = 1
# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")
for raster in rasters:
    # start stopwatch
    s = time.time()

    # Set local variables
    inRaster1 = raster
    year = '{0}-01-01'.format(int(raster.split("_")[-2][1:5])) # 输入栅格的原始命名需满足类似"xxx_a2004001_xx.tif"
    d1 = datetime.datetime.strptime(year, '%Y-%m-%d')
    d2 = datetime.datetime.strptime('1970-01-01', '%Y-%m-%d')
    t = d1 - d2
    inRaster2 = t.days

    # Execute Minus
    outMinus = Minus(inRaster1, inRaster2)

    # Save the output
    newraster = os.path.join(out_path, "days" + "_" + "_".join(raster.split("_")[3:]))
    outMinus.save(newraster)

    # stop stopwatch
    e = time.time()
    arcpy.AddMessage("{0}/{1} | {2} Completed, time used {3}s".format(num, nums, raster, e-s))
    num += 1

可将上述代码保存为.py文件后在Arcmap中的python窗口加载该文件运行，具体操作步骤如下：
在这里插入图片描述

在这里插入图片描述
在这里插入图片描述
处理后的栅格文件如下图所示：
在这里插入图片描述

Salierib
关注
  • 19
    点赞
  • 踩
  • 108
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 9
    评论
专栏目录
GEE：植被生长季节开始和结束时间DOY（MCD12Q2数据集Greenup_1 和 Dormancy_1波段）
养乐多的博客
04-06 2297
本文将介绍Google Earth Engine（GEE）平台上的MODIS/006/MCD12Q2数据集的Greenup_1和Dormancy_1波段的数据转换方法和代码。数据的像素值是从1970年开始到至今的天数，需要转换到指定年份的Day of Year（DOY）才可以使用，本文记录了转换方法和代码。MODIS/006/MCD12Q2数据集的Greenup_1和Dormancy_1波段分别反映了地表植被的绿度增长时间和休眠时间。
MODIS 2010年中国1km植被指数（NDVI）空间分布数据集
03-18
该植被指数数据来源于美国NASA定期发布的MODIS数据系列，原始数据集为MOD13A3。MOD13A3的空间分辨率为1km，时间分辨率为monthly。经过提取子数据集、拼接、投影栅格、换算单位、裁剪得到逐月1km的NDVI数据后，再利用最大合成法得到对应年份的NDVI数据。 地区：中国 时间分辨率：年 空间分辨率：1km 投影坐标系：Albers conical equal area 椭球：WGS84 变形比例：1.0 中央经线：105 标准纬线：25 标准纬线：47 引用：Didan, K. (2015). MODIS/Terra Vegetation Indices Monthly L3 Global 1km SIN Grid V006, NASA EOSDIS Land Processes DAAC.https://doi.org/10.5067/MODIS/MOD13A3.006该植被指数数据来源于美国NASA定期发布的MODIS数据系列，原始数据集为MOD13A3。MOD13A3的空间分辨率为1km，时间分辨率为monthly。经过提取子数据集、拼接、投影栅格、换算单位、
【RS&GIS数据资源】2000‐2019 年全国三大粮食作物 1km 物候数据集
最新发布
Lee大哲博客
04-18 1107
2000-2019 年全国三大粮食作物 1km 物候数据集（ChinaCropPhen1km）是基于 GLASS LAI 遥感产品，利用农业气象站长期观测的物候记录，根据作物物候曲线特征定义了每种作物的三个关键物候期的拐点和峰值，结合拐点法和阈值法识别出 LAI 特征曲线对应的拐点和峰值，最终得到每个作物像元的关键物候期在一年中出现的日期（儒略日，DOY）
Google Earth Engine ——MCD12Q2 V6土地覆盖动态产品（非正式地称为MODIS全球植被表征产品）提供全球范围内的植被表征时间估计
qq_31988139的博客
10-10 2209
The MCD12Q2 V6 Land Cover Dynamics product (informally called the MODIS Global Vegetation Phenology product) provides estimates of the timing of vegetation phenology at global scales. Additionally, it provides information related to the range and summation
GEE：植被生长季长度（MCD12Q2数据集）
养乐多的博客
04-06 1439
本文将介绍Google Earth Engine（GEE）平台上的MODIS/006/MCD12Q2数据集的Greenup_1和Dormancy_1波段计算植被生长季长度的代码，原始数据的像素值是从1970年开始到至今的天数，需要转换到指定年份的Day of Year（DOY）才可以使用，本文记录了转换方法和代码。MODIS/006/MCD12Q2数据集的Greenup_1和Dormancy_1波段分别反映了地表植被的绿度增长时间和休眠时间。 它们的差表示生长季长度length。
GEE开发之MODIS_MCD12Q1数据分析和获取
等待着冬天的风的博客
05-09 3769
Terra 和 Aqua 组合的中分辨率成像光谱仪 (MODIS) 土地覆盖类型 (MCD12Q1) 版本 6.1 数据产品每年提供全球土地覆盖类型。MCD12Q1 版本 6.1 数据产品是使用 MODIS Terra 和 Aqua 反射率数据的监督分类得出的。土地覆盖类型源自国际地圈-生物圈计划 (IGBP)、马里兰大学 (UMD)、叶面积指数 (LAI)、生物群落-生物地球化学循环 (BGC) 和植物功能类型 (PFT) 分类方案。
MODIS数据知识积累
GeoSuper的博客
10-17 7825
@[TOPIC] (MODIS数据积累) MODIS数据使用的投影 sinusoidal 投影
MODIS 及其 数据、产品
qq_40292148的博客
04-06 2587
MODIS 是美国宇航局研制的一种大型空间，全称为 MoDerate Resolution Imaging Spectroradiometer。它搭载在两颗卫星上：Terra（EOS AM-1）和 Aqua（EOS PM-1）。
MODIS 2023年中国250m植被指数（NDVI）空间分布数据集
03-24
该植被指数数据来源于美国NASA定期发布的MODIS数据系列，原始数据集为MOD13Q1。MOD13Q1的空间分辨率为250m，时间分辨率为16d。经过提取子数据集、拼接、投影栅格、换算单位、裁剪得到16d 250m的NDVI数据后，再利用...
MODIS 2020年逐月中国1km植被指数（NDVI）空间分布数据集
03-18
该植被指数数据来源于美国NASA定期发布的MODIS数据系列，原始数据集为MOD13A3。MOD13A3的空间分辨率为1km，时间分辨率为monthly。经过提取子数据集、拼接、投影栅格、换算单位、裁剪得到逐月1km的NDVI数据。 地区：...
MODIS 2020年中国250m植被指数（NDVI）空间分布数据集
04-27
该植被指数数据来源于美国NASA定期发布的MODIS数据系列，原始数据集为MOD13Q1。MOD13Q1的空间分辨率为250m，时间分辨率为16d。经过提取子数据集、拼接、投影栅格、换算单位、裁剪得到16d 250m的NDVI数据后，再利用...
MCD12Q2_Collection6_UserGuide.pdf
05-15
MODIS土地覆盖数据 MCD12Q2用户指南。资料格式：pdf。发布时间：2019年5月15日。包括内容：数据概况，方法，数据格式和投影，数据说明等。
Modis MCD12Q1数据提取工具
03-17
批量提取MODIS土地覆盖数据MCD12Q1中的数据集。使用方法：选择输入目录下的一个hdf文件，选择要提取的字段，指定输出目录，点击运行即可将输入目录下所有hdf文件中MCD12Q1数据集中指定的数据集提取出来，结果按按照原文件名，保存为tiff格式。参考：http://blog.csdn.net/giselite/article/details/21081297
MODIS 2020年中国1km植被指数（NDVI）空间分布数据集
03-18
该植被指数数据来源于美国NASA定期发布的MODIS数据系列，原始数据集为MOD13A3。MOD13A3的空间分辨率为1km，时间分辨率为monthly。经过提取子数据集、拼接、投影栅格、换算单位、裁剪得到逐月1km的NDVI数据后，再利用...
MODIS 2019年中国1km植被指数（NDVI）空间分布数据集
03-18
该植被指数数据来源于美国NASA定期发布的MODIS数据系列，原始数据集为MOD13A3。MOD13A3的空间分辨率为1km，时间分辨率为monthly。经过提取子数据集、拼接、投影栅格、换算单位、裁剪得到逐月1km的NDVI数据后，再利用...
【ArcGIS自定义脚本工具】批量提取子数据集
Salierib的博客
09-02 2687
文章目录一、功能介绍二、脚本代码三、工具参数四、工具界面 一、功能介绍 二、脚本代码 #!/usr/bin/python # -*- coding: UTF-8 -*- import os import arcpy import time in_path = arcpy.GetParameterAsText(0) out_path = arcpy.GetParameterAsText(1) subdataset_index = arcpy.GetParameterAsText(2) suffix
使用MATLAB批量将MCD12Q2物候产品数据初始值转换为儒略日
qq_44996109的博客
07-13 256
适用于已进行了预处理的影像
基于MODIS数据的2000-2021年中国植被生物量数据集
qq_31988139的博客
09-22 759
中国植被生物量分布图集由航天宏图实验室提供，利用MODIS地表反射率数据和植被指数（MCD43A4、MCD12Q1），与美国农业部FIA项目提供的地上生物量数据结合进行随机森林机器学习，训练得到各种植被类型的地上生物量估算模型。模型的总体预测精度为R2=0.72，RMSE=33.90 Mg/ha。将模型应用到中国地区，产生了覆盖全国的500米分辨率的年度植被地上生物量数据和相关图集。对森林火灾前后受灾区域的生物量分别进行估算和变化监测，能够有效评估火灾带来的植被生物量的损失，为火灾灾后评估提供有效的支撑。
[ArcGIS批处理] MOD04 | 影像分类（Arcpy）+抽取子集+投影+融合
geodoer
08-08 3789
【内容介绍】 将一个文件夹下的文件按照【一定的规则】分流到不同文件夹下 对一组数据进行批量抽取子集 对文件夹下的所有数据进行定义投影 对文件夹下的所有数据投影到指定坐标（UTM50N） 对文件夹下的所有栅格数据进行拼接 文章目录批处理【第一步】按规则分流数据【第二步】抽取子集并定义坐标系【第三步】批量投影【第四步】拼接（即ArcGIS中的镶嵌）【附录】工具详情[1] ClassfiyAnd...
MODIS/006/MCD12Q1 MODIS/061/MCD12Q1
12-28
MODIS/006/MCD12Q1和MODIS/061/MCD12Q1是MODIS的MCD12Q1产品的不同版本。MCD12Q1是MODIS的土地类型分类产品，用于对地表进行分类和监测。这些产品提供了全球范围内的土地类型信息，可以用于研究和监测地表的变化。 要使用MODIS/006/MCD12Q1和MODIS/061/MCD12Q1，您可以按照以下步骤进行操作： 1. 遥感影像查看：使用遥感影像查看工具（如ENVI、ArcGIS等）打开MODIS/006/MCD12Q1或MODIS/061/MCD12Q1数据集，以便查看地表的土地类型分类结果。 2. MCD12Q1分类介绍：了解MCD12Q1产品的分类方案和分类结果。MCD12Q1产品将地表划分为不同的土地类型，如森林、草地、农田等，并为每个土地类型分配一个唯一的标识码。 3. 年数据下载：根据您的需求，选择相应的年份和空间分辨率（500m）下载MODIS/006/MCD12Q1或MODIS/061/MCD12Q1的年度数据。这些数据可以从MODIS数据存档中获取，例如LP DAAC（Land Processes Distributed Active Archive Center）。 请注意，MODIS数据和产品通过LP DAAC获取后，对于后续的使用、销售或再分发没有任何限制。

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 9
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

Salierib

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付：¥1
获取中
扫码支付

您的余额不足，请更换扫码支付或充值

打赏作者

实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值