Extrahieren Sie Raster aus Raster mit Polygon-Shapefile in R

13

Ich bin neu in R und benutze das Raster-Paket. Beim Extrahieren von Polygonen aus einer vorhandenen Rasterdatei ist ein Problem aufgetreten. Wenn ich benutze

extract(raster, poly_shape)

Funktion auf dem Raster erstellt es immer eine Liste mit den Daten. Ich möchte wirklich eine weitere Raster-Datei extrahieren, die ich wieder mit ArcGIS laden kann. Nachdem ich ein bisschen mehr gelesen habe, denke ich, dass die Zuschneidefunktion genau das ist, was ich wirklich brauche. Aber wenn ich versuche, diese Funktion zu nutzen

crop(raster, poly_shape)

Ich bekomme diesen Fehler:

Error in .local(x, y, ...) : extents do not overlap
In addition: Warning message:
In intersect(extent(x), extent(y)) : Objects do not overlap

Die Dateien raster und poly_shape sind für beide Funktionen gleich. Können Sie mir sagen, was hier falsch sein könnte? Ist es überhaupt richtig, dass die Zuschneidefunktion ein anderes Raster und keine Liste erstellt?

BEARBEITEN : Die Funktion Extent () funktioniert bei mir nicht. Ich bekomme immer noch den gleichen Fehler. Ich bin mir aber sicher, dass sich die beiden Datensätze überschneiden! Mit dem

extract(raster, poly_shape)

Ich bekomme die richtigen Daten daraus. Nur als Liste und nicht als Raster, wie ich es haben möchte. Ich habe gerade die Datasets in ArcGIS geladen, und sie passen sehr gut, sodass ich die Projektion nicht überprüft habe. Jetzt habe ich es versucht

projection(raster) # "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs"
projection(poly_shape) # "+proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs"

und Sie können sehen, dass die Projektionen nicht passen. Die Extraktionsfunktion scheint in der Lage zu sein, die Dateien automatisch auf die richtige Weise zu transformieren. Ich weiß das, weil ich Folgendes getan habe:

  • Ich habe den genauen Teil des Polygons, das ich in R extrahiert habe, auch in ArcGIS ausgeschnitten
  • Ich habe die Summe aller Werte des extrahierten R-Polygons berechnet (Liste)
  • Ich habe die Summe aller Rasterzellen berechnet, die ich in ArcGIS ausgeschnitten habe

Die 2 haben genau das gleiche Ergebnis, daher sollte die Schlussfolgerung lauten, dass die Extraktionsfunktion korrekt funktioniert hat. Jetzt habe ich zwei Möglichkeiten:

  1. Ich brauche eine Möglichkeit, um einen Raster wieder aus der extrahierten Liste zu entfernen oder
  2. Die beiden Datensätze (Raster + Poly_Shape) müssen dieselbe Projektion verwenden und die Zuschneidefunktion sollte funktionieren

Was würden Sie hier vorschlagen?

Lars
quelle
Was ist, wenn es sich um ein 4-Band-RGBI-Raster handelt? Bands sind bis jetzt verloren ...
Doris
Willkommen bei der GIS SE! Nehmen Sie als neuer Benutzer unbedingt an der kurzen Tour teil . Überlegen Sie dann, ob Sie Ihre Antwort bearbeiten möchten, um zusätzliche Informationen und Verweise bereitzustellen. Siehe Wie schreibe ich eine gute Antwort? Für mehr Information.
Andy
Wenn Sie eine neue Frage haben, fragen Sie es bitte durch Klicken Frage stellen Taste. Fügen Sie einen Link zu dieser Frage hinzu, wenn dies zur Bereitstellung des Kontexts beiträgt. - Aus der Bewertung
Erik

Antworten:

19

Die Extraktionsfunktion verhält sich genau so, wie sie sollte. Sie können erzwingen, dass die Zuschneidefunktion die Ausdehnung des Polygons verwendet, und dann das Objekt maskieren, um das genaue Raster zurückzugeben, das den Polygonbereich darstellt. Wenn Sie den Fehler weiterhin erhalten, bedeutet dies, dass sich Ihre Daten tatsächlich nicht überschneiden.

Bitte denken Sie daran, dass R keine "on the fly" -Projektion durchführt. Überprüfen Sie daher Ihre Projektionen. Mit der Funktion "extend ()" können Sie überprüfen, ob sich Ihre Extents überschneiden.

Hier ist ein Beispiel für das Beschneiden mithilfe der Polygonausdehnung und das anschließende Maskieren des resultierenden Rasters mithilfe des "gerasterten" Polygons.

# Add required packages
require(raster)
require(rgeos)
require(sp)

# Create some data using meuse 
data(meuse)
  coordinates(meuse) <- ~x+y
    proj4string(meuse) <- CRS("+init=epsg:28992")
data(meuse.grid)
  coordinates(meuse.grid) = ~x+y
    proj4string(meuse.grid) <- CRS("+init=epsg:28992")
      gridded(meuse.grid) = TRUE    
        r <- raster(meuse.grid) 
          r[] <- runif(ncell(r))

