将hdf文件读入R并将其转换为geoTIFF栅格

Question

我正在尝试将MODIS 17数据文件读入R,操纵它们(裁剪等),然后将它们保存为geoTIFF.数据文件.hdf格式化,似乎没有一种简单的方法将它们读入R.

与其他主题相比,没有很多建议,其中大部分都是几年之久.其中一些还建议使用其他程序,但我想坚持使用R.

人们用什么包来处理.hdfR中的文件？

Answer 1

好的,所以我的MODIS hdf文件是hdf4而不是hdf5格式.很难发现这一点,MODIS没有在他们的网站上提到它,但在各种博客和堆栈交换帖中有一些提示.最后我不得不下载HDFView以确定.

R不执行hdf4文件,几乎所有的包(如rgdal)都只支持hdf5文件.有一些关于下载驱动程序和从源代码编译rgdal的帖子,但它们看起来相当复杂,帖子是针对MAC或Unix而且我使用的是Windows.

对于任何想要在R中使用hdf4文件的人来说,基本上gdal_translate从gdalUtils包中获得了优势.它将hdf4文件转换为geoTIFF而不将其读入R.这意味着你无法操纵它们,例如通过裁剪它们,所以它的价值获得最小的瓷砖(通过类似Reverb的MODIS数据),以最大限度地缩短计算时间.

这是代码的例子:

library(gdalUtils)

# Provides detailed data on hdf4 files but takes ages

gdalinfo("MOD17A3H.A2000001.h21v09.006.2015141183401.hdf")

# Tells me what subdatasets are within my hdf4 MODIS files and makes them into a list

sds <- get_subdatasets("MOD17A3H.A2000001.h21v09.006.2015141183401.hdf")
sds

[1] "HDF4_EOS:EOS_GRID:MOD17A3H.A2000001.h21v09.006.2015141183401.hdf:MOD_Grid_MOD17A3H:Npp_500m"   
[2] "HDF4_EOS:EOS_GRID:MOD17A3H.A2000001.h21v09.006.2015141183401.hdf:MOD_Grid_MOD17A3H:Npp_QC_500m"

# I'm only interested in the first subdataset and I can use gdal_translate to convert it to a .tif

gdal_translate(sds[1], dst_dataset = "NPP2000.tif")

# Load and plot the new .tif

rast <- raster("NPP2000.tif")
plot(rast)

# If you have lots of files then you can make a loop to do all this for you

files <- dir(pattern = ".hdf")
files

 [1] "MOD17A3H.A2000001.h21v09.006.2015141183401.hdf" "MOD17A3H.A2001001.h21v09.006.2015148124025.hdf"
 [3] "MOD17A3H.A2002001.h21v09.006.2015153182349.hdf" "MOD17A3H.A2003001.h21v09.006.2015166203852.hdf"
 [5] "MOD17A3H.A2004001.h21v09.006.2015099031743.hdf" "MOD17A3H.A2005001.h21v09.006.2015113012334.hdf"
 [7] "MOD17A3H.A2006001.h21v09.006.2015125163852.hdf" "MOD17A3H.A2007001.h21v09.006.2015169164508.hdf"
 [9] "MOD17A3H.A2008001.h21v09.006.2015186104744.hdf" "MOD17A3H.A2009001.h21v09.006.2015198113503.hdf"
[11] "MOD17A3H.A2010001.h21v09.006.2015216071137.hdf" "MOD17A3H.A2011001.h21v09.006.2015230092603.hdf"
[13] "MOD17A3H.A2012001.h21v09.006.2015254070417.hdf" "MOD17A3H.A2013001.h21v09.006.2015272075433.hdf"
[15] "MOD17A3H.A2014001.h21v09.006.2015295062210.hdf"

filename <- substr(files,11,14)
filename <- paste0("NPP", filename, ".tif")
filename

[1] "NPP2000.tif" "NPP2001.tif" "NPP2002.tif" "NPP2003.tif" "NPP2004.tif" "NPP2005.tif" "NPP2006.tif" "NPP2007.tif" "NPP2008.tif"
[10] "NPP2009.tif" "NPP2010.tif" "NPP2011.tif" "NPP2012.tif" "NPP2013.tif" "NPP2014.tif"

i <- 1

for (i in 1:15){
  sds <- get_subdatasets(files[i])
  gdal_translate(sds[1], dst_dataset = filename[i])
}

现在,您可以使用例如raster栅格包中的.tif文件读取到R中并正常工作.我已经使用QGIS手动转换的几个文件检查了结果文件,它们匹配,所以我相信代码正在做我认为的事情.感谢LoïcDutrieux和这个帮助!

Answer 2

现在您可以使用terra包含 HDF 文件的软件包

要么获取子数据集

 library(terra)
 s <- sds("file.hdf")
 s

可以像这样提取为 SpatRasters

s[1]

或者像这样创建所有子数据集的 SpatRaster

r <- rast("file.hdf")