【備忘録】RにおけるCRSの設定について

はじめに

Rを用いてGIS分析を行うにあたってCRSの設定が必要です。CRSは位置を記述するための標準的な方法です。地理データを記述するために、多くの異なるCRSが使用されています。選択されたCRSは、データが収集された時期、データの地理的範囲、データの目的などに応じて選択されます。

RでGISの分析を行う際には異なるCRSを持つデータを組み合わせる場合、共通のCRSに変換することが重要です。

今回はCRSの確認、設定についての備忘録になります。
記事内で使用している神奈川県のマップデータについては以下をご参照ください。

【前回記事】
【Rstats事始め】JPGIS2.1をマップ表示する

コード全体像

library(raster)
library(sf)
library(ggplot2)
# Package for colorpallet
library(viridis) 
library(cptcity)

#--------------------------
# DRAFT: CHECK CRS
#--------------------------
# Open raster layer
kanagawa = raster("[FilePath]/ASTGTMV003_N35E139_dem.tif")
kanagawa_mat = raster_to_matrix(kanagawa)
st_crs(kanagawa)

# Import Vector layer
## readOGR
kanagawa_geoData= rgdal::readOGR("[FilePath]/kanagawa_N03-20220101_14_GML/N03-22_14_220101.shp")
## CRSの確認
crs(kanagawa_geoData)
## CRSの変換
kanagawa_geoData <- st_transform(kanagawa_geoData, crs = proj4string(kanagawa))

## read_sf
kanagawa_sf <- sf::read_sf("[FilePath]/kanagawa_N03-20220101_14_GML/N03-22_14_220101.shp",
                        options= c("ENCODING=CP932"),
                        quiet = FALSE)
## CRSの確認
crs(kanagawa_sf)

#--------------------------
# DRAFT: TRANSFORM CRS
#--------------------------
## CRSの変換
kanagawa_sf <- st_transform(kanagawa_sf, crs = proj4string(kanagawa))

#--------------------------
# DRAFT: CREATE MAP
#--------------------------
## crop
kanagawa_all <- crop(kanagawa, extent(kanagawa_sf))

## Check that it worked
plot(kanagawa_all)

# ggplotで描画するため、tibble形式に変換
kanagawa_all_dem <-
  kanagawa_all %>%   # ラスターDEM
  as.data.frame( xy = TRUE ) %>% # rasterをデータフレームに変換
  tibble::as_tibble() %>% # tipple形式に変換
  dplyr::rename("Elevation" = ASTGTMV003_N35E139_dem)  # 標高値をElevationにリネーム

# Map出力
ggplot() +
  theme_void() +
  geom_raster(data = kanagawa_all_dem ,   # 標高ラスタ
              aes(x, y, fill = Elevation),
              hjust = 0, vjust = 0) +
  scale_fill_gradientn(colours = cptcity::cpt(pal = "tp_tpushum" , n = 50)) + # カラーパレットにcptを利用
  geom_contour(data = kanagawa_all_dem, aes(x, y, z = Elevation), col = "gray70", size = 0.2) +  # 等高線描画
  geom_sf(data = kanagawa_sf , col = "gray40" , fill = NA, lwd = 0.3) +
  theme_void() +
  scale_colour_viridis(option = "cividis" , discrete =TRUE ) +
  guides(fill="none") 

CRSの確認と設定

以下のサイトを参考に各ファイルのCRS設定を確認する。
CRS設定の確認については以下のサイトを参考にしました。

STEP1 Tiffファイルの準備

ラスターデータであるtifファイルは以下から取得した神奈川県付近のものを使用します。
ASTER GDEM

一方、ベクターデータは以前作成した国土数値情報からJPGIS2.1をダウンロードしたshpファイルを使用しました。

参考記事

STEP2 ファイルの読み込み

ラスターデータ

上記のサイトよりダウンロードしたTIFFファイルを用いる。

# Open raster layer
kanagawa = raster("/Users/masa/Downloads/Download_N35E139/ASTGTMV003_N35E139/ASTGTMV003_N35E139_dem.tif")
kanagawa_mat = raster_to_matrix(kanagawa)

# 画像データの縮小
kanagawa_small = resize_matrix(kanagawa_mat,0.25)

