将hdf文件读入R并将其转换为geoTIFF栅格 [英] Reading hdf files into R and converting them to geoTIFF rasters

问题描述

我正在尝试将MODIS 17数据文件读入R，对其进行处理(裁剪等)，然后将其另存为geoTIFF.数据文件采用.hdf格式，似乎没有一种简单的方法可以将它们读入R.

I'm trying to read MODIS 17 data files into R, manipulate them (cropping etc.) and then save them as geoTIFF's. The data files come in .hdf format and there doesn't seem to be an easy way to read them into R.

与其他主题相比，这里没有太多建议，而且大多数建议已有数年历史.其中一些建议还建议使用其他程序，但是我想坚持使用R.

Compared to other topics there isn't a lot of advice out there and most of it is several years old. Some of it also advises using additional programmes but I want to stick with just using R.

人们使用什么包处理R中的.hdf文件?

What package/s do people use for dealing with .hdf files in R?

推荐答案

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

Ok, so my MODIS hdf files were hdf4 rather than hdf5 format. It was surprisingly difficult to discover this, MODIS don't mention it on their website but there are a few hints in various blogs and stack exchange posts. In the end I had to download HDFView to find out for sure.

R不执行hdf4文件，并且几乎所有软件包(如rgdal)仅支持hdf5文件.有几篇关于从源代码下载驱动程序和编译rgdal的文章，但这些看起来都相当复杂，这些文章是针对MAC或Unix的，而我正在使用Windows.

R doesn't do hdf4 files and pretty much all the packages (like rgdal) only support hdf5 files. There are a few posts about downloading drivers and compiling rgdal from source but it all seemed rather complicated and the posts were for MAC or Unix and I'm using Windows.

基本上，gdalUtils软件包中的gdal_translate是想要在R中使用hdf4文件的任何人的节省之选.它可以将hdf4文件转换为geoTIFF，而无需将其读入R.这意味着您不能在全部例如通过裁剪它们，因此值得获得最小的切片(通过诸如Reverb之类的MODIS数据获取MODIS数据)以最小化计算时间.

Basically gdal_translate from the gdalUtils package is the saving grace for anyone who wants to use hdf4 files in R. It converts hdf4 files into geoTIFFs without reading them into R. This means that you can't manipulate them at all e.g. by cropping them, so its worth getting the smallest tiles you can (for MODIS data through something like Reverb) to minimise computing time.

这里是代码示例:

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和此的帮助！

Now you can read your .tif files into R using, for example, raster from the raster package and work as normal. I've checked the resulting files against a few I converted manually using QGIS and they match so I'm confident the code is doing what I think it is. Thanks to Loïc Dutrieux and this for the help!

这篇关于将hdf文件读入R并将其转换为geoTIFF栅格的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！