MODIS数据批量下载与处理-pymodis的使用
关注公众号,分享GIS知识、ArcGIS教程、SCI论文与科研日常等
MODIS数据,下载量小的话,可以从NASA官网 根据产品、时间、位置进行筛选自己要的,之后下载。当下载量比较大时…肯定是编程来的快了…之前一直在想怎么编程,后来接触python后,得知pymodis工具…一个函数搞定。后续拼接处理也有相应的函数,可以调用MRT(官网上好像已经下架,替代产品还未了解)软件和GDAL库
PS.这种方法找不到需要下载的产品,可以试一下这个方法
关于pymodis
pymodis是基于python语言的开源modis数据处理库。pymodis 提供按照用户提供的时间范围的modis产品的批量下载函数、对modis瓦片数据进行镶嵌以及投影、将hdf格式转化为其他格式,并提取数据质量信息。
具有如下特点:
1. 用于下载大量modis hdf/xml 文件以及解析xml文档以获取有关HDF文件的信息
2. 使用MRT或GDAL将modis文件转为geotiff格式
3. 使用MRT或GDAL将modis瓦片数据拼接、重投影
4. 用于下载大量modis hdf/xml 文件以及解析xml文档以获取有关HDF文件的信息
5. 创建融合后影像的xml格式元数据
6. 从modis字节编码质量评估图层提取具体信息
包括:
pymodis脚本工具:
modis_download.py
modis_download_from_list.py
modis_parse.py
modis_multiparse.py
modis_mosaic.py
modis_convert.py
modis_quality.py
pymodis库:
downmodis module
parsemodis module
convertmodis module
convertmodis_gdal module
qualitymodis module
optparse_required module
optparse_gui module
脚本工具使用
warning:无论用脚本工具还是模块编程,下载数据是从NASA的官方网站上下载,所以在使用下载前要提前在网站上注册,记住用户名和密码
脚本可以在终端调用,通过下面命令方式:
modis_download.py [options] destination_folder
其中[options]为可选命令,由下面参数组成,参数都设置了默认值,一般要根据自己需求修改。
-h --help shows the help message and exit
-u --url http/ftp server url [default=https://e4ftl01.cr.usgs.gov]
-I --input insert user and password from standard input
-P --password password to connect
-U --username username to connect
-t --tiles string of tiles separated by comma
[default=none] for all tiles
-s --source directory on the http/ftp server
[default=MOLT]
-p --product product name as on the http/ftp server
[default=MOD11A1.005]
-D --delta delta of day starting from first day [default=10]
-f --firstday the day to start download, if you want change
data you have to use this format YYYY-MM-DD
([default=none] is for today)
-e --endday the day to finish download, if you want change
data you have to use this format YYYY-MM-DD
([default=none] use delta option)
-x useful for debugging the download
[default=False]
-j download also the jpeg files [default=False]
-O download only one day, it sets delta=1 [default=False]
-A download all days, useful for initial download of a
product. It overwrites the 'firstday' and 'endday'
options [default=False]
-r remove files with size same to zero from
'destination_folder' [default=False]
如下载过去15天的Aqua卫星LST数据
modis_download.py -I -r -s MOLA -p MYD11A1.005 -t h18v03,h18v04 -D 15 lst_aqua/
还是要彻底了解一下各个参数的输入及意义,有的看的一头雾水,不了解变量很难变通。解析一下前面的代码,其中modis_download.py就是脚本命令,-I -r -s MOLA -p MYD11A1.005 -t h18v03,h18v04 -D 15这部分都是可选输入参数部分,lst_aqua为下载数据的存储文件夹
其中,可选输入参数部分,每个都是变量标识后跟相应参数,为空则是默认变量值,其中 -s 是网站上的目录,为 MOLA;-p是产品名称,为MYD11A1.005;-t是要下载产品的行列号,多个用逗号隔开,在这里下载h18v03,h18v04,即18行3列和18行4列的两个数据瓦片;-D表示从下载日期开始,往历史搜寻的天数,在这个时间段的数据即下载下来,这里下载从今天开始过去15天内的数据。其中下载开始日期/结束日期都可以分别通过 -f 和 -e设定,开始日期默认为今天。其中,刚开始比较迷惑的是如何确定网站下的产品所在的那一级目录、产品名称,这里可以参考这个网站。
MYD11A1.005是在MOLA目录下的。库的函数参数的确定与此相同。总结来说,明确工具或函数的各个变量含义才能灵活使用。
再比如modis_convert.py脚本,他的运行命令为:
modis_convert.py [options] hdf_file
他的可选参数输入部分如下:此处调用了MRT软件或GDAL库,分别给出了针对这两种方法的输入参数
-h, --help show this help message and exit
Required options:
-s SUBSET, --subset=SUBSET
(Required) a subset of product's layers. The string
should be similar to: ( 1 0 )
-o OUTPUT_FILE, --output=OUTPUT_FILE
(Required) the prefix of output file
-g RESOLUTION, --grain=RESOLUTION
the spatial resolution of output file
-r RESAMPLING_TYPE, --resampl=RESAMPLING_TYPE
the method of resampling. -- mrt methods:
'NEAREST_NEIGHBOR', 'BICUBIC', 'CUBIC_CONVOLUTION',
'NONE' -- gdal methods: 'AVERAGE', 'BILINEAR',
'CUBIC', 'CUBIC_SPLINE', 'LANCZOS', 'MODE',
'NEAREST_NEIGHBOR' [default=NEAREST_NEIGHBOR]
Options for GDAL:
-f OUTPUT_FORMAT, --output-format=OUTPUT_FORMAT
output format supported by GDAL [default=GTiff]
-e EPSG, --epsg=EPSG
EPSG code for the output
-w WKT, --wkt_file=WKT
file or string containing projection definition in WKT
format
-v, --vrt Read from a GDAL VRT file.
--formats print supported GDAL formats
Options for MRT:
-m MRT_PATH, --mrt=MRT_PATH
the path to MRT software
-d DATUM, --datum=DATUM
the code of datum. Available: 'NODATUM', 'NAD27',
'NAD83', 'WGS66', 'WGS72', 'WGS84' [default=WGS84]
-t PROJECTION_SYSTEM, --proj_type=PROJECTION_SYSTEM
the output projection system. Available: 'AEA', 'GEO',
'HAM', 'IGH', 'ISIN', 'LA', 'LCC', 'MOL', 'PS', 'SIN',
'TM', 'UTM', 'MERCAT' [default=GEO]
-p PROJECTION_PARAMETERS, --proj_parameters=PROJECTION_PARAMETERS
a list of projection parameters, for more info check
the 'Appendix C' of MODIS reprojection tool user's
manual https://lpdaac.usgs.gov/content/download/4831/2
2895/file/mrt41_usermanual_032811.pdf [default=( 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 )]
-u UTM_ZONE, --utm=UTM_ZONE
the UTM zone if projection system is UTM
the UTM zone if projection system is UTM
例如,将modis数据的前三个图层转化为500米分辨率UTM32投影:
modis_convert.py -s "( 1 1 1 )" -o OUTPUT_FILE -g 500 -e 32632 FILE
结构同上,每个变量标识代表的含义也都有。
downmodis module批量下载modis数据
downModis类初始化参数如下:
def __init__(self, destinationFolder, password=None, user=None,
url="https://e4ftl01.cr.usgs.gov", tiles=None, path="MOLT",
product="MOD11A1.005", today=None, enddate=None, delta=10,
jpg=False, debug=False, timeout=30, checkgdal=True):
各参数含义:
destinationFolder: 下载数据存储的目标文件夹
password: 用户注册密码
user: 用户名
url: 下载地址,这个参数不变
tiles: 下载条带号,多条用逗号分隔,或数组形式
path: 在下载地址与产品目录之间的那层目录,参照网站上内容
product: 产品名称
today: 开始下载的日期,yyyy-mm-dd形式
enddate: 下载结束日期,形式同上,结束日期应该在开始下载日期之前,downmodis是在时间轴上反向下载数据。
delta: 下载天数,表示从下载日期开始向历史追溯的天数。enddate和delta设定一个,当都有参数时,以enddate为准。
jpg: 是否下载JPG数据
debug=False, timeout=30, checkgdal=True,不重要可默认。
python代码示例:
用pymodis库中的downmodis类中的相关函数可以实现快速下载,包括的变量及含义与运行脚本一致。因为变量比较多,最好将函数变量的变量名带着,易读且不用考虑变量顺序问题,以免造成困扰。
convertmodis_gdal module批量拼接投影等预处理
借助于GDAL的modis预处理模块,主要包括createMosaicGDAL 和 convertModisGDAL 等类。
createMosaicGDAL 实现拼接;convertModisGDAL 将hdf中指定波段导出转换为TIFF,并可对TIFF影响重采样、投影变换等。
createMosaicGDAL 和 convertModisGDAL怎么同时用呢,还没搞懂,createMosaicGDAL 可以将不同条带号的拼接起来,但是不能进行投影重采样,当处理多个波段,拼接的结果,也是一个TIF文件多个波段集的形式。而convertModisGDAL虽然可以能够将不同波段的数据分别导出并进行投影重采样,但是不能实现拼接。两者的输入数据都是modis的hdf数据。可行方案:createMosaicGDAL拼接导出hdf格式,在用convertModisGDAL将拼接后的hdf进行接下来的投影重采样。但是如果用convertModisGDAL输出结果格式设置为hdf,计算速度超慢。因此,可以只用 convertModisGDAL将需要的波段(如NDVI)导出,同时进行了投影变换和重采样,虽然结果是一个个瓦片,没有实现拼接,但是可以之后再用arcpy或者GDAL等在进行TIF影响的融合/拼接操作,也是一个可行方案
convertModisGDAL类 的初始化方法:
def __init__(self, hdfname, prefix, subset, res, outformat="GTiff",
epsg=None, wkt=None, resampl='NEAREST_NEIGHBOR', vrt=False):
各参数含义:
hdfname: 输入的hdf文件(路径)
prefix: 保存数据的前缀。这里比较奇怪,相当于确定输入的路径和保存名称。“D:\User Files\caishuohao\MODIS_Data\MOD13A3_convert\MOD13A3_A2019152_h27v07”,如果只处理了NDVI波段,最后保存的文件就是"D:\User Files\caishuohao\MODIS_Data\MOD13A3_convert\MOD13A3_A2019121_h27v05_1 km monthly NDVI.tif"。前缀包含的是存储文件夹位置和一部分文件名。而1 km monthly NDVI.tif 则是处理时候自己生成的,如果选择多个波段也以这种形式保存。
subset: 处理波段,形如"(1 0)",0,1 代表是否处理,位序代表第几个波段。如"(1 0)" 其中1在第一位 表示第一个波段,值1表示处理;0在第二位 表示第二个波段 0代表不处理。
outformat: 输出格式,GDAL的保留格式
epsg: 投影或地理坐标系的EPSG码,epsg查询网址
wkt: 描述空间参照系统(转换)的文本标记语言,可以看什么是wkt。 也就是指定地理或投影坐标系的意思。epsg和wkt两者都是指定投影或地理坐标系,有其中之一即可,若两者都有,则用epsg,其实如果只指定epsg也是,程序也是通过epsg获取wkt。wkt比较复杂,这里有个简单的确定输入参数的方法,即指定地理投影坐标系的arcGIS的.prj文件即可。
resampl: 重采样方法
vrt:
python代码示例:
convertmodis = pymodis.convertmodis_gdal.convertModisGDAL(
hdfname=hdf_name, # 输入的hdf(路径)
prefix=_perfix, # 输入前缀:路径+半个文件名(后半个程序自己加)
subset="( 1 )", # 只处理第一个波段 NDVI
res=0.01, # 分辨率
outformat='GTiff', # 输入格式
epsg=None, # 我输入wgs84的epsg code:4326 有问题,还不知道原因
wkt=r"D:\User Files\caishuohao\MODIS_Data\HKHshp\HKHbianjie.prj", # 这里输入的是个投影文件,程序自己会解析
resampl='NEAREST_NEIGHBOR', # 重采样方法
vrt=False
)
其他方法:Python调用MRT批处理
可联系