5.5 座標参照系（CRS: coordinate reference system）

  • 地理座標系 (geographic coordinate system): 緯度経度で位置座標を定義
    • WGS84（World Geodesic System 1984）: 米国が構築した標準的な地理座標系
  • 投影座標系 (projected coordinate system): 平面に投影してkmなどで位置座標を定義する

5.6 sfパッケージによる空間データ処置の実例

パッケージ

# install.packages(c("sf", "NipponMap", "RColorBrewer"))

library(sf)
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1; sf_use_s2() is TRUE
library(NipponMap)
library(RColorBrewer)

データ

  • system.fileはパッケージ内部のファイルを探すための関数
shp <- system.file("shapes/jpn.shp", package = "NipponMap")

pref <- read_sf(shp)

pref
## Simple feature collection with 47 features and 5 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 127.6461 ymin: 26.0709 xmax: 148.8678 ymax: 45.5331
## CRS:           NA
## # A tibble: 47 × 6
##    SP_ID jiscode name      population region                            geometry
##    <chr> <chr>   <chr>          <dbl> <chr>                            <POLYGON>
##  1 1     01      Hokkaido     5506419 Hokkaido ((139.7707 42.3018, 139.8711 42.…
##  2 2     02      Aomori       1373339 Tohoku   ((140.8727 40.48187, 140.6595 40…
##  3 3     03      Iwate        1330147 Tohoku   ((140.7862 39.85982, 140.8199 39…
##  4 4     04      Miyagi       2348165 Tohoku   ((140.2802 38.01415, 140.2802 38…
##  5 5     05      Akita        1085997 Tohoku   ((140.7895 39.86026, 140.8253 39…
##  6 6     06      Yamagata     1168924 Tohoku   ((140.2802 38.01415, 140.2799 37…
##  7 7     07      Fukushima    2029064 Tohoku   ((140.2799 37.9721, 140.2799 37.…
##  8 8     08      Ibaraki      2969770 Kanto    ((139.7322 36.08471, 139.7058 36…
##  9 9     09      Tochigi      2007683 Kanto    ((139.4263 36.33584, 139.3703 36…
## 10 10    10      Gunma        2008068 Kanto    ((139.3836 36.35887, 139.6541 36…
## # ℹ 37 more rows
  • CRSの確認
st_crs(pref)
## Coordinate Reference System: NA
  • CRSをWGS84（EPSG:4326）に設定
st_crs(pref) <- 4326

st_crs(pref)
## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["World"],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]

投影変換

  • 距離の計算などを行うために、投影変換が必要。
  • 平面直角座標系第９系（EPSG: 6677）に変換
pref_tr <- st_transform(pref, crs = 6677)
  • 各県の重心をst_centroid(pref_tr)で計算。
  • メートル単位の位置座標をst_coordinates()で示す。
st_coordinates(st_centroid(pref_tr))
## Warning: st_centroid assumes attributes are constant over geometries
##                   X            Y
##  [1,]   196224.2824  810081.2252
##  [2,]    84655.7095  532108.5734
##  [3,]   131703.1173  399919.8404
##  [4,]    96017.9602  273316.1844
##  [5,]    48881.5628  416347.1935
##  [6,]    23755.3191  271399.0833
##  [7,]    34266.1553  153839.0704
##  [8,]    44215.9529   37128.0710
##  [9,]     -871.9157   76479.3353
## [10,]   -76310.4859   56936.7880
## [11,]   -43533.5526    -553.8963
## [12,]    34413.5114  -51484.9921
## [13,]   -34852.2698  -33855.7322
## [14,]   -43884.2433  -65626.1280
## [15,]   -79717.9893  162827.0732
## [16,]  -230935.0153   73201.4063
## [17,]  -274961.0384   86240.8624
## [18,]  -325340.5584  -12355.4692
## [19,]  -110559.2326  -41748.5331
## [20,]  -161840.0905   15641.9943
## [21,]  -251072.6811  -20199.7179
## [22,]  -136000.3959 -109305.3331
## [23,]  -238755.3609 -102430.2375
## [24,]  -317901.0856 -155946.0509
## [25,]  -337500.5177  -79848.9249
## [26,]  -399277.1317  -74798.3965
## [27,]  -396392.6852 -144278.6496
## [28,]  -456294.9490  -90138.1773
## [29,]  -365480.6626 -180323.9483
## [30,]  -399112.6329 -223922.1975
## [31,]  -542967.9991  -53780.8032
## [32,]  -670574.0069  -85899.8114
## [33,]  -550683.7390 -106026.7720
## [34,]  -645936.7624 -128587.4491
## [35,]  -764207.8268 -167439.4100
## [36,]  -515993.9098 -217340.9137
## [37,]  -539725.3675 -183320.0065
## [38,]  -647794.7435 -243432.0426
## [39,]  -600797.2096 -267587.4156
## [40,]  -852329.4713 -235718.1884
## [41,]  -906344.6457 -258580.0314
## [42,]  -928055.9891 -293559.2984
## [43,]  -844978.0139 -336080.8817
## [44,]  -783834.7892 -280424.3105
## [45,]  -806696.3552 -391493.7123
## [46,]  -876460.2743 -447720.1656
## [47,] -1189252.4588 -999282.1039

5.7 sfパッケージによる地図化の実例

  • 沖縄を除く都道府県データの作成
pref2 <- pref[pref$name != "Okinawa",]

都道府県人口の地図化（1）

plot(pref2[, "population"])

都道府県人口の地図化（2）

  • パレットの一覧
display.brewer.all()

  • 色分け数 nc = 7, nbreaks = nc と設定
  • 配色 RdYlGn (Red-Yellow-Green)
  • axes = True: 緯度経度を表示
nc <- 7 #色分け数 nc = 7
pal <- brewer.pal(nc, "RdYlGn") # 配色 RdYlGn (Red-Yellow-Green)
plot(pref2[,"population"], pal = pal, axes = TRUE, nbreaks = nc)

都道府県人口の地図化（3）

  • 色分けの閾値を設定
breaks <- c(0, 1000000, 2000000, 3000000, 5000000, 8000000, max(pref2$population) )
  • 色の数を閾値の数-1となる
 nc <- length(breaks) - 1
  • 配色を反転させる
pal <- rev(brewer.pal(nc, "RdYlGn"))
  • 地図
plot(pref2[,"population"], pal = pal, axes = TRUE, breaks = breaks)

  • 地図（凡例の位置を下部に移動、小さくする）
plot(pref2[,"population"], pal = pal, axes = TRUE, breaks = breaks,
     key.pos = 1, key.length = 0.8)

出典

  • 村上大輔 (2022)『Rではじめる地理空間データの統計解析入門』講談社サイエンティフィク。
  • サポートページ