十亩洞人分享 http://blog.sciencenet.cn/u/shimudongren 全球变化

博文

利用R中的MODIS包获取和处理相关遥感数据

已有 1745 次阅读 2020-9-2 10:32 |系统分类:科研笔记

一、题外话

MODIS数据在资源环境领域应用广泛(产品介绍网: https://modis.gsfc.nasa.gov/data/dataprod/ )。目前国内获取MODIS数据比较常规的途径:一是去官方网站下载。二是遥感界新秀GEE(但是由于网络限制,需要买VPN翻墙实现)。(据说航天宏图也整了一个类似的,http://engine.piesat.cn/,刚上市可以测试)。


目前NASA将其公开的遥感资源都集中在EARTHDATA上面了(https://earthdata.nasa.gov/),注册后可以进行下载,不再是以前的FTP下载方式(ftp方式在2018年4月20日就就退休啦)。


改成了HTTPS协议的下载方式。下图就是EARTHDATA的用户注册和登陆界面,自己按照要求注册,就可以在里面搜索数据和下载啦。

官方网站也提供了wget、Linux shel、Perl和Python下载数据的脚本案例。https://ladsweb.modaps.eosdis.nasa.gov/tools-and-services/data-download-scripts/

二、MODIS包

对于经常用R的童鞋,也不想整什么高大上的GEE,MODIS包就是很大的福音/利器。

MODIS包也在github上,上面有常见问题模块,可以向开发者提问和提供bug.

https://github.com/MatMatt/MODIS

亲测问过一个问题,1天左右回复,比较及时。

主要的使用步骤如下:

1,在EARTHDATA注册用户。

到官方网站注册https://earthdata.nasa.gov/,比如你的用户名和密码分别是:用户名:james;密码:lileihanmeimei123

2,安装MODIS包

install.packages("MODIS")

安装成功了。

3,下载数据流程

3.1 配置用户名和密码

第一个是earthdata的账户配置,在MODIS包中 EarthdataLogin函数来配置,假如你的earthdata账户和密码是下面的例子:

EarthdataLogin(usr = 'james', pwd = 'lileihanmeimei123')

在R中运行上面的代码后,会在电脑上生成一个不可见的密码文件。

3.2 配置下载和处理数据目录

lap = "D:/R_MODIS/MODIS_ARC" # '设置MODIS下载到本地的位置' 

odp = file.path(lap, "PROCESSED") # '设置MODIS处理数据的位置' MODISoptions(localArcPath = lap, outDirPath = odp) # 配置文件

配置好后,运行 MODISoptions(),控制台出现下面的信息

STORAGE: 

_______________ 

localArcPath : D:/R_MODIS/MODIS_ARC

outDirPath : D:/R_MODIS/MODIS_ARC/PROCESSED 

DOWNLOAD: _______________ 

MODISserverOrder : LPDAAC, LAADS 

dlmethod : auto 

stubbornness : high 

wait : 0.5 

quiet : TRUE 

PROCESSING: _______________ 

GDAL : 3.0.4 

MRT : Enabled 

pixelSize : asIn 

outProj : asIn 

resamplingType : NN 

dataFormat : GTiff 

cellchunk : 1

3.3 下载数据

(1)先看看可以下载哪些产品

运行 getProduct() 出现了详细的产品列表,目前用166种MODIS产品可以通过这个包下载。详细的介绍可以看官方网站:https://modis.gsfc.nasa.gov/data/dataprod/

(2)利用getHdf函数进行下载

getHdf函数是下载数据的主要函数,简单说明如下:

getHdf(

product, #产品名称,就是上面getProduct()的结果

HdfName,# hdf文件的名称,可以直接提供

begin = NULL, #开始的时间,一般提供day of year

end = NULL, #结束的时间,一般提供day of year

tileH, # MDOIS的行幅号

tileV, # MDOIS的列幅号

extent, #下载的范围

collection = NULL,

checkIntegrity = TRUE,

forceDownload = TRUE,

...

)


  • 直接使用hdf名字下载

a <- getHdf(HdfName = c("MYD11A1.A2009001.h18v04.006.2015363221538.hdf",

"MYD11A1.A2009009.h18v04.006.2015364055036.hdf",

"MYD11A1.A2009017.h18v04.006.2015364115403.hdf"))


  • 提供开始和结束时间等信息下载

c <- getHdf(product = "MOD09GA", begin = "2020165", end = "2020165",

tileH = 27, tileV = 5)


  • 提供范围下载

Austria <- extent(9.2, 17.47, 46.12, 49.3)

b2 <- getHdf(product = "MOD13A2", begin = "2020001", end = "2020031", extent = Austria)

  • 查看下载的文件

每次下载结束后,到配置好的文件夹中去查看,比如我的:

  • 查看行列幅号

先在控制台getTile(),右边点击需要查看的“图块”,然后DONE。