# Create a polygon
f <- gBuffer(meuse[10,], byid=FALSE, id=NULL, width=250, 
                         joinStyle="ROUND", quadsegs=10)   

# Plot full raster and polygon                       
plot(r)
  plot(f,add=T)

# Crop using extent, rasterize polygon and finally, create poly-raster
#          **** This is the code that you are after ****  
cr <- crop(r, extent(f), snap="out")                    
  fr <- rasterize(f, cr)   
    lr <- mask(x=cr, mask=fr)

# Plot results
plot(lr)
  plot(f,add=T)
Jeffrey Evans
quelle
4
extract () führt eine sofortige Neuprojektion durch, crop () jedoch nicht. Das könnte für einige Verwirrung
sorgen
@Jefferey crop () und mask () schneiden das Raster nur entsprechend der rechteckigen Ausdehnung des Polygons ab, nicht jedoch innerhalb der Begrenzung des Polygons. Haben Sie eine Idee, mit welchen Befehlen das Raster innerhalb der Grenze des angegebenen Polygons abgeschnitten werden kann?
csheth
1
@Chintan Sheth, damit die Maske innerhalb des Polygons untergeordnet werden kann, muss ein Raster vorhanden sein, das die Werte innerhalb des Polygons darstellt. Aus diesem Grund rastern Sie das Teilmengenpolygon und maskieren es anschließend. Der Zuschneideschritt besteht darin, die Ausdehnung des Rasters auf das gleiche Maß wie das Polygon der Teilmenge zu reduzieren, das in der Vergangenheit, wenn es übersprungen wird, einen Ausdehnungsfehlanpassungsfehler hervorruft.
Jeffrey Evans
spTransformaus dem spPaket (das manchmal automatisch mit anderen räumlichen R-Paketen geladen wird) ermöglicht die Neuprojektion, so dass beide Dateien in der gleichen Projektion sind, z. good_poly=spTransform(spolygon, CRSobj=crs(raster_file))
user3386170
@ user3386170, nicht wahr? Ich bin mir nicht sicher, worauf du hinaus willst. Diese Frage trat zu einem Zeitpunkt auf, als das Rasterpaket innerhalb einiger seiner Funktionen "on the fly projection" hinzufügte. Diese Funktionen hatten zuvor einen Fehler ausgelöst, als die Projektionen nicht übereinstimmten, aber dieser Beitrag stammte aus dem Jahr 2014. Sie sollten auch beachten, dass rgdal immer geladen wird, wenn Sie sp :: spTransform () verwenden, da dies zusätzliche wichtige Funktionen hinter den Kulissen hinzufügt.
Jeffrey Evans
8

Was ich eigentlich gesucht habe, war die mask()Funktion.

mask(raster, poly_shape)

funktioniert fehlerfrei und gibt das zurück, wonach ich gesucht habe.

Lars
quelle
2
Projizieren Sie Ihre Daten so, dass sie sich im selben Projektionsbereich befinden. Sogar in ArcGIS, wo die On-the-Fly-Projektion automatisch erfolgt, ist es eine sehr schlechte Praxis, Analysen in verschiedenen Projektionen durchzuführen. Mit Daten in verschiedenen Projektionen teilen sich Ihre Daten keine gemeinsamen Bereiche. Dies ist der Fehler, den Sie erhalten.
Jeffrey Evans
Verwenden Sie projectExtent (), um genau die Ausdehnung zum Beschneiden des Rasters abzurufen.
mdsumner
Um dem Q & A-Format der Site zu entsprechen, sollte dies als Bearbeitung / Aktualisierung in den Hauptteil der Frage eingefügt werden (und dann ein Kommentar zu der Antwort hinzugefügt werden, auf die unter "Antworten" geantwortet wird, damit sie wissen, dass mehr zu betrachten ist).
Matt Wilkie
@mattwilkie Entschuldigung, dass das Format nicht passt, aber mein Text war zu lang, um ihn hier als Kommentar zu posten. @JeffreyEvans Ich habe versucht , die folgenden: projection(raster) = projection(poly_shape)und umgekehrt , projection(poly_shape) = projection(raster)aber beide Wege produzieren den gleichen Fehler: Error in .local(x, y, ...) : extents do not overlap In addition: Warning message: In intersect(extent(x), extent(y)) : Objects do not overlap. Gibt es eine Möglichkeit, wie ich mithilfe der Funktion extract () sehen kann, welche Projektion im laufenden Betrieb verwendet wird (weil diese offensichtlich funktioniert)?
Lars
1
Was ich eigentlich gesucht habe war die mask () Funktion. mask(raster, poly_shape)funktioniert fehlerfrei und gibt das zurück, wonach ich gesucht habe.
Lars
-1

Die Ausdehnung funktioniert einwandfrei ... Ich denke, die Xmin, Xmax, Ymin und Ymax Ihrer Ausdehnung unterscheiden sich von den X- und Y-Werten Ihres Rasters, dh sie sind entgegengesetzt gesetzt

ToNoY
quelle