Wie kann ich mit R Daten in Form von lat, lon, value in eine Rasterdatei konvertieren?

40

Ich habe einen Datensatz mit Werten über ein km-Raster in den kontinentalen USA. Die Spalten lauten "Breite", "Länge" und "Beobachtung", z.

"lat"    "lon"     "yield"
 25.567  -120.347  3.6 
 25.832  -120.400  2.6
 26.097  -120.454  3.4
 26.363  -120.508  3.1
 26.630  -120.562  4.4

oder als R-Datenrahmen:

mydata <- structure(list(lat = c(25.567, 25.832, 26.097, 26.363, 26.63), 
lon = c(-120.347, -120.4, -120.454, -120.508, -120.562), 
yield = c(3.6, 2.6, 3.4, 3.1, 4.4)), .Names = c("lat", 
"lon", "yield"), class = "data.frame", row.names = c(NA, -5L))

(der vollständige Datensatz kann hier als csv heruntergeladen werden )

Die Daten werden von einem Erntemodell (vorgesehen) in einem Raster von 30 km x 30 km (von Miguez et al. 2012 ) ausgegeben .

Bildbeschreibung hier eingeben

Wie kann ich diese in eine Rasterdatei mit GIS-bezogenen Metadaten wie Kartenprojektion konvertieren?

Idealerweise ist die Datei eine Textdatei (ASCII?), Da ich möchte, dass sie plattform- und softwareunabhängig ist.

raster convert r Abe
quelle
Als CSV, das bereits ist eine „Textdatei“ in ASCII. Da überhaupt keine Projektion verwendet wird, müssen möglicherweise nur wenige relevante Metadaten hinzugefügt werden (Datum, meistens). Könnten Sie etwas genauer beschreiben, welche Art von Output Sie suchen und was Sie damit vorhaben?
Whuber
Ich möchte es jemandem so einfach wie möglich machen, die Daten mit einer Vielzahl von Kartierungssoftware (ArcGIS, Google Maps, Grass, R usw.) zu verwenden, um die Wiederverwendung zu erleichtern, z. B. indem keine zusätzlichen Konvertierungsschritte erforderlich sind. Basierend auf der Wikipedia-Seite von GIS-Dateiformaten schließe ich, dass 1) eine "Raster" -Datei Rownamen mit Breiten- und Spaltennamen von Längengrad wie ein Bild haben sollte und dass 2) Metadaten geografische Informationen enthalten sollten (Position einer Ecke, abgedeckter Bereich) nach Daten).
Abe
Dies ist eine der besten Referenzen, auf die ich bei R und GIS gestoßen bin. Vielen Dank! Kannst du bitte eine andere csv mit lat und long mit korrektem proj4string versehen? Ich werde das wirklich zu schätzen wissen.
@Nandini nicht sicher , was die richtige proj4string ist, vermute ich , lambert konform:proj +proj=lcc +lat_1=50.0 +lat_2=50.0 +units=km +lon_0=-145.5 +lat_0=1.0 . Ich bin mir nicht sicher, wonach Sie in Bezug auf eine andere CSV-Datei fragen - wie würde sie sich von der in der Frage verlinkten unterscheiden, oder würde sie von der akzeptierten Antwort erzeugt werden?
Abe
für mich geht nicht Ich weiß nicht, was ich "x" und "y" an "Koordinaten (Punkte) = ~ x + y" setzen soll

Antworten:

44

Mehrere Schritte erforderlich:

  1. Sie sagen, es ist ein reguläres 1-km-Raster, aber das bedeutet, dass die Lat-Long-Werte nicht regulär sind. Zuerst müssen Sie es in ein reguläres Gitterkoordinatensystem umwandeln, damit die X- und Y-Werte einen regelmäßigen Abstand haben.

    ein. Lesen Sie es als Datenrahmen in R ein, wobei die Spalten x, y, ergeben.

    pts = read.table("file.csv",......)

    b. Konvertieren Sie den Datenrahmen mit dem sp-Paket in einen SpatialPointsDataFrame.

    library(sp)
library(rgdal)
coordinates(pts)=~x+y

    c. Stellen Sie auf Ihr reguläres km-System um, indem Sie zuerst angeben, um welches CRS es sich handelt, und dann spTransform zum Ziel.

    proj4string(pts)=CRS("+init=epsg:4326") # set it to lat-long
pts = spTransform(pts,CRS("insert your proj4 string here"))

    d. Sagen Sie R, dass dies gerastert ist:

    gridded(pts) = TRUE

    An diesem Punkt erhalten Sie eine Fehlermeldung, wenn Ihre Koordinaten nicht auf einem schönen regelmäßigen Gitter liegen.

  2. Verwenden Sie nun das Raster-Paket, um in ein Raster zu konvertieren und dessen CRS festzulegen:

    r = raster(pts)
