获取MOD16产品的实际蒸散发数据

月蒸散发数据下载链接：http://www.ntsg.umt.edu/project/modis/mod16.php

利用python中的pyModis模块对数据进行处理

以下为处理单个MOD16数据的代码

import pymodis
import gdal
import numpy as np

## download the monthly evapotranspiration data from http://www.ntsg.umt.edu/project/modis/mod16.php


inputfile = r'MOD16A2.A2007M06.h25v05.105.2013121062506.hdf'
output = r'Q:\ET'

## dataset transformation and reprojection
# hdfname ---> name of input data
# prefix ----> prefix for output data
# subset ----> which subset is read. for MOD16, 1--->ET, 2---->
# res ----> output resolution
# outformat -----> output format, it is possible to use all the support GDAL format,
#                  however, I meet error when using other GDAL format except GTiff
# epsg ----> the EPSG code for the preojection of output file
# wkt ----> the WKT string for the preojection of output file (Available: 'NODATUM', 'NAD27',
#           'NAD83', 'WGS66', 'WGS72', 'WGS84', [default = WGS84])
# resampl ----> the resampling method to use. methods:['NEAREST_NEIGHBOR', 'BICUBIC', 'CUBIC_CONVOLUTION',
#                'NONE',---gdal methods: 'AVERAGE', 'BILINEAR', 'CUBIC', 'CUBIC_SPINE', 'LANCZOS', 'MODE',
#                 'NEAREST_NEIGHBOR']
aa = pymodis.convertmodis_gdal(hdfname = inputfile, prefix = 'ET', subset = '1', res = 1000,
                               outformat = 'GTiff', epsg = '4326', wkt = 'WGS84',
                               resampl = 'NEAREST_NEIGHBOR' )

## run the transformation and reprojection
aa.run()

## read the GTIFF data
dataset = gdal.Open(r'Q:\ET_ET_1km.tif')
# get the x and y size
xsize = dataset.RasterXSize
ysize = dataset.RasterYSize

# get the information of orginal coordinate (upper left point) and range of each grid
# coordinate[0] -----> longitude
# coordinate[1] -----> longitude range of one grid
# coordinate[3] -----> latitute
# coordinate[5] -----> latitute range of one grid
coordinate = dataset.GetGeoTransform()

# read the ET data
data = dataset.ReadAsArray()

# create the grid mesh using numpy
xrange = range(0, xsize)
yrange = range(0, ysize)
x, y = np.meshgrid(xrange, yrange)
lon = coordinate[0] + x*coordinate[1] + y*coordinate[2]
lat = coordinate[3] + x*coordinate[4] + y*coordinate[5]