R: tworzenie mapy wybranych kanadyjskich prowincji i stanów USA

Question

Próbuję utworzyć mapę wybranych kanadyjskich prowincji/terytoriów i wybranych stanów USA. Jak dotąd najpiękniejsze są mapy wygenerowane danymi GADM: http://www.gadm.org/

Jednak nie byłem w stanie wykreślić USA i Kanady na tej samej mapie lub działce tylko wybranych prowincji/terytoriów i stanów. Na przykład interesują mnie Alaska, Yukon, NWT, Kolumbia Brytyjska, Alberta i Montana.

Ponadto wydaje się, że mapa USA jest podzielona wzdłuż międzynarodowej linii danych.

Czy ktoś mógłby mi pomóc:

1. działki wyżej wymienione prowincje/Terytoria i stwierdza na jednej mapie
2. uniknąć USA podzielona wzdłuż Linia Zmiany Daty
3. nakładki siatce szerokość-długość geograficzna
4. wybierz konkretną projekcję, być może polikoniczną.

Może spplot nie pozwala użytkownikom na określenie rzutów. Nie widziałem opcji wyboru projekcji na stronie pomocy spplot. Wiem, jak wybrać projekcje za pomocą funkcji map w pakiecie map, ale te mapy nie wyglądały tak ładnie i nie mogłem wykreślić pożądanego podzbioru prowincji/terytoriów i stanów z tą funkcją.

Nie wiem, jak zacząć dodawać siatkę szerokości i szerokości geograficznej. Jednak sekcja 3.2 pliku "sp.pdf" wydaje się poruszać ten temat.

Poniżej znajduje się kod, który wymyśliłem do tej pory. Załadowałem każdy pakiet związany z mapą, na które natknąłem się i skomentowałem dane GADM z wyjątkiem granic prowincjonalnych/terytorialnych lub państwowych.

Niestety, do tej pory udało mi się jedynie do wykreślenia mapy

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) 

dzięki za pomoc Kanadzie lub Stanach Zjednoczonych. Do tej pory nie zrobiłem żadnego postępu w Kroku 2 - 4 powyżej. Być może proszę o zbyt wiele. Być może powinienem po prostu przejść na ArcGIS i wypróbować to oprogramowanie.

Znam ten StackOverflow wpis:

Can R be used for GIS?

EDIT

Mam teraz pożyczył elektroniczną kopię 'Applied przestrzennego Analiza danych z R' Bevand et al. (2008) i pobrać (lub siedzibę) związany kod R i dane ze strony internetowej książki:

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

Znalazłem też jakieś fajne wyglądające związanych z GIS-kod R tutaj:

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

Jeśli i kiedy nauczę się, jak osiągnąć pożądane cele, zamieszczę tutaj rozwiązania. Chociaż w końcu mogę przejść do ArcGIS, jeśli nie mogę osiągnąć celów w R.

Answer 1

Aby wykreślić wiele obiektów SpatialPolygons na tym samym urządzeniu, jedną z metod jest określenie zakresu geograficznego, który ma zostać wykreślony jako pierwszy, a następnie użycie plot(..., add=TRUE). Spowoduje to dodanie do mapy tylko tych punktów, które są interesujące.

kreślenie przy projekcji, (na przykład polyconic projekcji) wymaga najpierw przy użyciu funkcji spTransform() w opakowaniu rgdal, czy wszystkie te warstwy są w tym samym rzucie.

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

Dodawanie innych funkcji do mapy, takie jak wyróżnianie jurysdykcje interesów, można łatwo zrobić za pomocą tego samego podejścia:

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

Oto wynik:

Dodanie linii siatki, gdy używana jest projekcja, jest wystarczająco skomplikowane, że wymaga innego postu, jak sądzę. Wygląda na to, że @Mark Miller dodał go poniżej!

Answer 2

Poniżej zmodyfikowałem zaległą odpowiedź PaulG, aby wyświetlić siatkę szerokości i szerokości geograficznej. Siatka jest grubsza niż chciałbym, ale może być wystarczająca. Używam Zjednoczonego Królestwa z poniższym kodem. Nie wiem, jak zawrzeć wynik w tym poście.

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

Rezultat wygląda następująco: