Viele Raster in R effizient lesen und neu klassifizieren?

8

Ich wurde beauftragt, eine Eignungsanalyse der Wellenbedingungen im Golf von Mexiko zu erstellen. Ich habe ungefähr zweitausend Rasterdateien mit jeweils ca. 8 MB (2438 Spalten, 1749 Zeilen, 1 km Zellengröße). Der Parameter in den Rasterdateien ist die Wellenperiode, und ich möchte alle Raster so umklassifizieren, dass if 4<= wave period <=9 then make cell = 1, else cell = 0. Dann alle Raster zu einem endgültigen Raster zusammenfassen und durch die Gesamtzahl der Raster dividieren, um a zu erhalten Gesamtprozentsatz geeigneter Beobachtungen und Exportergebnisse in ein ESRI-kompatibles Format ... möglicherweise etwas, das Floats bei Bedarf unterstützen kann. Ich habe weder mit Python noch mit R viel gearbeitet, aber nach der Online-Suche scheint es sinnvoll, diesen Prozess in einer dieser Sprachen durchzuführen. Ich habe mir bisher einen Code in R ausgedacht, bin aber verwirrt darüber, wie das funktioniert.

library(rgdal)
raster_data <- list.files(path=getwd())    #promt user for dir containing raster files
num_files <- length(raster_data)
for (i in raster_data) {                   #read in rasters
   my_data <- readGDAL(raster_data[i])

An diesem Punkt bin ich verwirrt, ob ich auch die Daten innerhalb dieser Schleife neu klassifizieren und summieren soll oder nicht. Meine Vermutung wäre ja, da ich sonst denke, dass mir möglicherweise der Speicher auf meinem Computer ausgehen würde, aber ich bin mir einfach nicht sicher. Ich bin mir auch nicht sicher, wie ich die Daten neu klassifizieren soll.

Bei Online-Recherchen denke ich, ich würde verwenden reclass(my_data, c(-Inf,3,0, 4,9,1, 10,Inf,0)), aber sieht das richtig aus?

Und zum Summieren würde ich das verwenden sum(stack(my_data))und irgendwie ausgeben. Auch ... wenn dies effizienter ausgeführt oder in Python geschrieben werden könnte, wäre ich auch dafür offen. Ich bin wirklich ein Anfänger, wenn es um Programmierung geht.

Nigel
quelle
Verwenden Sie einfach Python-GDAL. Es wird in Ihrem Fall viel effizienter sein.
SS_Rebelious
Danke, Rebelious. Nur neugierig zu wissen, warum Python-GDAL in dieser Situation effizienter ist? Wäre es auch möglich, einige der Schritte im Code zu sehen, die dazu erforderlich wären? Es ist verwirrend herauszufinden, wie der Code so geschrieben wird, dass er die Daten einliest, verarbeitet, aus dem Speicher entfernt und Fahren Sie dann mit dem nächsten Raster fort.
Nigel
Ich kann Ihnen nicht genau sagen, warum, aber die allgemeine Ursache ist, dass R für andere Zwecke entwickelt wurde und bekanntermaßen bei Zyklen langsam arbeitet. Wenn zum Beispiel kein Code bereitgestellt wird, werde ich einen in ungefähr 10 Stunden mit Ihnen teilen, wenn ich Zugriff auf den Computer erhalte, auf dem das entsprechende Skript gespeichert ist.
SS_Rebelious

Antworten:

8

Dies ist eine kurze Möglichkeit, dies in R zu tun --- hier ohne Zwischendateien:

library(raster)
raster_data <- list.files(path=getwd())    #promt user for dir containing raster files
s <- stack(raster_data)
f <- function(x) { rowSums(x >= 4 & x <= 9) }
x <- calc(s, f, progress='text', filename='output.tif')
Robert Hijmans
quelle
1
+1 Dies ist gut für kleine Probleme, aber lassen Sie uns diesbezüglich rechnen: 2438 Spalten mal 1749 Zeilen mal 8 Bytes / Wert mal zweitausend Gitter = 63,6 GB, die alle Rzum Erstellen im RAM bleiben müssen s. (Wahrscheinlich wird doppelt so viel RAM benötigt, da es swahrscheinlich nicht ersetzt wird raster_data.) Ich hoffe, Sie haben viel RAM! Ihre Lösung kann praktikabel gemacht werden, indem Sie die 2000 Gitter in kleinere Gruppen aufteilen, die Berechnung für jede Gruppe durchführen und diese Berechnungen dann kombinieren.
whuber
2
@whuber: 's' ist ein kleines Objekt, nur ein paar Zeiger auf die Dateien. Die Calc-Funktion lädt wie andere Funktionen im Rasterpaket nicht alle Werte in den Speicher. es wird sie in Stücken verarbeiten. Das heißt, die Aufteilung in Gruppen, wie Sie vorschlagen, erfolgt automatisch hinter den Kulissen. Die Blockgröße kann für die über rasterOptions () verfügbare RAM-Größe optimiert werden.
Robert Hijmans
1
Danke, dass du das erklärt hast! Ich hatte dies ohne Überprüfung angenommen stackund calcarbeitete wie die meisten anderen RFunktionen, indem ich zuerst alle Daten in den RAM lud.
whuber
+1 Ich mag die Prägnanz von R gegenüber dem bereitgestellten Python-Beispiel ...
SlowLearner
5

Wie ich in den Kommentaren bemerkt habe, sollten Sie R im Allgemeinen aufgrund von Leistungsproblemen in bestimmten Aspekten für nicht statistische Zwecke vermeiden (das Arbeiten mit Zyklen ist ein Beispiel). Hier ist ein Codebeispiel für Sie in Pyhton (dank dieses Artikels ) zur Neuklassifizierung einer einzelnen Datei mit einem einzelnen Band. Sie können es problemlos für die Stapelverarbeitung ändern, wenn Sie wissen, wie Sie alle Dateien aus dem Verzeichnis abrufen . Beachten Sie, dass Raster als Arrays dargestellt werden. Sie können daher Array-Methoden verwenden, um die Leistung gegebenenfalls zu verbessern. Informationen zum Arbeiten mit Arrays in Python finden Sie in der Numpy-Dokumentation .

UPD: Der Code, den ich ursprünglich gepostet habe, war eine abgeschnittene Version eines benutzerdefinierten Filters, der pro Pixelverarbeitung benötigt wurde. Bei dieser Frage erhöht die Verwendung von Numpy jedoch die Leistung (siehe Code).

from osgeo import gdal
import sys
import numpy

gdalData = gdal.Open("path_to_file")
if gdalData is None:
  sys.exit("ERROR: can't open raster")

#print "Driver short name", gdalData.GetDriver().ShortName
#print "Driver long name", gdalData.GetDriver().LongName
#print "Raster size", gdalData.RasterXSize, "x", gdalData.RasterYSize
#print "Number of bands", gdalData.RasterCount
#print "Projection", gdalData.GetProjection()
#print "Geo transform", gdalData.GetGeoTransform()


raster = gdalData.ReadAsArray()
xsize = gdalData.RasterXSize
ysize = gdalData.RasterYSize

#print xsize, 'x', ysize

## per pixel processing
#for col in range(xsize):
#  for row in range(ysize):
#    # if pixel value is 16 - change it to 7
#    if raster[row, col] == 16:
#      raster[row, col] = 7

#reclassify raster values equal 16 to 7 using Numpy
temp = numpy.equal(raster, 16)
numpy.putmask(raster, temp, 7)

# write results to file (but at first check if we are able to write this format)
format = "GTiff"
driver = gdal.GetDriverByName(format)
metadata = driver.GetMetadata()
if metadata.has_key(gdal.DCAP_CREATE) and metadata[gdal.DCAP_CREATE] == "YES":
  pass
else:
  print "Driver %s does not support Create() method." % format
  sys.exit(1)
if metadata.has_key(gdal.DCAP_CREATECOPY) and metadata[gdal.DCAP_CREATECOPY] == "YES":
  pass
else:
  print "Driver %s does not support CreateCopy() method." % format
  sys.exit(1)

# we already have the raster with exact parameters that we need
# so we use CreateCopy() method instead of Create() to save our time
outData = driver.CreateCopy("path_to_new_file", gdalData, 0)
outData.GetRasterBand(1).WriteArray(raster)
SS_Rebelious
quelle
4
Das "R ist langsam, wenn Zyklen (Schleifen) ausgeführt werden" wird häufig als Grund für die Vermeidung von R missbraucht. Ja, wenn Sie die Zellen eines Rasters in R durchlaufen, ist es langsam, aber das Rasterpaket funktioniert auf ganze Raster gleichzeitig und hat viel C-Code und läuft daher mit C-Geschwindigkeit. Für ein Raster dieser Größe würde der größte Teil der Arbeit bei C-Geschwindigkeiten ausgeführt, der Schleifenaufwand wäre unbedeutend.
Spacedman
@Spacedman, ja, 'Raster' ist ein nützliches Paket (und ich mag es), aber ich war nie zufrieden mit seiner Leistung, selbst wenn keine Schleifen beteiligt waren.
SS_Rebelious
2
Okay, vergleichen Sie die Zeit, die in R benötigt wird, mit der Zeit, die in Python benötigt wird. Können Sie nicht ein ganzes numpy-Array bearbeiten, anstatt eine Schleife zu erstellen?
Spacedman
@Spacedman, ich habe gerade die Antwort aktualisiert.
SS_Rebelious
Vielen Dank an euch beide. Ich werde versuchen, sowohl mit dem von Ihnen bereitgestellten Python-Code als auch mit etwas R zu basteln und zu sehen, was ich erreichen kann. Ich werde mit Ergebnissen oder Problemen aktualisieren.
Nigel
2
  1. Verwenden Sie nicht readGDAL. Es liest in ein Spatial * -Objekt ein, was möglicherweise keine gute Idee ist.

  2. Verwenden Sie das rasterPaket. Es kann GDAL-Dinge in Raster-Objekte einlesen. Das ist eine gute Sache. r = raster("/path/to/rasterfile.tif")werde es in lesen r.

  3. Ihre Klassifizierung ist dann t = r > 4 & r <= 9

  4. Die große Frage ist, ob diese in neue Rasterdateien ausgegeben und dann der Zusammenfassungsschritt in einer anderen Schleife ausgeführt werden sollen. Wenn Sie den Speicher haben, würde ich sie in neue Dateien schreiben, nur weil Sie erneut starten müssen, wenn Ihre Schleife fehlschlägt (weil eine dieser 2000 Dateien Junk ist). Verwenden writeRasterSie diese Option, um Raster mit Schwellenwerten zu erstellen, wenn Sie sich dazu entschließen.

  5. Ihre Schleife ist dann nur so etwas wie

diese.

count = raster(0,size of your raster)
for(i in 1:number of rasters){
  r = raster(path to binary raster file 'i')
  count = count + r
}
return(count)

Die Speicherverwaltung von R könnte Sie hier treffen - wenn Sie dies tun, erstellt count=count+rR möglicherweise eine neue Kopie von count. Das sind also möglicherweise 2000 Kopien davon - aber die Speicherbereinigung wird aktiviert und aufgeräumt, und hier war die Schleife von R früher sehr schlecht.

In Bezug auf das Timing dauert auf meinem 4-jährigen PC der Schwellenwert op bei einem Raster dieser Größe von Zufallszahlen zwischen null und zwanzig etwa 1,5 Sekunden. Times 2000 = Sie rechnen. Erstellen Sie wie immer einen kleinen Testsatz für die Entwicklung des Codes, starten Sie ihn dann mit Ihren realen Daten und trinken Sie einen Kaffee.

Spacedman
quelle
Ich bin gespannt, was Sie unter "Schwelle op" verstehen. Auf meinem System (mit R2.15.0) dauert das Lesen eines 1,6-Megapixel-Rasters (im nativen ESRI-Gleitkommaformat) 0,19 Sekunden und das Ausführen der fünf Schleifenoperationen - zwei Vergleiche, eine Konjunktion, eine Addition und ein Speicher - dauert eine weitere 0,09 Sekunden für 0,28 Sekunden pro Iteration. Wenn stattdessen die Schleife ausgeführt wird, indem jeweils 100 Gitter in einem Array zwischengespeichert werden und rowSumsdie Summen verwendet werden, steigt natürlich die RAM-Auslastung. Ebenso offensichtlich gibt es jedoch keine Speicherbereinigung. Das Timing sinkt nur auf 0,26 Sekunden pro Iteration.
whuber