データ分析メモと北欧生活

旧Untitled Note. データ分析、計量経済・統計とR、水産管理、英語勉強、海外生活などについて備忘録や自分の勉強のOutputの場所として

MENU

【R】sfパッケージでCRS情報のないファイルをあるファイルと一緒に使いたい

RでGISといえば最近はsfパッケージらしいですね。
spとsfの違いとか、そもそもCRSのとは?みたいなレベルながら空間情報を扱いたくて、少しずつ使い始めています。

今回やりたいのは、すでにある地図になにか緯度経度の情報のみのデータを重ねたいときに、ちょっと躓いて調べたので、その備忘録です。

まず、きちんとした CRS情報があるデータがあるとします。今回は私が使ってるノルウェーの県別のデータなんですが、一般公開データです。

library(sf)
library(ggplot2)

# Norwegian Map
fylkermap = st_read("Fylker18.geojson") 

st_crs(fylkermap)


st_crsはsfオブジェクトのCRSを調べる関数です。

> st_crs(fylkermap)
Coordinate Reference System:
  EPSG: 25833 
  proj4string: "+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

と、このようにCRS情報が含まれています。

ちなみにこのオブジェクトで地図を描くとこういう感じ。

ggplot(fylkermap) +
  geom_sf(aes(fill = navn))

f:id:keita43a:20181023015231p:plain

県別に色をつけています。

さて、今回はこれにノルウェー近海の漁業区画のデータを重ねたいのです。
漁業区画のデータはシェープファイルで用意してあります。
これも一般公開データです。
https://yggdrasil.fiskeridir.no/

このマップデータですが、CRS情報がなく、各行に漁業区画のセルのgeometryが緯度経度で記録されているだけです。

map = st_read("494.shp")

st_crs(map)
> st_crs(map)
Coordinate Reference System: NA

CRS を指定しないまま描画するとこんな感じ。

ggplot(map) + geom_sf()

f:id:keita43a:20181023224737p:plain

一応縦横に緯度経度通り並べてある感じ。


ここでつまづいていたのですが、どうやら緯度経度を基準とするCRSの情報をまず渡してあげるとよいみたい。
以下のようにCRSを指定します。

map2 = map

st_crs(map2) = "+proj=longlat +datum=WGS84"

ggplot(map2) + geom_sf()

CRS情報を渡してから描画すると、もう少し縦横の緯度経度の感じを考慮したアウトプットになる。しかし緯度経度がそのままXとY軸になってるのは変わらず。


このCRS情報をもった緯度経度の図を、ノルウェーの地図に重ねるには同じCRSに揃えないと行けない。
ノルウェーの地図の方にCRSを揃えたいので、CRSを保存してから漁業区画の方のCRSを変換する。

crs1 = st_crs(fylkermap)

map3 = st_transform(map2, crs = crs1)

変換した区画マップはCRS情報がノルウェーマップとおなじになっている。

> st_crs(map3)
Coordinate Reference System:
  EPSG: 25833 
  proj4string: "+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
> st_crs(fylkermap_simple)
Coordinate Reference System:
  EPSG: 25833 
  proj4string: "+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

あとはggplotのgeom_sfを使って重ねて描画する。

ggplot() +
  geom_sf(data = fylkermap,fill = "lightblue") +
  geom_sf(data = map3,fill = NA)

f:id:keita43a:20181023225539p:plain


これで2つを重ねて描画する事ができた。
地図の穴が開いてるところは、アイスランドとイギリス・アイルランドだと思うけど、それらのマップは手元にないのでまた今度・・・。


緯度経度情報のみでCRS情報を合わせたいときにまずCRS情報を付与してから変換するとうまくいくよ、という話でした。