Rでラスターをまるく切り抜きたい

#目的
地点から半径Xmのバッファー作って、それでラスタを切り抜く。(バッファーが丸くなくてもできるよ)

#環境

• R version 3.6以上

#手順
国土数値情報の土地利用細分メッシュ（ラスタ版)の"L03-b-14_5339"を使ってやってみました。

1.バッファーをつくる

buffer.r

library(sf)
library(dplyr)
library(raster)
#地点の緯度経度(ここでは座標参照系JSD2000)
  p0=c(139.25,35.75) 
#平面直角座標系に再投影(epsg:4126→3100)
  p1=st_point(p0) %>% st_sfc() %>% st_set_crs(4126) %>% st_transform(3100) 
#地点p1を中心に半径1000mのバッファーを作成して、座標系をp1と同じに設定。
  buf0=st_buffer(p1,dist=1000,endCapStyle = "ROUND") %>% st_transform(crs=crs(p1,asText=TRUE))#半径lmのバッファーポリゴン

2.ラスタファイル読み込みして、バッファーで切り抜き

• projectRasterのmethod
• 土地利用は、ラスタの値がカテゴリ的なものを示しているので、セルの値を最近傍ngbにしました
• 標高や気温のように連続値ならば、bilinear

raster-crop.r

#ラスタファイル読み込み
  r1=raster("L03-b-14_5339.tif") #JSD2000
#平面直角座標系に再投影、p1と同じにしている
  r2=projectRaster(r1,crs=crs(p1,asText=TRUE),method="ngb")

#バッファーで切り抜き(crop)
  rc1=crop(r2,st_as_sf(buf0),snap="out")

ただしこの段階では切り抜いたラスタは四角い。四隅がじゃまです。

この四隅のセルにNAを入れることで消します。

　3.ラスタのセル番号をもつpointを作る

raster.r

  rc2<-rc1 #ラスタのコピーを作って
  names(rc2)<-"no"
  values(rc2)[1:ncell(rc2)]<-c(1:ncell(rc2)) #中身をセル番号にします

#ラスタの各セルの中心の座標を取り出して、セル番号-ポイントのsfクラスのファイルに変換
  pc1=rasterToPoints(rc2) %>% as.data.frame() %>% st_as_sf(coords=c("x","y"),crs=st_crs(buf0)) 

こういうpointのデータができます。

res1.r

Simple feature collection with 418 features and 1 field
geometry type:  POINT
dimension:      XY
bbox:           xmin: 340771.1 ymin: 3956646 xmax: 342805.1 ymax: 3958586
CRS:            +proj=utm +zone=54 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
First 10 features:
   no                 geometry
1   1 POINT (340771.1 3958586)
2   2 POINT (340884.1 3958586)
3   3 POINT (340997.1 3958586)
4   4 POINT (341110.1 3958586)
##以下略##

4.pointをバッファーで切り抜く
バッファーの外にあるpointと対応するラスタのセルNAを入れます。

point-crop.r

  pc2=st_intersection(pc1,buf0) #警告メッセージ:出たけど気にしない
#pc2はバッファーの中にあるラスターのセル番号をもつ。
#バッファー外側にあるラスターのセルにNAを入れる
 values(rc1)[-pc2$no]<-NA #pc2$noはバッファー内のセル番号
 rc1

バッファーの中のpointと重なってないセルにNAを入れたわけです

出来上がり

ちなみに、QGISで作業したものとの比較。ラスタの再投影ところから少し違うのね。
左：R、右:QGIS

#全部まとめてみた

もうちょっと行数減らしてみた

buffer-crop.r

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

#ポイント
  p0=c(139.25,35.75) 

#平面直角座標系に再投影(epsg:4126→3100) → 半径lmのバッファーポリゴン
  buf0=st_point(p0) %>% st_sfc() %>% st_set_crs(4126) %>% st_transform(3100) %>% st_buffer(dist=1000,endCapStyle = "ROUND")
#ラスタファイル読み込み → バッファーで切り抜き(crop)
  r1=raster("L03-b-14_5339.tif") %>% projectRaster(crs=crs(buf0,asText=TRUE),method="ngb") %>% crop(st_as_sf(buf0),snap="out")

  rc1<-r1 #ラスタのコピーを作って
  names(rc1)<-"no"
  values(rc1)[1:ncell(r1)]<-c(1:ncell(r1)) #中身をセル番号にします

#ラスタの各セルの中心の座標を取り出して、セル番号-ポイントのsfクラスのファイルに変換 → バッファーで切り抜き
  pc1=rasterToPoints(rc1) %>% as.data.frame() %>% st_as_sf(coords=c("x","y"),crs=st_crs(buf0)) %>% st_intersection(buf0)
 #警告メッセージ:出たけど気にしない

#バッファー外側にあるラスタのセルにNAを入れる
 values(r1)[-pc1$no]<-NA 
 r1