kanagawa_small %>% 
  height_shade() %>% 
  plot_map()

出力されたデータ

CRSの確認

# CRSの確認
st_crs(kanagawa)

出力

> st_crs(kanagawa)
# 座標参照系
Coordinate Reference System:
  User input: WGS 84 (with axis order normalized for visualization) 
  wkt:
# 地理座標系 a geographic coordinate system.
GEOGCRS["WGS 84 (with axis order normalized for visualization)",
    # DATUM:測地系
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    # PRIMEM: 本初子午線 グリニッジ標準時
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    # Ellipsoidal coordinate systems: 楕円体
    CS[ellipsoidal,2],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]],
    REMARK["Axis order reversed compared to EPSG:4326"]]

ベクターデータ

以前作成したものを使用する。
作成の流れについては以下の記事を参照ください。
【Rstats事始め】JPGIS2.1をマップ表示する

readOGRで読み込んだ場合

# Import Vector layer
kanagawa_geoData= rgdal::readOGR("[filePath]/kanagawa_N03-20220101_14_GML/N03-22_14_220101.shp")

CRSの確認

crs(kanagawa_geoData)

出力

> crs(kanagawa_geoData)
# 座標参照系
Coordinate Reference System:
Deprecated Proj.4 representation:
 +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs 
WKT2 2019 representation:
# 地理座標系 a geographic coordinate system.
GEOGCRS["JGD2011",
    # DATUM:測地系
    DATUM["Japanese Geodetic Datum 2011",
        ELLIPSOID["GRS 1980",6378137,298.257222101,
            LENGTHUNIT["metre",1]]],
    # PRIMEM: 本初子午線 グリニッジ標準時
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    # Ellipsoidal coordinate systems: 楕円体
    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["Horizontal component of 3D system."],
        AREA["Japan - onshore and offshore."],
        BBOX[17.09,122.38,46.05,157.65]],
    ID["EPSG",6668]] 

sf::read_sfで読み込んだ場合

kanagawa_sf <- sf::read_sf("[FilePath]/kanagawa_N03-20220101_14_GML/N03-22_14_220101.shp",
                        options= c("ENCODING=CP932"),
                        quiet = FALSE)

CRSの確認

crs(kanagawa_sf)

出力

> crs(kanagawa_sf)
# 座標参照系
Coordinate Reference System:
Deprecated Proj.4 representation:
 +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs 
WKT2 2019 representation:
# 地理座標系 a geographic coordinate system.
GEOGCRS["JGD2011",
    # DATUM:測地系
    DATUM["Japanese Geodetic Datum 2011",
        ELLIPSOID["GRS 1980",6378137,298.257222101,
            LENGTHUNIT["metre",1]]],
    # PRIMEM: 本初子午線 グリニッジ標準時
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    # Ellipsoidal coordinate systems: 楕円体
    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["Horizontal component of 3D system."],
        AREA["Japan - onshore and offshore."],
        BBOX[17.09,122.38,46.05,157.65]],
    ID["EPSG",6668]] 

STEP3 CRSの設定

CRSの変換はst_transfromを用います。

'新データ' <- st_transform('元データ', crs = proj4string('CRSを合わせるデータ'))

（1） readOGRで読み込んだ場合

> kanagawa_geoData <- st_transform(kanagawa_geoData, crs = proj4string(kanagawa))
 UseMethod("st_transform") でエラー: 
   'st_transform' をクラス "c('SpatialPolygonsDataFrame', 'SpatialPolygons', 'Spatial', 'SpatialVector')" のオブジェクトに適用できるようなメソッドがありません 

ジオメタリタイプによってはエラーが出力される。

（2） sf::read_sfで読み込んだ場合

> kanagawa_sf <- st_transform(kanagawa_sf, crs = proj4string(kanagawa))
> crs(kanagawa_sf)
Coordinate Reference System:
Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs 
WKT2 2019 representation:
GEOGCRS["unknown",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ID["EPSG",6326]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8901]],
    CS[ellipsoidal,2],
        AXIS["longitude",east,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]],
        AXIS["latitude",north,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]] 

問題なく、変換できました。

STEP4 マップの出力

無事、異なるCRSをもつデータを重ね、マップを作成することができました。