projection(r) = CRS("insert your proj4 string here")

  3. Jetzt schau mal:

    plot(r)

  4. Schreiben Sie es nun mit dem Raster-Paket als geoTIFF-Datei:

    writeRaster(r,"pts.tif")

Dieses geoTIFF sollte in allen gängigen GIS-Paketen lesbar sein. Das offensichtlich fehlende Teil hier ist der proj4-String, in den konvertiert werden muss: Dies wird wahrscheinlich eine Art UTM-Referenzsystem sein. Schwer zu sagen ohne weitere Daten ...

Raumfahrer
quelle
+1 Vielen Dank für die Gestaltung des Workflows. Beachten Sie, dass die Daten unter dem in der Frage angegebenen Link verfügbar sind: Schauen Sie rein. Sie werden leider feststellen, dass einige Ihrer Annahmen über sie falsch sind. (Insbesondere gejagt ich für jede Dokumentation über die Projektion verwendet , um das Raster zu erstellen , aber keine gefunden und es ist eine seltsame Projektion, wie Sie durch Auftragen der Punkte sehen können..)
whuber
Es ist sehr nahe daran, ein UTM-System zu sein, aber keines der Systeme, die ich ausprobiert habe, ist nahe genug an einem regulären Raster, damit R sie rastert. Ich bin halb versucht, die gesamte epsg-Datenbank von R zu durchlaufen ...
Spacedman
Das wäre eine echte Tour de Force, wenn Sie die Projektion auf diese Weise entdecken könnten! Der Schlüssel besteht darin, ein effektives und effizientes Kriterium zu finden, um zu bestimmen, wann diese über 7.000 Punkte nahe genug an einem regulären Raster liegen (da es möglich ist, dass sie in keiner Standardprojektion ein perfektes Raster bilden). Für einen schnellen Durchlauf der Datenbank sollte es ausreichen, eine kleine Anzahl von Entfernungen zu vergleichen, z. B. eine Ost-West-Entfernung im nördlichen Teil des Rasters mit einer Ost-West-Entfernung im südlichen Teil. Das sollte die große Mehrheit der Kandidaten schnell beseitigen.
Whuber
3
Ich habe alle von Mathematica 8 unterstützten (Standard-) Projektionen durchgearbeitet . Es wurde eine Projektion gefunden, bei der die Punkte tatsächlich auf ein Raster zu fallen scheinen: Alaska State Plane (1983) Zone 10! Dies ist eine Lambert Conformal Conic-Projektion. Ich glaube, es ist EPSG 26940 . Wenn Sie dies so ändern, dass es ungefähr bei -106 Längengrad zentriert ist, bilden die Punkte ein ziemlich gutes Gitter.
whuber
1
Abe, meinst du die Webseite zu lesen? Es war r = Import[ "https://ebi-forecast.igb.illinois.edu/bety/miscanthusyield.csv", "Data"];. Anschließend können Sie über data = Rest[r]; ListPlot[data[[;; , {3, 2}]]](oder ListPointPlot3D[data[[;; , {3, 2, 4}]]]) eine schnelle Darstellung der Punkte erhalten . Beginnen Sie bei Neuprojektionen mit der Hilfe von GeoGridPositionund machen Sie dann einige intelligente Vermutungen und Querverweise, um herauszufinden, was los ist :-). Übrigens, die Erklärung von @ Spacedman ist wirklich relevant: Die metrische Verzerrung von 25 bis 49 Grad entspricht cos (25) / cos (49) = 1,38; das ist erheblich.
whuber
29

Seit die Frage zuletzt beantwortet wurde, gibt es eine viel einfachere Lösung, indem die rasterFromXYZFunktion des Rasterpakets verwendet wird, in der alle erforderlichen Schritte (einschließlich der Angabe der CRS-Zeichenfolge) zusammengefasst sind. 

library(raster)
rasterFromXYZ(mydata)
Lucas Fortini
quelle
1
Entschuldigung an den unermüdlichen @Spacedman, der mir oft geholfen hat, aber ich denke, diese Antwort verdient es, die lustige grüne Zecke zu erben.
Geotheory
@ Theorie Ich würde diese Antwort auswählen, es ist eine großartige Funktion, aber es scheint sehr langsam auf dem Datensatz, den ich verwende (verknüpft mit im op)
Abe
1
... in der Tat erstickte es, weil es meine ~ 400KB-Datei nahm und eine Datei /tmp/mit ~ 19GB erstellte, als mir der Speicherplatz ausgegangen war.
Abe
Es gibt wahrscheinlich irgendwo einen n-Quadrat-Prozess. Möglicherweise können Sie die Punktedaten nach einem breiten Raster gruppieren, jede Gruppe einzeln rastern und dann merge()die Ergebnisse zusammenfassen.
Geotheorie
Bei allem Respekt, aber diese Antwort ist viel besser als die von Spacedman.
Ghost