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

結果如下：