R: creación de un mapa de provincias canadienses seleccionadas y estados de EE. UU.

Question

Estoy intentando crear un mapa de provincias/territorios canadienses seleccionados y estados de EE. UU. Seleccionados. Hasta ahora, los mejores mapas parecen ser los generados con datos de GADM: http://www.gadm.org/

Sin embargo, no he podido trazar los Estados Unidos y Canadá en el mismo mapa o trazar solo las provincias/territorios y estados seleccionados. Por ejemplo, estoy interesado en Alaska, Yukon, NWT, Columbia Británica, Alberta y Montana, entre otros.

Además, el mapa de EE. UU. Parece estar dividido en la fecha internacional.

por favor alguien puede ayudar a:

1. parcela de las mencionadas provincias/territorios y estados en un solo mapa
2. evitar tener los EE.UU. dividido a lo largo de la línea de fecha internacional
3. superposición de una cuadrícula de latitud y longitud
4. seleccione una proyección específica, tal vez la policónica.

Tal vez spplot no permite a los usuarios especificar las proyecciones. No vi una opción para seleccionar una proyección en la página de ayuda de spplot. Sé cómo seleccionar proyecciones con la función de mapa en el paquete de mapas, pero esos mapas no parecían tan bonitos y tampoco pude trazar el subconjunto deseado de provincias/territorios y estados con esa función.

No sé cómo comenzar a agregar una cuadrícula de latitud y longitud. Sin embargo, la Sección 3.2 del archivo 'sp.pdf' parece abordar el tema.

A continuación se muestra el código que he encontrado hasta ahora. He cargado todos los paquetes relacionados con mapas que he encontrado y comentado los datos de GADM a excepción de los límites provinciales/territoriales o estatales.

Por desgracia, hasta ahora sólo han logrado trazar mapas de Canadá o los EE.UU.

library(maps) 
library(mapproj) 
library(mapdata) 
library(rgeos) 
library(maptools) 
library(sp) 
library(raster) 
library(rgdal) 

# can0<-getData('GADM', country="CAN", level=0) # Canada 
    can1<-getData('GADM', country="CAN", level=1) # provinces 
# can2<-getData('GADM', country="CAN", level=2) # counties 

plot(can1)  
spplot(can1, "NAME_1") # colors the provinces and provides 
         # a color-coded legend for them 
can1$NAME_1   # returns names of provinces/territories 
# us0 <- getData('GADM', country="USA", level=0) 
    us1 <- getData('GADM', country="USA", level=1) 
# us2 <- getData('GADM', country="USA", level=2) 
plot(us1)    # state boundaries split at 
         # the dateline 
us1$NAME_1    # returns names of the states + DC 
spplot(us1, "ID_1") 
spplot(us1, "NAME_1") # color codes states and 
         # provides their names 
# 
# Here attempting unsuccessfully to combine U.S. and Canada on one map. 
# Attempts at selecting given states or provinces have been unsuccessful. 
# 
plot(us1,can1) 
us.can1 <- rbind(us1,can1) 

Gracias por cualquier ayuda. Hasta ahora no he progresado con los pasos 2 a 4 anteriores. Tal vez estoy pidiendo demasiado. Tal vez debería simplemente cambiar a ArcGIS y probar ese software.

He leído este post StackOverflow:

Can R be used for GIS?

EDITAR

ahora me han prestado una copia electrónica de 'Applied espacial Análisis de datos con R' Bevand et al. (2008) y descargado (o localizado) R código y datos relacionados del sitio web del libro:

http://www.asdar-book.org/

También encontré algo de código R de buen aspecto afín a los SIG aquí:

https://sites.google.com/site/rodriguezsanchezf/news/usingrasagis

Si y cuando aprendo cómo lograr los objetivos deseados, publicaré soluciones aquí. Aunque eventualmente puedo moverme a ArcGIS si no puedo lograr los objetivos en R.

Answer 1

