2016-04-21 16 views
4

Ich versuche MODIS 17-Dateien in R zu lesen, manipuliere sie (Beschneiden usw.) und speichere sie als geoTIFFs. Die Datendateien kommen in .hdf Format und es scheint nicht eine einfache Möglichkeit zu sein, sie in R zu lesen.hdf-Dateien in R lesen und in geoTIFF-Raster konvertieren

Im Vergleich zu anderen Themen gibt es nicht viele Ratschläge da draußen und das meiste davon ist mehrere Jahre alt . Einige davon rät auch, zusätzliche Programme zu verwenden, aber ich möchte nur mit R bleiben.

Welche Paket/s verwenden Leute für den Umgang mit .hdf Dateien in R?

+0

Haben Sie [* diesen Artikel *] (http://www.r-bloggers.com/working-with-hdf-files-in-r-example-pathfinder-sst-data/) auf R-Bloggern überprüft ? Es bietet ein funktionierendes Beispiel, ob die angebotene Lösung * einfach * ist, es wäre eine Frage des individuellen Urteils. In Bezug auf Ihre zweite Frage, bietet die Methode ['writeRater'] (http://www.inside-r.org/packages/cran/raster/docs/writeRaster) im [**' Raster' **] (http://www.inside-r.org/packages/cran/raster) Paket von Nutzen für Sie? – Konrad

+0

@Konrad Ja, ich stieß auf diesen Artikel und verbrachte eine Weile damit, es zum Laufen zu bringen. Problem ist, dass ich den folgenden Fehler bekomme, wenn ich versuche, die h5ls-Funktion für meine Daten zu verwenden: Fehler in H5Fopen (Datei, "H5F_ACC_RDONLY"): HDF5. Dateizugriffsfähigkeit. Datei kann nicht geöffnet werden. Ich benutze das Raster-Paket sehr, es ist das Einlesen und Bearbeiten von .hdf-Dateien Ich interessiere mich für die Meinungen von Leuten, da das ist, was mich Probleme verursacht. – James

+0

Bitte denken Sie daran, Ihr Problem reproduzierbar zu machen in Bezug auf die in [dieser Diskussion] geäußerten Vorschläge (http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example?lq= 1); Es wird einfacher sein, dir zu helfen. – Konrad

Antwort

10

Ok, also meine MODIS-hdf-Dateien waren hdf4 statt hdf5-Format. Es war überraschend schwierig, dies zu entdecken, MODIS erwähnen es nicht auf ihrer Website, aber es gibt ein paar Hinweise in verschiedenen Blogs und Stack Exchange Posts. Am Ende musste ich HDFView herunterladen, um es herauszufinden.

R macht keine hdf4-Dateien und so ziemlich alle Pakete (wie rgdal) unterstützen nur hdf5-Dateien. Es gibt ein paar Posts zum Herunterladen von Treibern und zum kompilieren von rgdal aus der Quelle, aber alles schien ziemlich kompliziert und die Posts waren für MAC oder Unix und ich benutze Windows.

gdal_translate aus dem gdalUtils Paket ist die Rettung für alle, die hdf4-Dateien in R verwenden möchten. Es konvertiert hdf4-Dateien in geoTIFFs, ohne sie in R zu lesen. Dies bedeutet, dass Sie sie überhaupt nicht manipulieren können. Indem man sie beschneidet, lohnt es sich, die kleinsten Kacheln zu bekommen (für MODIS-Daten durch etwas wie Reverb), um die Rechenzeit zu minimieren.

Hier und Beispiel für den Code: Verwendung zum Beispiel raster aus dem Raster-Paket und arbeitet als normales

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.2015113.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]) 
} 

Jetzt können Sie Ihre TIF-Dateien in R lesen. Ich habe die resultierenden Dateien mit ein paar manuell über QGIS konvertierten Dateien verglichen und sie stimmen überein, so dass ich sicher bin, dass der Code tut, was ich denke. Danke an Loïc Dutrieux und this für die Hilfe!

+0

Ich habe das nicht für eine große ~ 290 funktionieren lassen .hdf-Dateien. Hat jemand sonst Probleme? – Evan

+0

@James. Ich habe die Antwort unten gepostet, da diese Konvertierung im Vergleich zu ArcGIS zu anderen Ergebnissen führt. Irgendwelche Tipps? – MIH

1

Dieses Skript war sehr nützlich und ich habe es geschafft, einen Stapel von 36 Dateien zu konvertieren. Mein Problem ist jedoch, dass die Konvertierung nicht korrekt erscheint. Wenn ich es mit ArcGIS 'NetCDF-Raster-Layer-Werkzeug erstellen' erhalte ich unterschiedliche Ergebnisse + Ich kann die Zahlen von Kelvin nach C konvertieren mit einer einfachen Formel: RasterValue * 0.02 - 273.15. Mit den Ergebnissen der R-Konvertierung bekomme ich nach der Konvertierung nicht die richtigen Ergebnisse, was zu der Annahme führt, dass die ArcGIS-Konvertierung gut ist und die R-Konvertierung einen Fehler zurückgibt.

library(gdalUtils) 
library(raster) 

setwd("D:/Users/szynisze/Google Drive/Gates Project/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]) 
} 
-2

Folgendes funktionierte für mich. Es ist ein kurzes Programm und nimmt nur den Namen des Eingabeordners auf. Stellen Sie sicher, dass Sie wissen, welche Subdaten Sie möchten. Ich interessierte mich für Subdaten 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")) 
} 
+0

Es funktionierte für mich. Es ist ein kurzes Programm und nimmt nur den Namen des Eingangsordners auf. Stellen Sie sicher, dass Sie wissen, welche Subdaten Sie möchten. Ich war an Subdaten interessiert 1. –