3次メッシュデータをラスタにしたい

目的

国土数値情報の平年値(気候)メッシュなどのメッシュデータがシェイプファイルになっているものをラスタにしたい。

手順

国土数値情報の平年値メッシュの"5034のデータ(G02-12_5034-jgd.shp)"を使ってやってみました。
1月降水量(G02_002)をラスタ化します。

1. シェイプファイルの読み込み

shp01.r

library(sf)
library(raster)
#shpファイル読み込み
  s1 <- read_sf("G02-12_5034-jgd.shp")
#shpファイルのcrs設定
  s1 <- st_set_crs(s1,4612)
#bounding boxの値を取得
  s1_b <- st_bbox(s1)

shpファイルの中身。85列にメッシュ番号と気候値が入っています。
st_bboxはbounding boxの値、つまり地理的な範囲、東西南北の端の緯度経度の値を取得します。

shp02.r

> s1
Simple feature collection with 3005 features and 84 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 134 ymin: 33.33333 xmax: 134.825 ymax: 34
Geodetic CRS:  JGD2000
# A tibble: 3,005 × 85
   G02_001  G02_002 G02_003 G02_004 G02_005 G02_006 G02_007

> s1_b
     xmin      ymin      xmax      ymax 
134.00000  33.33333 134.82500  34.00000 

2. メッシュの数を計算してラスタ作成

raster::rasterは、列数・行数・地理的な範囲が分かればラスタが作れる。
メッシュデータなので、メッシュ1個の大きさが分かれば南北・東西方向に何個並んでいるか(列数と行数)わかる。

shp03.r

s2 <- s1[1,] #メッシュを1個だけ取り出す
s2_b <- st_bbox(s2) #メッシュを1個の地理的な範囲、すなわちmaxとminの差がメッシュのサイズ

ここから、
- 東西方向のメッシュ数
(s1_b[3]-s1_b[1])/(s2_b[3]-s2_b[1])
- 南北方向のメッシュ数(南北方向のメッシュの長さは割り切れない数なので、四捨五入する)
round((s1_b[4]-s1_b[2])/(s2_b[4]-s2_b[2]),0)
で計算できます。

南北・東西のメッシュ数(nrows, ncols)・地理的な範囲(s1_b)・CRSを入力して、セル(グリッド)のサイズがshpファイルと同じになるラスタを作成。
番号の入ったラスタと、何も値が入っていないラスタを作ります。
何も値が入っていないラスタは、海や広い湖沼が範囲内にあるデータをラスタ化するときに便利(海や湖沼上は気候値がない)。

rast01.r

  r1 <- raster(nrows=round((s1_b[4]-s1_b[2])/(s2_b[4]-s2_b[2]),0),
               ncols=(s1_b[3]-s1_b[1])/(s2_b[3]-s2_b[1]),
               xmn=s1_b[1], xmx=s1_b[3], ymn=s1_b[2], ymx=s1_b[4], 
               crs=CRS("+init=epsg:4612"))
#nrows: 南北方向のメッシュ数
#ncols: 東西方向のメッシュ数
#xmn,xmx,ymn,ymx: shpファイルの地理的範囲
#crs: crsの設定、shpファイルと同じにする

##警告メッセージが出るかもしれないが、うまくできてれば無視##

#セルの値には、1から全セル数までを順番に入れておく：セル番号
  values(r1) <- 1:ncell(r1)
#作ったラスタをコピーして、セルに値が入っていないラスタも作成
  r0 <- r1
  values(r0) <-NA

ラスタr1の中身

rast02.r

> r1
class      : RasterLayer 
dimensions : 80, 66, 5280  (nrow, ncol, ncell)
resolution : 0.0125, 0.008333333  (x, y)
extent     : 134, 134.825, 33.33333, 34  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs 
source     : memory
names      : layer 
values     : 1, 5280  (min, max)

3. メッシュのデータを同じ位置のラスタのセルに入れていく

セル番号をr1から取得して、気候値データをr0に入れていく。

rast03.r

#shpファイルの方のメッシュの中心座標を取得
  s1_c=st_centroid(s1)  
#メッシュの中心(s1_c)と重なるセルのセル番号(no)をr1から取って、気候データと付けて、データフレームにする
  s1_d=data.frame(as.data.frame(s1_c),no=extract(r1,s1_c))

#メッシュ１つ１つの気候データを、r0の対応する位置のセルの値として入れていく
  for(i in 1:nrow(s1_d)){
    values(r0)[s1_d$no[i]] <- s1_d$G02_002[i]/10
}
##G02_002には1月降水量のデータが整数で入っているので、全て10で割る

出来上がり

rast04.r

> r0
class      : RasterLayer 
dimensions : 80, 66, 5280  (nrow, ncol, ncell)
resolution : 0.0125, 0.008333333  (x, y)
extent     : 134, 134.825, 33.33333, 34  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs 
source     : memory
names      : layer 
values     : 22.1, 106.6  (min, max)

全部まとめてみた

mesh.r

library(sf)
library(raster)
library(dplyr)

s1 <- read_sf("G02-12_5034-jgd.shp") %>% st_set_crs(4612)
s1_b <- st_bbox(s1)

s2 <- s1[1,]
s2_b <- st_bbox(s2)

r1 <- raster(nrows=round((s1_b[4]-s1_b[2])/(s2_b[4]-s2_b[2]),0),
             ncols=(s1_b[3]-s1_b[1])/(s2_b[3]-s2_b[1]),
             xmn=s1_b[1], xmx=s1_b[3], ymn=s1_b[2], ymx=s1_b[4],
             crs=CRS("+init=epsg:4612"))
  values(r1) <- 1:ncell(r1)
  r0 <- r1
  values(r0) <- NA

s1_d <- st_centroid(s1) %>% data.frame(no=extract(r1,.))

  for(i in 1:nrow(s1_d)){
    values(r0)[s1_d$no[i]] <- s1_d$G02_002[i]/10
   }

出力してQGISで重ねて、うまくできているか確認してみると良いかも。