【MODIS】数据下载与处理

目录：

目录

方法一：Data Pool 直接下载

方法二：NASA Earthdata Search

方法三：python下载

处理数据：

情况一：瓦片

情况二：全球月均求年均

---

以MOD11C3 v.061数据为例，这是月均地表温度数据

在MODIS官网的Data栏进行搜索   https://lpdaac.usgs.gov/

进入数据页面之后又ACCESS DATA，里面有几个选项：

方法一：Data Pool 直接下载

进去会给出数据文件list，要哪个点哪个，数据量大会比较麻烦

方法二：NASA Earthdata Search

需要注册一个账号，右上角login可以注册

进去后继续选择：

然后选择日期，数据可能是瓦片，也可能是每个月的数据

然后继续选择Download

然后会出现以下界面，一定要等到进度条100%。之后普遍有两种选择。一是Download Files底下直接点击链接下载，量多麻烦。第二就是有Linux系统的，从Download Scrip下载脚本。

# 赋予权限
chomd 777 XX.sh 
# 运行
./XX.sh
# 随后会出现输入用户名和密码
# 如果不想每次都输入，可以在脚本里做如下修改
prompt_credentials() {
    echo "Enter your Earthdata Login or other provider supplied credentials"
    #read -p "Username (用户名): " username
    #username=${username:-用户名}
    #read -s -p "Password: " password
    username="用户名"
    password="密码"
    echo "machine urs.earthdata.nasa.gov login $username password $password" >> $netrc
    echo
}

修改脚本

然后就可以进行下载了。

方法三：python下载

待探索

处理数据：

分两种情况处理1 下载的全球瓦片 2 全球的月均数据求年均（python+arcgis处理）

情况一：瓦片

1 hdf转tif

import numpy as np
import arcpy
arcpy.CheckOutExtension("Spatial")
arcpy.gp.overwriteOutput=1

yr = np.arange(2001,2023,1)
for year in yr:
    year = str(year)
    inDir=r'E:/MODIS/MOD11C3 v.061/' + year +'/' #H:\HDFFILES\ALLHDFS
    outDir=r'E:/MODIS/MOD11C3 v.061/' + year+'/'
    
    arcpy.env.workspace=inDir
    rasters=arcpy.ListRasters("*","hdf")
    print(year)
    for raster in rasters:
        print(raster)
        outName=outDir+'\\'+raster[:-4]+'.tif'
        arcpy.ExtractSubDataset_management(raster,outName)
        print(outName)

2 合并tif

from arcpy import env
for year in yr:
    print(year)
    year = str(year)
    # set environment setting
    env.workspace = 'E:/LAND_DATA/MCD12QV061/'+ year+'/'
    # set output path
    output = "E:/LAND_DATA/MCD12QV061/merge/"
    
    # set local variables
    rasters = arcpy.ListRasters('*', 'tif')
    
    mosaic_rasters = ""
    for raster in rasters:
       mosaic_rasters = mosaic_rasters + raster + ';'
    #print(mosaic_rasters)
    
    arcpy.CheckOutExtension('Spatial')
    arcpy.MosaicToNewRaster_management(mosaic_rasters, output, year+'.tif', "", "16_BIT_UNSIGNED", "", '1', "", "")

3 批量转投影

先用arcgis转一个出来作为模板，之后的都可以照着转换

import arcpy
inDir=r'E:/LAND_DATA/MCD12QV061/merge'
outDir=r'E:/LAND_DATA/MCD12QV061/merge/reproject'
originReferenceFile=r'E:/LAND_DATA/MCD12QV061/merge/2001.tif'
referenceFile=r'E:/LAND_DATA/MCD12QV061/merge/reproject/2001_Reproject.tif'

arcpy.CheckOutExtension("spatial")
arcpy.gp.overwriteOutput=1

arcpy.env.workspace=inDir
rasters=arcpy.ListRasters("*","tif")

for raster in rasters:
    if ('2001' in raster)  : #continue
        print(raster)
        outName=outDir+'\\'+raster[:-4]+'_Reproject1.tif'
        print(outName)
        arcpy.ProjectRaster_management(raster, outName,referenceFile, "#",\
                                       '#','#','#',originReferenceFile)

情况二：全球月均求年均

1 读取hdf

2 求平均

# 处理LST数据需要乘0.02来转变成正常值，之前属于16位无符号整型。
from pyhdf.SD import SD
import os
files = [f for f in os.listdir('E:/MODIS/MOD11C3 v.061/2013/') if f[-3:] == 'hdf']
Whole_Day = 0
for file in files :
    df = SD('E:/MODIS/MOD11C3 v.061/2013/'+file,)
    #print(df.info())
    datasets_dic = df.datasets()
    a = (df.select('LST_Day_CMG').get()*0.02-272.15  + df.select('LST_Night_CMG').get()*0.02-272.15 )/2
    Whole_Day += a

Whole_Day = Whole_Day/len(files)

3 输出tif

lonmin,lonmax,latmin,latmax = (-180,180,-90,90) # 根据实际更改
nx,ny = (7200,3600) # 经度、纬度
f = Whole_Day
# 创建一个tif ,写入数据   
tiffpath = r'E:/MODIS/MOD11C3 v.061/2013/2013.tif'
from osgeo import gdal,osr
driver = gdal.GetDriverByName("GTiff")
New_YG_dataset = driver.Create(tiffpath, nx, ny, 1, gdal.GDT_Float32)
lon_res = np.round((lonmax-lonmin)/nx,2)
lat_res = np.round((latmax-latmin)/ny,2)
geotransform = (lonmin,lon_res, 0, latmax, 0, -lat_res)
New_YG_dataset.SetGeoTransform(geotransform)
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326) # 定义输出的坐标系为"WGS 84"，AUTHORITY["EPSG","4326"]
New_YG_dataset.SetProjection(srs.ExportToWkt()) # 给新建图层赋予投影信息
New_YG_dataset.GetRasterBand(1).WriteArray(f)
New_YG_dataset.FlushCache() # 将数据写入硬盘
del New_YG_dataset

 注意：LST数据不能直接用，要成0.02才能变成开尔文，-272.15变成摄氏度