r - R：選択したカナダの州と米国の州の地図を作成する

Question

選択したカナダの州/準州と選択した米国の州の地図を作成しようとしています。これまでのところ、最も優れたマップはGADMデータで生成されたマップのようです：http：//www.gadm.org/

ただし、米国とカナダを同じ地図にプロットしたり、選択した州/準州および州のみをプロットしたりすることはできませんでした。たとえば、アラスカ、ユーコン、NWT、ブリティッシュコロンビア、アルバータ、モンタナなどに興味があります。

また、米国の地図は国際日付変更線に沿って分割されているようです。

誰かが私を助けてくれますか：

1. 前述の州/準州および州を単一の地図にプロットします
2. 米国が国際日付変更線に沿って分割されることは避けてください
3. 緯度経度グリッドをオーバーレイする
4. 特定の投影法、多円錐図法を選択します。

たぶん、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の投稿を読みました：

RはGISに使用できますか？

編集

私は今、「Rによる応用空間データ分析」の電子コピーを借りました。（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)。これにより、関心のあるポイントのみがマップに追加されます。

多円錐図法などの投影法を使用してプロットするには、最初にパッケージ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)

結果は次のとおりです。

投影法を使用するときにグリッド線を追加するのは十分に複雑なので、別の投稿が必要になると思います。以下に追加された@MarkMillerのように見えます！

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

結果は次のようになります。