Para trazar múltiples objetos SpatialPolygons en el mismo dispositivo, un enfoque es especificar la extensión geográfica que desea trazar primero, y luego usar plot(..., add=TRUE). Esto agregará al mapa solo los puntos que sean de interés.

Para trazar utilizando una proyección, (por ejemplo, una proyección policónica) se requiere primero usar la función spTransform() en el paquete rgdal para asegurarse de que todas las capas estén en la misma proyección.

## Specify a geographic extent for the map 
## by defining the top-left and bottom-right geographic coordinates 
mapExtent <- rbind(c(-156, 80), c(-68, 19)) 

## Specify the required projection using a proj4 string 
## Use http://www.spatialreference.org/ to find the required string 
## Polyconic for North America 
newProj <- CRS("+proj=poly +lat_0=0 +lon_0=-100 +x_0=0 
      +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs") 

## Project the map extent (first need to specify that it is longlat) 
mapExtentPr <- spTransform(SpatialPoints(mapExtent, 
        proj4string=CRS("+proj=longlat")), 
        newProj) 

## Project other layers 
can1Pr <- spTransform(can1, newProj) 
us1Pr <- spTransform(us1, newProj) 

## Plot each projected layer, beginning with the projected extent 
plot(mapExtentPr, pch=NA) 
plot(can1Pr, border="white", col="lightgrey", add=TRUE) 
plot(us1Pr, border="white", col="lightgrey", add=TRUE) 

Adición de otras características para el mapa, como destacando las jurisdicciones de interés, se puede hacer fácilmente utilizando el mismo enfoque:

## Highlight provinces and states of interest 
theseJurisdictions <- c("British Columbia", 
         "Yukon", 
         "Northwest Territories", 
         "Alberta", 
         "Montana", 
         "Alaska") 

plot(can1Pr[can1Pr$NAME_1 %in% theseJurisdictions, ], border="white", 
    col="pink", add=TRUE) 

plot(us1Pr[us1Pr$NAME_1 %in% theseJurisdictions, ], border="white", 
    col="pink", add=TRUE) 

aquí está el resultado:

Agregar líneas de cuadrícula cuando se utiliza una proyección es lo suficientemente complejo como para requerir otra publicación, creo. ¡Parece como si @Mark Miller lo hubiera agregado a continuación!

Answer 2

A continuación he modificado la respuesta destacada de PaulG para mostrar una cuadrícula de latitud y longitud. La grilla es más gruesa de lo que quisiera, pero podría ser adecuada. Yo uso el Reino Unido con el siguiente código. No sé cómo incluir el resultado en esta publicación.

library(rgdal) 
library(raster) 

# define extent of map area 
mapExtent <- rbind(c(0, 62), c(5, 45)) 

# BNG is British National Grid  
newProj <- CRS("+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.999601271625 
      +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs") 

mapExtentPr <- spTransform(SpatialPoints(mapExtent, 
        proj4string=CRS("+proj=longlat")), 
        newProj) 

# provide a valid 3 letter ISO country code 
# obtain a list with: getData("ISO3") 
uk0 <- getData('GADM', country="GBR", level=0) # UK 
uk1 <- getData('GADM', country="GBR", level=1) # UK countries 
uk2 <- getData('GADM', country="GBR", level=2) # UK counties 

# United Kingdom projection 
uk1Pr  <- spTransform(uk1, newProj) 

# latitude-longitude grid projection 
grd.LL  <- gridlines(uk1, ndiscr=100) 
lat.longPR <- spTransform(grd.LL, newProj) 

# latitude-longitude text projection 
grdtxt_LL <- gridat(uk1) 
grdtxtPR <- spTransform(grdtxt_LL, newProj) 

# plot the map, lat-long grid and grid labels 
plot(mapExtentPr, pch=NA) 
plot(uk1Pr, border="white", col="lightgrey", add=TRUE) 
plot(lat.longPR, col="black", add=TRUE) 
text(coordinates(grdtxtPR), 
    labels=parse(text=as.character(grdtxtPR$labels))) 

El resultado será similar a: