【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))
県別に色をつけています。
さて、今回はこれにノルウェー近海の漁業区画のデータを重ねたいのです。
漁業区画のデータはシェープファイルで用意してあります。
これも一般公開データです。
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()
一応縦横に緯度経度通り並べてある感じ。
ここでつまづいていたのですが、どうやら緯度経度を基準とする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)
これで2つを重ねて描画する事ができた。
地図の穴が開いてるところは、アイスランドとイギリス・アイルランドだと思うけど、それらのマップは手元にないのでまた今度・・・。
緯度経度情報のみでCRS情報を合わせたいときにまずCRS情報を付与してから変換するとうまくいくよ、という話でした。