VNP22Q2 物候数据集简介及波段value转DOY代码

最新推荐文章于 2025-07-28 23:15:09 发布

LXZYJY

最新推荐文章于 2025-07-28 23:15:09 发布

阅读量566

点赞数 1

CC 4.0 BY-SA版权

文章标签： python

本文链接：https://blog.youkuaiyun.com/LLLJSJKD/article/details/133638373

本文介绍了如何使用Python处理VNP22Q2物候数据集，包括从GEE平台下载数据、波段值转换到DOY日期，以及针对EOS、SOS、POS和LOS等物候指标的计算和处理过程。作者提供了详细的步骤和代码示例，以帮助读者理解和应用这些技术。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

一、VNP22Q2 物候数据集简介

空间分辨率：500m

时间分辨率：年

波段数：每年38个波段

gee数据集链接：https://developers.google.com/earth-engine/datasets/catalog/NOAA_VIIRS_001_VNP22Q2

官方说明手册：

https://lpdaac.usgs.gov/documents/637/VNP22_User_Guide_V1.pdf

二、波段value转DOY的Python代码

本文的价值在于有关波段value转DOY的内容在gee平台无说明；而官方说明手册的相关说明不清晰，笔者也走了一圈弯路，故而在这里与大家分享(●'◡'●)

本文以EOS（对应下表 Onset_Greenness_Decrease）、 SOS（Onset_Greenness_Increase）、 POS（Onset_Greenness_Maximum）、LOS（Growing_Season_Length ）四个物候指标为例。其中LOS为SOS和EOS之差，故而不需要offset计算。因此重点是前三个指标。

上表是官方手册的截图，其中-366(given year-2000)书写十分迷惑，但经过查阅及核验发现其意思是：（直接套公式即可）

Scaled Data = Unscaled data – (given year-2000)*366

引用：Land Surface Phenology - End of growing season - First growing cycle - Datasets - "FAO catalog"

这里分享我的代码：

我的情况：从gee上将2013-2019的数据堆叠为一张图下载，总共38*7=266个波段，目标提取EOS、 SOS、 POS、LOS四个物候指标。

step1 首先在arcgis将图像裁剪

step2

from osgeo import gdal
import numpy as np

# 读写遥感数据
class Grid:
    # 读取图像文件
    def read_img(self, filename):
        data = gdal.Open(filename)  # 打开文件
        im_width = data.RasterXSize  # 读取图像行数
        im_height = data.RasterYSize  # 读取图像列数

        im_geotrans = data.GetGeoTransform()
        # 仿射矩阵，左上角像素的大地坐标和像素分辨率。
        # 共有六个参数，分表代表左上角x坐标；东西方向上图像的分辨率；如果北边朝上，地图的旋转角度，0表示图像的行与x轴平行；左上角y坐标；
        # 如果北边朝上，地图的旋转角度，0表示图像的列与y轴平行；南北方向上地图的分辨率。
        im_proj = data.GetProjection()  # 地图投影信息
        im_data = data.ReadAsArray(0, 0, im_width, im_height)  # 此处读取整张图像
        # ReadAsArray(<xoff>, <yoff>, <xsize>, <ysize>)
        # 读出从(xoff,yoff)开始，大小为(xsize,ysize)的矩阵。
        del data  # 释放内存，如果不释放，在arcgis，envi中打开该图像时会显示文件被占用

        return im_proj, im_geotrans, im_data

    def write_img(self, filename, im_proj, im_geotrans, im_data):
        # filename-创建的新影像
        # im_geotrans,im_proj该影像的参数，im_data，被写的影像
        # 写文件，以写成tiff为例
        # gdal数据类型包括
        # gdal.GDT_Byte,
        # gdal .GDT_UInt16, gdal.GDT_Int16, gdal.GDT_UInt32, gdal.GDT_Int32,
        # gdal.GDT_Float32, gdal.GDT_Float64

        # 判断栅格数据的类型
        if 'int8' in im_data.dtype.name:
            datatype = gdal.GDT_Byte
        elif 'int16' in im_data.dtype.name:
            datatype = gdal.GDT_UInt16
        else:
            datatype = gdal.GDT_Float32
        if len(im_data.shape) == 3:  # len(im_data.shape)表示矩阵的维数
            im_bands, im_height, im_width = im_data.shape  # （维数，行数，列数）
        else:
            im_bands, (im_height, im_width) = 1, im_data.shape  # 一维矩阵

        #         创建文件
        driver = gdal.GetDriverByName('GTiff')  # 数据类型必须有，因为要计算需要多大内存空间
        data = driver.Create(filename, im_width, im_height, im_bands, datatype)
        data.SetGeoTransform(im_geotrans)  # 写入仿射变换参数
        data.SetProjection(im_proj)  # 写入投影
        if im_bands == 1:
            data.GetRasterBand(1).WriteArray(im_data)  # 写入数组数据
        else:
            for i in range(im_bands):
                data.GetRasterBand(i + 1).WriteArray(im_data[i])
        del data


a = Grid()
#读取物候数据
im_proj, im_geotrans,VNP = a.read_img(r'VNP22Q2_2013_2019_clip.tif')


# 原始数据频数分布检查
# 提取LOS EOS SOS POS波段
# LOS---波段8
# EOS---波段9
# SOS---波段10
# POS---波段11

# 初始化4个空数组来存储结果
LOS = []
EOS = []
SOS = []
POS = []

# 切片数组
for i in range(7):  # 对应2013-2019总共7年
    year = 2013 + i
    start_idx = i * 38
    end_idx = start_idx + 38
    sub_array = VNP[start_idx:end_idx].astype(np.float64)  # 提取子数组并减去天数差

    # 从子数组中提取第8、9、10、11个数组并添加到对应的结果数组中
    LOS.append(sub_array[7])
    EOS.append(sub_array[8])
    SOS.append(sub_array[9])
    POS.append(sub_array[10])

# 将结果列表转换为numpy数组
LOS = np.array(LOS)
EOS = np.array(EOS)
SOS = np.array(SOS)
POS = np.array(POS)

print(LOS.shape, EOS.shape, SOS.shape, POS.shape)

#%%
#去除背景值（我的背景值是65535.0，替换为你的即可）
LOS[LOS == 65535.0] = np.nan
EOS[EOS == 65535.0] = np.nan
SOS[SOS == 65535.0] = np.nan
POS[POS == 65535.0] = np.nan
#去除无值（我在下载图像时设置质量为good quality,且云着实多，因此区域内许多像素无值）
LOS[LOS == 0.0] = np.nan
EOS[EOS == 0.0] = np.nan
SOS[SOS == 0.0] = np.nan
POS[POS == 0.0] = np.nan

#%%
#去除偏差
#计算偏差值
offset_list = []
for i in range(7): 
    year = 2013 + i
    print(year)
    offset = 366*(year-2000)
    offset_list.append(offset)

#对EOS SOS POS逐个逐年去除偏差
for i in range(7): 
    EOS[i,:,:] = EOS[i,:,:] - offset_list[i]
    SOS[i, :, :] = SOS[i, :, :] - offset_list[i]
    POS[i, :, :] = POS[i, :, :] - offset_list[i]
    
#%%导出
a.write_img(r'VNP_EOS.tif', im_proj, im_geotrans, EOS)
a.write_img(r'VNP_SOS.tif', im_proj, im_geotrans, SOS)
a.write_img(r'VNP_POS.tif', im_proj, im_geotrans, POS)

撒花完结~