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

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

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

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

Кто-нибудь может помочь мне:

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

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

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

Ниже приведен код, который я привел до сих пор. Я загрузил каждый связанный с картой пакет, на который я наткнулся, и прокомментировал данные GADM, за исключением провинциальных/территориальных или государственных границ.

К сожалению, до сих пор мне удалось отобразить карты Канады или U.S.

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:

Может ли R использоваться для ГИС?

ИЗМЕНИТЬ

Теперь я заимствовал электронную копию "Прикладного анализа пространственных данных с R" Бевандом и др. (2008) и загрузил (или разместил) соответствующий R-код и данные с веб-сайта книги:

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

Я также нашел некоторый красивый G-код, связанный с GIS:

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

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

Ответ 1

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

Построение с использованием проекции (например, поликонической проекции) требует сначала использования функции spTransform() в пакете rgdal, чтобы убедиться, что все слои находятся в одной и той же проекции.

## 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 добавил это ниже!

Ответ 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)))

Результат выглядит так: