博主分享了如何利用Python和GDAL库将MODIS MOD16A2GFv061蒸散发数据从8天尺度转换为月尺度的过程。代码详细解释了从文件读取、处理跨月影像权重分配到计算月累计值的步骤，并提供了代码实现。处理过程中考虑了数据的有效值范围和时间权重，确保了转换的科学性。
最近处理一批长时间序列MODIS蒸散发数据，官方数据为8-days，想转成月尺度的数据，想用Python+GDAL批处理，本来觉得应该也不算多复杂的一个工作，想着偷懒在网上搜一下看看有没有现成的代码可以参考，结果中英文网页查完后居然都没有找到合适的代码，感觉不太能理解，现有的几篇博客介绍的处理方法跟我自己设想的差距有点远，都是直接找每个月的大概4景，然后做个平均值或者最大值合成，我个人感觉这样处理不是特别科学，因此零零碎碎花了大概2~3天时间自己写了一个代码，我自己写代码水平一般，很多地方可能写的也不是很简洁，有问题的话欢迎大家批评指正。
首先，简单说有下背景及思路，1. 我用的蒸散发数据是MOD16A2GF v061（LP DAAC - MOD16A2GF）版本；2. 原数据单景为每8天的和，一年46景，并且每年都是从1日开始（需要注意的是影像日期为一个8天的起始日期，如2003001代表的是2003年1月1日~2003年1月8日的蒸散发和，2003025代表的是2003年1月25日~2003年2月1日的蒸散发和，我写代码的时候这个地方还弄错了一次），另外我处理的数据是2003年1月~2014年12月，代码中可以根据需要修改；3. 原始的MODIS ET数据是HDF格式的，需要先转为tif格式，然后重投影、裁剪、重采样到500m分辨率，这个预处理工作比较简单，我自己先用的官方的HEG工具转的格式，然后用gdal.Warp处理的，因此我的输入数据是tif格式的；4. 我的思路是首先根据影像DOY转换为日期，把同一年相同月份的影像存放一起，其次如果有跨月的影像，找出该月份首尾可能的跨月影像，最后根据跨月影像在各个月所占天数分配权重到各月，然后求和得到月累计值（mm/month），如果需要月平均值（mm/day）直接用月累计值除以各月的天数就可以得到。具体代码如下，代码中我加的注释比较详细，就不再继续详细解释了。其实最开始我本来还想把数据的质量波段QC文件也加进来，把质量差的值剔除掉，但是最近事情太多了，也就没有去做，代码中是根据官方文档把有效值范围（-32767~32700）外的值直接设为背景值了，并且栅格值乘以尺度因子0.1换算成了实际的物理值。
写代码的时候，简单画了个图帮助理解，并附上官方文档关于该产品描述的截图。
# -*- coding: utf-8 -*-
import os
import re
import sys
import glob
import datetime
import numpy as np
from os import path
from osgeo import gdal
# Input Path(8-days)
et_path = r'D:\MHH\15_SWAT\3_Data\4_ET\MOD16A2GF_v061_GeoTIFF_Clip'
# Output Path(Monthly)
et_mon_path = r'D:\MHH\15_SWAT\3_Data\4_ET\MOD16A2GF_v061_GeoTIFF_Clip_Monthly'
# 栅格数据访问
class RasterTiff:
gdal.AllRegister()
# 读栅格文件
def read_img(self, filename):
dataset = gdal.Open(filename) # 打开文件
im_width = dataset.RasterXSize # 栅格矩阵的列数
im_height = dataset.RasterYSize # 栅格矩阵的行数
im_bands = dataset.RasterCount # 波段数
im_geotrans = dataset.GetGeoTransform() # 仿射矩阵，左上角像素的大地坐标和像素分辨率
im_proj = dataset.GetProjection() # 地图投影信息，字符串表示
im_data = dataset.ReadAsArray(0, 0, im_width, im_height) # 将数据写成数组，对应栅格矩阵
type(im_data), im_data.shape
del dataset
return im_height, im_width, im_bands, im_geotrans, im_proj, im_data
# 写栅格文件
def write_img(self, filename, im_proj, im_geotrans, im_data):
# gdal数据类型包括：gdal.GDT_Byte, gdal.GDT_UInt16, 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:
im_bands, im_height, im_width = im_data.shape
else:
im_bands, (im_height, im_width) = 1, im_data.shape
# 创建文件
driver = gdal.GetDriverByName('GTiff') # 数据类型必须有
dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)
dataset.SetGeoTransform(im_geotrans) # 写入仿射变换参数
dataset.SetProjection(im_proj) # 写入投影
if im_bands == 1:
dataset.GetRasterBand(1).WriteArray(im_data) # 写入数组数据
else:
for i in range(im_bands):
dataset.GetRasterBand(i + 1).WriteArray(im_data[i])
del dataset
# 某月第一天/最后一天/天数
def mon_statistics(year_para, month_para):
# 当月第一天
first_day_mon = datetime.date(year_para, month_para, day=1)
first_day_next_mon = datetime.date(year_para + 1, month=1, day=1) if month_para == 12 \
else datetime.date(year_para, month_para + 1, day=1)
# 当月最后一天
last_day_mon = first_day_next_mon - datetime.timedelta(days=1)
# 当月天数
days_mon = last_day_mon.day
return first_day_mon, last_day_mon, days_mon
# ET数据按月聚合
et_list = glob.glob(os.path.join(et_path, '*.tif'))
for y_idx in range(2003, 2015):
for m_idx in range(1, 13):
month_list = []
print('Date:', str(y_idx) + '-' + str(m_idx))
for et_name in et_list:
# 获取文件名中日期信息
name = re.findall(r'MOD16A2GF.A(\d{7}).h28v06', path.basename(et_name))[0]
date = datetime.datetime.strptime(name, '%Y%j')
year, month, day = date.year, date.month, date.day
if year == y_idx:
if month == m_idx:
# 同一年相同月份的影像存放一起并插入各景影像数据在各月份中所占权重
# 根据所占月份时间长短确定权重（暂时设为1）
month_list.append([et_name, 1.0])
# 检查开头日期（是否从1日开始，注：每年第一个月必然从1日开始，无需考虑当月跨月问题）
name_mon_first = re.findall(r'MOD16A2GF.A(\d{7}).h28v06', path.basename(month_list[0][0]))[0]
date_mon_first = datetime.datetime.strptime(name_mon_first, '%Y%j')
if date_mon_first.day != 1: # 当月第一景数据不是从1日开始，则加入上一景数据
# 往前推8天找到上一景影像
date_last_mon = date_mon_first - datetime.timedelta(days=8)
name_last_mon = datetime.datetime.strftime(date_last_mon, '%Y%j')
# 在影像列表中查找并插入上个月最后一景影像（权重暂时分配为默认值1）
for et_name_sup in et_list:
name_sup = re.findall(r'MOD16A2GF.A(\d{7}).h28v06', path.basename(et_name_sup))[0]
if name_sup == name_last_mon:
month_list.insert(0, [et_name_sup, 1.0])
## 为每月首尾可能跨月影像分配权重
name_mon_list_first = re.findall(r'MOD16A2GF.A(\d{7}).h28v06', path.basename(month_list[0][0]))[0]
name_mon_list_last = re.findall(r'MOD16A2GF.A(\d{7}).h28v06', path.basename(month_list[-1][0]))[0]
date_mon_list_first = datetime.datetime.strptime(name_mon_list_first, '%Y%j')
date_mon_list_last = datetime.datetime.strptime(name_mon_list_last, '%Y%j')
## 每月首尾景影像年月日
year_mon_list_first, month_mon_list_first, day_mon_list_first = date_mon_list_first.year, \
date_mon_list_first.month, \
date_mon_list_first.day
year_mon_list_last, month_mon_list_last, day_mon_list_last = date_mon_list_last.year, \
date_mon_list_last.month, \
date_mon_list_last.day
## 判断首景影像是否属于上月跨月影像（注：判断其是否从1日开始）
if day_mon_list_first != 1:
# 属于跨月影像，则获取上个月第一天/最后一天/天数
first_day_mon, last_day_mon, days_mon = mon_statistics(year_mon_list_first, month_mon_list_first)
last_day_mon_datetime = datetime.datetime.strptime(str(last_day_mon), '%Y-%m-%d')
# 求首景影像（跨月）在上个月所占天数
delta_last_mon_last = (last_day_mon_datetime - date_mon_list_first).days + 1
# 求首景影像（跨月）在当月所占天数并据此分配影像权重
delta_this_mon_first = 8 - delta_last_mon_last
month_list[0][1] = delta_this_mon_first / 8
## 当年当月第一天/最后一天/天数
first_day_this_mon, last_day_this_mon, days_this_mon = mon_statistics(y_idx, m_idx)
## 判断尾景数据是否属于下月跨月数据（注：判断尾景数据起始日期与当月最后一日是否相等，注：12月尾景不考虑跨月问题）
if (day_mon_list_last + 7) != days_this_mon and m_idx != 12:
last_day_this_mon_datetime = datetime.datetime.strptime(str(last_day_this_mon), '%Y-%m-%d')
# 获取尾景影像（跨月）当月所占天数并据此分配影像权重
delta_this_mon_last = (last_day_this_mon_datetime - date_mon_list_last).days + 1
month_list[-1][1] = delta_this_mon_last / 8
# ET月尺度转换
month_list_et = []
## 读取栅格数据
raster = RasterTiff()
for et_name_mon in month_list:
# 读数据
et_rows, et_columns, et_bands, et_geotrans, et_proj, et_data = raster.read_img(et_name_mon[0])
# 数组转换为浮点型，便于后续剔除无效值
et_data = et_data.astype(np.float32)
# 获取影像栅格矩阵并添加到新列表
et_name_mon.extend([et_rows, et_columns, et_bands, et_geotrans, et_proj, et_data])
# 列表参数顺序：影像名，权重，et_rows, et_columns, et_bands, et_geotrans, et_proj, et_data
month_list_et.append(et_name_mon)
## 根据各景影像权重进行月尺度转换
et_data_mon_aggregation = np.zeros(shape=(month_list_et[0][2], month_list_et[0][3]), dtype=np.float32)
for et_data_mon in month_list_et:
# <-32767以及>32700设为空值，不参与计算
et_data_mon[-1][np.where((et_data_mon[-1] < -32767) | (et_data_mon[-1] > 32700))] = np.nan
et_data_mon_aggregation += et_data_mon[1] * et_data_mon[-1]
# 月尺度ET（月累计值，单位mm/month）
# 如果需要换算成月平均值，需要除以各月的天数，单位：mm/day
## 尺度因子
Scale_Factor = 0.1
et_data_mon_aggregation = et_data_mon_aggregation * Scale_Factor
raster.write_img(et_mon_path + '\\' + 'MOD16A2GF_ET_Mon_' + str(y_idx) + '_' + str(m_idx) + '.tif',
month_list_et[0][6], month_list_et[0][5], et_data_mon_aggregation)
print('ET 8-days aggregation to monthly done!')
https://download.csdn.net/download/u011534341/82018871
4607
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
