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

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

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

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

**【前回記事】**
[【Rstats事始め】JPGIS2.1をマップ表示する](https://qiita.com/m37335/items/df6f54b97c2ab85299d1)

# コード全体像
```r
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設定の確認については以下のサイトを参考にしました。

https://qiita.com/iwasaki_kenichi/items/e2f0c6777263a93c7230

## STEP1 Tiffファイルの準備
ラスターデータであるtifファイルは以下から取得した神奈川県付近のものを使用します。
[ASTER GDEM](https://gdemdl.aster.jspacesystems.or.jp/ "ASTER GDEM")
![ダウンロード｜ASTER_GDEM.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/867981/6cc925c6-1c28-14d3-80d9-1d91bdea2d6a.png)

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

**参考記事**

https://qiita.com/m37335/items/df6f54b97c2ab85299d1

## STEP2 ファイルの読み込み
### ラスターデータ
上記のサイトよりダウンロードしたTIFFファイルを用いる。
```r
# 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()
```
**出力されたデータ**
![kanagawa_DEM.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/867981/8baa9ec8-85d4-4dbc-4f34-27e1ea23fcb7.png)

#### CRSの確認
```r
# CRSの確認
st_crs(kanagawa)
```
**出力**
```r
> 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をマップ表示する](https://qiita.com/m37335/items/df6f54b97c2ab85299d1 "【Rstats事始め】JPGIS2.1をマップ表示する")


![kanagawa.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/867981/ab8f2634-32c8-9482-a18f-5091cb1f5761.png)

### readOGRで読み込んだ場合
```r
# Import Vector layer
kanagawa_geoData= rgdal::readOGR("[filePath]/kanagawa_N03-20220101_14_GML/N03-22_14_220101.shp")
```

#### CRSの確認
```r
crs(kanagawa_geoData)
```
**出力**
```r
> 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で読み込んだ場合
```r
kanagawa_sf <- sf::read_sf("[FilePath]/kanagawa_N03-20220101_14_GML/N03-22_14_220101.shp",
                        options= c("ENCODING=CP932"),
                        quiet = FALSE)
```

#### CRSの確認
```r
crs(kanagawa_sf)
```
**出力**
```r
> 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を用います。
```r
'新データ' <- st_transform('元データ', crs = proj4string('CRSを合わせるデータ'))
```

### （1） readOGRで読み込んだ場合
```r
> kanagawa_geoData <- st_transform(kanagawa_geoData, crs = proj4string(kanagawa))
 UseMethod("st_transform") でエラー: 
   'st_transform' をクラス "c('SpatialPolygonsDataFrame', 'SpatialPolygons', 'Spatial', 'SpatialVector')" のオブジェクトに適用できるようなメソッドがありません 
```
ジオメタリタイプによってはエラーが出力される。

### （2） sf::read_sfで読み込んだ場合
```r
> 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をもつデータを重ね、マップを作成することができました。
![kanagawa_map_DEM.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/867981/72520e5a-21d4-e043-4703-9760d3fc8e15.png)