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

#目的
地点から半径Xmのバッファー作って、それでラスタを切り抜く。(バッファーが丸くなくてもできるよ)
![Rplot0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/917414/980327dc-3b37-12d5-feb6-268f5b9ab9cb.png)

#環境
 - R version 3.6以上

#手順
[国土数値情報の土地利用細分メッシュ（ラスタ版)](https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-L03-b_r.html)の"L03-b-14_5339"を使ってやってみました。
 
![L03-b-14_5339のコピー.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/917414/032a6ba8-8d5d-533a-247b-2950d82fb0a9.png)

 1.バッファーをつくる

```r: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`

```r: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")
```
ただしこの段階では切り抜いたラスタは四角い。四隅がじゃまです。
![Rplot1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/917414/2ecb8ef7-f78b-833f-b569-4c01b7800d6d.png)
この四隅のセルにNAを入れることで消します。

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

```r: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のデータができます。

```r: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を入れます。
![Rplot00.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/917414/2682222b-1bba-b813-cf62-01e9b0040f9f.png)

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

![Rplot7.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/917414/e3c1cea2-a31a-8b8a-7a92-ccd39c8f16e3.png)

出来上がり:ramen:

![Rplot5.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/917414/79c786bb-9da7-e75b-8787-95f706eef95f.png)

ちなみに、QGISで作業したものとの比較。ラスタの再投影ところから少し違うのね。
`左：R`、`右:QGIS`
![QGIS0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/917414/a9b0fe59-a1b8-54b6-0a37-b7b26fb2950c.png)

#全部まとめてみた

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

```r: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
```