Modis数据下载及后处理

最近需要用大量的modis数据，这里记录一下最近踩的坑吧

下载modis的地址：

https://ladsweb.modaps.eosdis.nasa.gov/

点击上面的find data就可以愉快找数据了。我一开始只知道上面这个网址，后面还会给别的选择。

首先我想要MOD16A2和MOD43A3的全球数据，都是500m分辨率的，所以数据量还是很大的，如果直接下载这两个原始数据，发现hdf文件中有很多别的波段，比如MOD43A3竟然有28个波段，主要针对不同的电磁波波段，但是我真的需要的只是其中很少的3个数据，然后就发现在提交订单的时候其实是可以更改的，另外还可以直接选择投影的方式，也可以选择输出tif

由于hdf文件拼接需要用到modis官方推荐的工具HEG（之前用过MRT不过已经忘完了），我觉得好难用，于是就选择了tif输出加地理坐标，自己后期进行拼接。

在下载完所有数据后，我用arcpy进行了拼接，这里我拼接了三张图，每张都有12个月份，代码如下：

from __future__ import print_function
import sys
arcpy_path = [r'C:\Python27\ArcGIS10.3\Lib\site-packages',
              r'C:\Program Files (x86)\ArcGIS\Desktop10.3\arcpy',
              r'C:\Program Files (x86)\ArcGIS\Desktop10.3\bin',
              r'C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcToolbox\Scripts']
sys.path.extend(arcpy_path)
import arcpy
from arcpy.sa import *
import numpy as np
import os

raster = ['ET','LE','ET_QC']
path_work = 'H:/modis/tif/MOD16A2/'
arcpy.env.workspace = (path_work)
arcpy.env.overwriteOutput = True
imon=0
iraster = 0
for imon in range(12):
    for iraster in range(3):
        pathin_tif = 'H:/modis/tif/MOD16A2/{:0>2d}/'.format(imon+1)
        #---------------拼接--------------------
        file0 = os.listdir(pathin_tif)
        mosaic_rasters=""
        for f in file0:
            f0 = f.split('-')[-1]
            if f0==raster[iraster]+'_500m.tif':
                mosaic_rasters=mosaic_rasters+pathin_tif+f+";"
        arcpy.MosaicToNewRaster_management(mosaic_rasters, "Mosaic2New","{}_{:0>2d}.tif".format(raster[iraster],imon+1), "",\
                                   "16_BIT_SIGNED","" , "1", "MINIMUM","FIRST")

最后一步MosaicToNewRaster_management里选择MINIMUM是因为之前的投影是双曲余弦（modis的奇怪投影），转成地理坐标后就会有很多filled value，这些值很大相比于我们需要的值，所有用最小值填充才是对的。

经过了上面一通操作，其实工作量还是蛮大的，我发现原来上面这里可以直接在另一个USGS的网站内部完成qaq，下面就是这个方便的网址：

https://lpdaacsvc.cr.usgs.gov/appeears/task/area

这里只需自己新建一个mask就好了，mask可以用arcgis进行操作，新建一个polygon，然后edit，随便在地图上点4个点画成一个矩形，然后在旁边修改经纬度就好了

之后需要把shp文件打包成zip：(.zip including .shp, .dbf, .prj, and .shx files)

这里也可以直接选择地理坐标输出

大概处理全球的3张图花了半小时，速度还是很可以的。

最后再用gdal重采样就好了，我对比了两个方法基本上结果是一致的，只是有大概一个象元的差距，是可以理解的。

另外如果是8天的数据，但是分析的时候想要逐月分析，可以在选择时间的时候选择single date，输入15号的日子，系统就会自动筛选出覆盖了月中的8天，当然全把数据下载下来进行加权也可以，只是比较麻烦，数据量还更大了。