Beschneiden eines Rasters in R

32

Ich baue eine Karte für den Nordosten der USA. Der Kartenhintergrund muss entweder eine Höhenkarte oder eine mittlere Jahrestemperaturkarte sein. Ich habe zwei Raster von Worldclim.org, die mir diese Variablen geben, aber ich muss sie auf das Ausmaß der Zustände zuschneiden, an denen ich interessiert bin. Vorschläge, wie man das macht. Das habe ich bisher:

#load libraries
library (sp)
library (rgdal)
library (raster)
library (maps)
library (mapproj)


#load data
state<- data (stateMapEnv)
elevation<-raster("alt.bil")
meantemp<-raster ("bio_1.asc")

#build the raw map
nestates<- c("maine", "vermont", "massachusetts", "new hampshire" ,"connecticut",
  "rhode island","new york","pennsylvania", "new jersey",
  "maryland", "delaware", "virginia", "west virginia")

map(database="state", regions = nestates, interior=T,  lwd=2)
map.axes()

#add site localities
sites<-read.csv("sites.csv", header=T)
lat<-sites$Latitude
lon<-sites$Longitude

map(database="state", regions = nestates, interior=T, lwd=2)
points (x=lon, y=lat, pch=17, cex=1.5, col="black")
map.axes()
library(maps)                                                                  #Add axes to  main map
map.scale(x=-73,y=38, relwidth=0.15, metric=T,  ratio=F)

#create an inset map

 # Next, we create a new graphics space in the lower-right hand corner.  The numbers are proportional distances within the graphics window (xmin,xmax,ymin,ymax) on a scale of 0 to 1.
  # "plt" is the key parameter to adjust
    par(plt = c(0.1, 0.53, 0.57, 0.90), new = TRUE)

  # I think this is the key command from http://www.stat.auckland.ac.nz/~paul/RGraphics/examples-map.R
    plot.window(xlim=c(-127, -66),ylim=c(23,53))

  # fill the box with white
    polygon(c(0,360,360,0),c(0,0,90,90),col="white")

  # draw the map
    map(database="state", interior=T, add=TRUE, fill=FALSE)
    map(database="state", regions=nestates, interior=TRUE, add=TRUE, fill=TRUE, col="grey")

Die Elevations- und Meanemp-Objekte müssen an der Flächenausdehnung des Nestates-Objekts abgeschnitten werden. Jede Eingabe würde helfen

I Del Toro
quelle
2
Gibt es eine Möglichkeit, dies für andere reproduzierbar zu machen, indem Sie möglicherweise Raster aus zufälligen Daten mit dem gleichen Umfang und der gleichen Auflösung erstellen?
Spacedman

Antworten:

38

Ich würde das mapsPaket verwenden und ein Status-Shapefile finden. Laden Sie das dann mit in R rgdal, und führen Sie dann einige Polygon-Overlay-Arbeiten aus.

library(raster)
# use state bounds from gadm website:
# us = shapefile("USA_adm1.shp")
us <- getData("GADM", country="USA", level=1)
# extract states (need to uppercase everything)
nestates <- c("Maine", "Vermont", "Massachusetts", "New Hampshire" ,"Connecticut",
         "Rhode Island","New York","Pennsylvania", "New Jersey",
         "Maryland", "Delaware", "Virginia", "West Virginia")

ne = us[match(toupper(nestates),toupper(us$NAME_1)),]


# create a random raster over the space:        
r = raster(xmn=-85,xmx=-65,ymn=36,ymx=48,nrow=100,ncol=100)
r[]=runif(100*100)

# plot it with the boundaries we want to clip against:
plot(r)
plot(ne,add=TRUE)

# now use the mask function
rr <- mask(r, ne)

# plot, and overlay:
plot(rr);plot(ne,add=TRUE)

Wie ist das? Das Gadm-Shapefile ist ziemlich detailliert. Vielleicht möchten Sie stattdessen ein verallgemeinertes finden.

Raumfahrer
quelle
Prost Robert, nette Bearbeitung. Ich glaube, ich hatte die Maske vergessen.
Spacedman
32

Hier ist ein Ansatz extract()aus dem rasterPaket. Ich habe es mit Höhen- und Durchschnittstemperaturdaten von der WorldClim- Website getestet (ich beschränke dieses Beispiel auf die Höhe, die Temperatur funktioniert ähnlich), und ein entsprechendes Shapefile der USA mit Staatsgrenzen finden Sie hier . Laden Sie einfach die .zip-Daten herunter und dekomprimieren Sie sie in Ihr Arbeitsverzeichnis.

Sie müssen rgdalund rasterBibliotheken laden , um fortzufahren.

library(rgdal)
library(raster)

Importieren wir nun das US-Shapefile mit readOGR(). Nachdem ich das CRS des Shapefiles festgelegt habe, erstelle ich eine Untermenge mit den gewünschten Zuständen. Achten Sie auf Groß- und Kleinschreibung!

state <- readOGR(dsn = path.data, layer = "usa_state_shapefile")
projection(state) <- CRS("+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs")

# Subset US shapefile by desired states
nestates <- c("Maine", "Vermont", "Massachusetts", "New Hampshire" ,"Connecticut",
             "Rhode Island","New York","Pennsylvania", "New Jersey",
             "Maryland", "Delaware", "Virginia", "West Virginia")

state.sub <- state[as.character(state@data$STATE_NAME) %in% nestates, ]

Importieren Sie anschließend die Rasterdaten mit raster()und beschneiden Sie sie mit der Ausdehnung der zuvor generierten Statusuntermenge.

elevation <- raster("/path/to/data/alt.bil")

# Crop elevation data by extent of state subset
elevation.sub <- crop(elevation, extent(state.sub))

Als letzten Schritt müssen Sie die Pixel Ihres Höhenrasters identifizieren, die innerhalb der Grenzen der angegebenen Zustandspolygone liegen. Verwenden Sie dazu die Funktion 'mask'.

elevation.sub <- mask(elevation.sub, state.sub)

Hier ist eine sehr einfache Darstellung der Ergebnisse:

plot(elevation.sub)
plot(state.sub, add = TRUE)

DEM der nordöstlichen US-Bundesstaaten

Prost,
Florian

fdetsch
quelle
Woher hast du das staatliche Shapefile?
I Del Toro
@IDelToro, ich habe es von Geocommons bekommen .
Fdetsch
Warum dauert dies so lange (>> 15 Minuten, möglicherweise Stunden), wenn Sie mit einer ~ 11-MB-Rasterebene und einer Datei mit einer einzelnen Polygonform arbeiten? Gibt es eine effizientere Methode?
ecologist1234
@ ecologist1234, kannst du ein Beispiel geben?
Fdetsch