logo标签列表

将hdf文件读入R并转换为geoTIFF栅格图像。

rgisgdalgeotiffhdf
19

我正在尝试将MODIS 17数据文件读入R中，对它们进行操作（裁剪等），然后将它们保存为geoTIFF文件。这些数据文件以.hdf格式提供，似乎没有一种简单的方法将它们读入到R中。

与其他主题相比，相关建议并不多，而且大部分都是几年前的。其中一些还建议使用其他程序，但我想坚持只使用R。

人们在R中处理.hdf文件时会使用哪些软件包？

你有看过R-bloggers上的这篇文章吗？它提供了一个工作示例，是否提供的解决方案容易，这将是个人判断的问题。关于你的第二个问题，writeRater方法在raster包中是否对你有用？ - Konrad
@Konrad 是的，我看到了那篇文章，并花了一段时间尝试让它工作。问题是，当我尝试在我的数据上使用h5ls函数时，我不断收到以下错误：Error in H5Fopen(file, "H5F_ACC_RDONLY") : HDF5. File accessability. Unable to open file.我经常使用raster包，我对人们对读取和操作.hdf文件的看法很感兴趣，因为这正是导致我问题的原因。 - James
请考虑按照此讨论中提出的建议，使您的问题可重现；这将更容易帮助您解决问题。 - Konrad
2
如果未来读者遇到此问题，R中的“MODISTools”软件包提供了很好的功能，可将MODIS数据导入R中（但不是处理HDF4文件的通用工具）。 - Jacob Socolar
5个回答
27

好的，我的MODIS hdf文件是hdf4格式，而不是hdf5格式。很难发现这一点，因为MODIS没有在他们的网站上提到它，但在一些博客和Stack Exchange帖子中有一些提示。最后，我不得不下载HDFView来确定。

R不能处理hdf4文件，几乎所有包（如rgdal）只支持hdf5文件。有一些关于下载驱动程序和从源代码编译rgdal的帖子，但这似乎相当复杂，并且这些帖子都是针对MAC或Unix，而我使用的是Windows。

基本上，gdalUtils包中的gdal_translate对于想在R中使用hdf4文件的任何人来说都是救命稻草。它将hdf4文件转换为geoTIFFs而不会将其读入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从raster包中读取您的.tif文件，并像往常一样处理。我已经将生成的文件与我手动使用QGIS转换的一些文件进行了比较，它们是匹配的，因此我对代码执行的操作非常有信心。感谢Loïc Dutrieux和这里的帮助！

我还没有成功处理大约290个.hdf文件。有其他人遇到麻烦吗？ - Evan
@James。我在下面发布了答案，因为这种转换与ArcGIS产生了不同的结果。有什么建议吗？ - MIH
嘿，感谢您的回答。我已经尝试使用这段代码，在250m和500m的MODIS产品上运行良好，但在1km产品上返回15个波段的TIF文件而不是36个。此外，get_subdatasets函数似乎没有返回文件中的所有36个数据波段。有什么想法吗？ - Ollie Perkins
Linux的file命令可以轻松告诉您HDF版本。 - Rodrigo
12

现在你可以使用terra软件包处理HDF文件

获取子数据集

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

这可以像这样提取为SpatRasters 

s[1]

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

r <- rast("file.hdf")
2
Robert Hijmans。您是众人之王。 - Ollie Perkins
2
2
这个脚本非常有用，我已经成功地使用它来转换36个文件。然而，我的问题是转换似乎不正确。当我使用ArcGIS的'Make NetCDF Raster Layer'工具时，我得到了不同的结果，并且可以使用简单的公式将数字从开尔文转换为摄氏度：RasterValue * 0.02 - 273.15。使用R转换后，我得不到正确的结果，这让我相信ArcGIS的转换是正确的，而R的转换返回了错误。
library(gdalUtils)
library(raster)

setwd("D:/Data/Climate/MODIS")

# Get a list of sds names
sds <- get_subdatasets('MOD11C3.A2009001.006.2016006051904.hdf')
# Isolate the name of the first sds
name <- sds[1]
filename <- 'Rasterinr.tif'

gdal_translate(sds[1], dst_dataset = filename)
# Load the Geotiff created into R
r <- raster(filename)

# Identify files to read:
rlist=list.files(getwd(), pattern="hdf$", full.names=FALSE)


# Substract last 5 digits from MODIS filename for use in a new .img filename
substrRight <- function(x, n){
        substr(x, nchar(x)-n+1, nchar(x))
}

filenames0 <- substrRight(rlist,9)
# Suffixes for MODIS files for identyfication:
filenamessuffix <- substr(filenames0,1,5)

listofnewnames <- c("2009.01.MODIS_","2009.02.MODIS_","2009.03.MODIS_","2009.04.MODIS_","2009.05.MODIS_",
                    "2009.06.MODIS_","2009.07.MODIS_","2009.08.MODIS_","2009.09.MODIS_","2009.10.MODIS_",
                    "2009.11.MODIS_","2009.12.MODIS_",
                    "2010.01.MODIS_","2010.02.MODIS_","2010.03.MODIS_","2010.04.MODIS_","2010.05.MODIS_",
                    "2010.06.MODIS_","2010.07.MODIS_","2010.08.MODIS_","2010.09.MODIS_","2010.10.MODIS_",
                    "2010.11.MODIS_","2010.12.MODIS_",
                    "2011.01.MODIS_","2011.02.MODIS_","2011.03.MODIS_","2011.04.MODIS_","2011.05.MODIS_",
                    "2011.06.MODIS_","2011.07.MODIS_","2011.08.MODIS_","2011.09.MODIS_","2011.10.MODIS_",
                    "2011.11.MODIS_","2011.12.MODIS_")

# Final new names for converted files:
newnames <- vector()
for (i in 1:length(listofnewnames)) {
        newnames[i] <- paste0(listofnewnames[i],filenamessuffix[i],".img")
}

# Loop converting files to raster from NetCDF
for (i in 1:length(rlist)) {
        sds <- get_subdatasets(rlist[i])
        gdal_translate(sds[1], dst_dataset = newnames[i])
}
你解决这个问题了吗？我在转换MODIS数据时也遇到类似的问题；似乎丢失了比例因子，但数值仍然不正确。 - SpatialProgramming
2
以下是我所使用的方法。这是一个简短的程序，只需输入文件夹名称即可。确保您知道要使用哪个子数据。我对子数据1感兴趣。
library(raster)
library(gdalUtils)

inpath <- "E:/aster200102/ast_200102"

setwd(inpath)

filenames <- list.files(,pattern=".hdf$",full.names = FALSE)

for (filename in filenames)
{
     sds <- get_subdatasets(filename)
     gdal_translate(sds[1], dst_dataset=paste0(substr(filename, 1, nchar(filename)-4) ,".tif"))
}
它对我起作用了。这是一个简短的程序，只需要输入文件夹名称。确保你知道你想要哪个子数据。我对子数据1感兴趣。 - Rajasweta Datta
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接