Recortar un ráster en R

32

Estoy construyendo un mapa para el noreste de EE. UU. El fondo del mapa debe ser un mapa de altitud o un mapa de temperatura anual promedio. Tengo dos rasters de Worldclim.org que me dan estas variables, pero necesito recortarlas en la medida de los estados en los que estoy interesado. Cualquier sugerencia sobre cómo hacerlo. Esto es lo que tengo hasta ahora:

#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")

Los objetos de elevación y significados son los que deben recortarse a la extensión del área del objeto anidado. Cualquier entrada ayudaría

    
pregunta I Del Toro 20.05.2013 - 05:07

2 respuestas

37

Yo soltaría usar el paquete maps y encontraría un shapefile de estado. Luego, cárguelo en R usando rgdal , y luego haga un trabajo de superposición de polígonos.

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)

¿Cómo es eso? El archivo shape de gadm es bastante detallado, es posible que desee encontrar uno más generalizado.

    
respondido por el Spacedman 20.05.2013 - 13:31
31

Aquí hay un enfoque que usa extract() del paquete raster . Lo probé con altitud y mean temperature datos de WorldClim sitio web (limito este ejemplo a la altitud, la temperatura funciona de manera similar), y se puede encontrar un shapefile apropiado de los límites estatales del estado de EE. UU. aquí . Simplemente descargue los datos .zip y descomprímalos en su directorio de trabajo.

Es necesario cargar las bibliotecas rgdal y raster para poder continuar.

library(rgdal)
library(raster)

Ahora importemos el shapefile de EE. UU. usando readOGR() . Después de configurar el CRS del shapefile, creo un subconjunto que contiene los estados deseados. ¡Presta atención al uso de mayúsculas y minúsculas iniciales!

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([email protected]$STATE_NAME) %in% nestates, ]

A continuación, importe los datos ráster utilizando raster() y recórtelos con la extensión del subconjunto de estados previamente generado.

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

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

Como paso final, necesitas identificar esos píxeles de tu ráster de elevación que se encuentran dentro de los bordes de los polígonos de estado dados. Utilice la función 'máscara' para eso.

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

Aquí hay una trama muy simple de los resultados:

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

Saludos,
Florian

    
respondido por el fdetsch 22.05.2013 - 17:02

Lea otras preguntas en las etiquetas