Чтение hdf-файлов в R и преобразование их в растры geoTIFF

Я пытаюсь прочитать файлы данных MODIS 17 в R, манипулировать ими (обрезка и т. Д.), А затем сохранить их как geoTIFF. Файлы данных приходят в .hdf формат, и, кажется, нет простого способа прочитать их в R.

По сравнению с другими темами там не так много советов, и большинству из них уже несколько лет. Некоторые из них также советуют использовать дополнительные программы, но я хочу придерживаться только использования R.

С какими пакетами люди работают? .hdf файлы в R?

Источник

5 ответов

Решение

Итак, мои файлы hdf MODIS были hdf4, а не hdf5. Это было на удивление трудно обнаружить, MODIS не упоминают об этом на своем веб-сайте, но есть несколько советов в различных блогах и обмене сообщениями в стеке. В конце концов мне пришлось скачать HDFView, чтобы узнать наверняка.

R не делает файлы hdf4 и почти все пакеты (например, rgdal) поддерживает только файлы hdf5. Есть несколько постов о загрузке драйверов и компиляции rgdal из исходного кода, но все это казалось довольно сложным, и посты были для MAC или Unix, и я использую Windows.

В принципе gdal_translate от gdalUtils Пакет является спасительной льготой для всех, кто хочет использовать файлы hdf4 в R. Он преобразует файлы hdf4 в geoTIFF, не считывая их в R. Это означает, что вы не можете манипулировать ими вообще, например, обрезая их, поэтому стоит получить наименьшее плитки вы можете (для данных 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])
}

Теперь вы можете читать ваши файлы.tif в R, используя, например, raster из растрового пакета и работай как обычно. Я проверил полученные файлы по нескольким, которые я конвертировал вручную, используя QGIS, и они совпадают, поэтому я уверен, что код делает то, что я думаю. Спасибо Loïc Dutrieux и этому за помощь!

Источник

В наши дни вы можете использовать terra пакет с файлами HDF

Либо получить поднаборы данных

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

Это можно извлечь как SpatRasters, как это

      s[1]

Или создайте SpatRaster для всех подобных наборов данных

      r <- rast("file.hdf")
Источник

Следующее сработало для меня. Это короткая программа, которая берет только имя входной папки. Убедитесь, что вы знаете, какие субданные вы хотите. Я был заинтересован в данных суб 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"))
}
Источник

Этот скрипт был очень полезен, и мне удалось преобразовать пакет из 36 файлов, используя его. Однако моя проблема в том, что преобразование не выглядит правильным. Когда я делаю это, используя ArcGIS "Инструмент создания растрового слоя NetCDF", я получаю разные результаты + я могу конвертировать числа в К из Кельвина, используя простую формулу: 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])
}
Источник

Используйте инструментарий HEG, предоставленный NASA, чтобы преобразовать ваш hdf-файл в geotiff, а затем используйте любой пакет (например, "растровый"), чтобы прочитать файл. Я делаю то же самое для старых и новых файлов hdf.

Вот ссылка: https://newsroom.gsfc.nasa.gov/sdptoolkit/HEG/HEGHome.html

Посмотрите на продукты NASA, которые поддерживаются здесь: https://newsroom.gsfc.nasa.gov/sdptoolkit/HEG/HEGProductList.html

Надеюсь это поможет.

Источник
Другие вопросы по тегам