R: создание карты отдельных провинций Канады и штатов США

Question

Я пытаюсь создать карту отдельных провинций Канады/территорий и выбранных штатов США. До сих пор наиболее красивыми были карты, сгенерированные с данными GADM: http://www.gadm.org/

Однако я не смог построить США и Канаду на одной карте или отобразить только выбранные провинции/территории и штаты. Например, меня интересуют Аляска, Юкон, СЗТ, Британская Колумбия, Альберта и Монтана.

Кроме того, карта США, как представляется, разделена по международной линии передачи данных.

Может кто-то пожалуйста, помогите мне:

1. участка вышеупомянутых провинций/территорий и государств на одной карте
2. избежать необходимости США расколоть по международной линии смены дат
3. наложения широты-долготы сетки
4. выберите конкретную проекцию, возможно, поликоническую.

Возможно, spplot не позволяет пользователям указывать проекции. Я не видел выбора для выбора проекции на странице справки spplot. Я знаю, как выбирать проекции с помощью функции карты в пакете карт, но эти карты выглядели не так красиво, и я не мог отобразить желаемое подмножество провинций/территорий и состояний с этой функцией.

Я не знаю, как начать добавление сетки широты-долготы. Однако в разделе 3.2 файла «sp.pdf», по-видимому, рассматривается тема.

Ниже приведен код, который я привел до сих пор. Я загрузил каждый связанный с картой пакет, на который я наткнулся, и прокомментировал данные 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?

EDIT

Я теперь заимствованы электронную копию 'Applied пространственного анализа данных с R' Bevand и др. (2008) и загрузить (или находится), связанный R-код и данные с сайта книги:

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

Я также нашел несколько хороших перспективных ГИС, связанные с R-код здесь:

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

Если и когда я узнаю, как достичь желаемых целей, я буду размещать решения здесь. Хотя я, возможно, в конечном итоге перейду в ArcGIS, если я не смогу выполнить цели в R.

Answer 1

Чтобы построить несколько объектов SpatialPolygons на одном устройстве, один из подходов состоит в том, чтобы указать географическую протяженность, которую вы хотите построить сначала, а затем с помощью plot(..., add=TRUE). Это добавит на карту только те точки, которые представляют интерес.

Для того чтобы убедиться, что все слои находятся в одной и той же проекции, необходимо использовать функцию 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) 

Добавление других функций на карте, например, выделение юрисдикций, представляющих интерес, можно легко сделать, используя тот же подход:

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

Вот результат:

Добавление линий сетки при использовании проекции достаточно сложно, поскольку для этого требуется другая запись. Похоже, что @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))) 

Результат выглядит следующим образом: