R：创建加拿大选定省份和美国州的地图

Question

我正在尝试创建加拿大选定省份/地区和选定美国州的地图。到目前为止，最好的地图似乎是用GADM数据生成的地图：http://www.gadm.org/

但是，我还没有能够在同一地图上绘制美国和加拿大地图，或只绘制了选定的省/地区和州。例如，我对阿拉斯加州，育空地区，西北地区，不列颠哥伦比亚省，艾伯塔省和蒙大拿州等地感兴趣。

此外，美国地图似乎沿着国际日期划分。

可有人请帮我：

1. 情节上述省/地区和国家在一张地图
2. 避免美国沿国际日期变更线
3. 覆盖经纬度网格划分上
4. 选择一个特定的投影，也许是polyconic。

也许spplot不允许用户指定投影。我没有看到在spplot帮助页面上选择投影的选项。我知道如何在地图包中选择具有地图功能的投影，但这些地图看起来看起来不太好，我也无法绘制所需的省份/地区子集以及具有该功能的州。

我不知道如何开始添加经纬度网格。然而，“sp.pdf”文件的第3.2节似乎涉及该主题。

下面是我到目前为止的代码。我已经加载了我所遇到的每个与地图有关的软件包，并且除了省/地区或州边界之外，还注释了GADM数据。

遗憾的是，到目前为止，我只设法绘制加拿大或美国

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) 

感谢任何帮助的地图。到目前为止，我还没有取得上述步骤2-4的进展。也许我要求太多。也许我应该简单地切换到ArcGIS并尝试使用该软件。

我已阅读本StackOverflow的帖子：

Can R be used for GIS?

编辑

我现在借的电子副本 '中的空间数据分析与R' Bevand等。 （2008年）和下载（或位于）从本书的网站有关的R-代码和数据：

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

我也发现了一些好看的GIS相关的R代码在这里：

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

如果当我学习如何完成预期的目标，我会在这里发布解决方案。尽管如果我无法完成R中的目标，我最终可能会转移到ArcGIS。

Answer 1

要在同一设备上绘制多个SpatialPolygons对象，一种方法是指定您希望首先绘制的地理范围，然后使用plot(..., add=TRUE)。这只会增加地图上那些有趣的点。

使用投影（例如多锥投影）进行绘图要求首先使用rgdal包中的spTransform()函数确保所有图层都处于相同投影中。

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

添加其他功能地图，如突出感兴趣的司法管辖区，可以方便地使用相同的方法来完成：

​​

下面是结果：

当使用投影时添加网格线非常复杂，以至于需要另一个帖子，我想。看起来好像@Mark Miller在下面添加它！

Answer 2

下面我修改了PaulG出色的答案来显示经纬度网格。网格比我想要的粗糙，但可能是足够的。我在下面的代码中使用英国。我不知道如何将结果包含在这篇文章中。

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

结果如下：