1 準備編

1.1 インストール

以下の環境を用意します。

  • Windows OS >=7
  • R >= 3.5.1
  • RStudio >= 1.0

使用するデータのダウンロードとテキストを用意します。

パッケージをインストールします。 RStudioを起動して、コンソールに以下を貼り付けて実行します。

install.packages(c("sf","geos","dplyr","tidyr","raster","ggplot2","ggrepel","mapview","jpndistrict","randomForest","kernlab"),dependencies=T,type="win.binary")

ライブラリが読み込めるか確認します。
RStudioのコンソールに以下を貼り付けて実行します。エラーメッセージが出なければOKです。
「～マスクされています：」というメッセージは問題ありません。

library(sf) #ベクタ処理
library(raster) #ラスタ処理
library(dplyr) #属性処理
library(purrr) #属性処理
library(tidyr) #属性処理
library(ggplot2) #地図表示
library(ggrepel) #ラベル表示処理
library(mapview) #地図表示（インタラクティブ）
library(jpndistrict) #日本の行政界データ
library(rpart) #決定木
library(randomForest) #ランダムフォレスト
library(kernlab) #SVM

1.2 プロジェクトの作成

RStudioで新規プロジェクトを作成します。
ダウンロードしたデータを解凍し「R_GIS_2018」フォルダを”日本語パス名ではない”場所に置きます。
File > New Projext… > Existing Directory > 「R_GIS_data2018」を指定し「Create Project」を押します。

1.3 基礎知識

1.3.1 パイプ

パイプ左側の関数の出力結果をパイプ右側の関数の第一引数にする

#head(1:5,3)をパイプで書くと
1:5 %>% head(3)
[1] 1 2 3
#パイプの出力を次のパイプの入力にできる
c(1,10,100) %>% tan %>% sin %>% cos
[1] 0.5403777 0.8231382 0.8504039
#.を使うと第1引数以外にできる
3 %>% head(1:5,.)
[1] 1 2 3

1.3.2 map関数

ベクトルかリストを入力し、各要素に関数を適用し、入力と同じ長さのリストを返す

#xを平均値とした正規乱数を5個作成する関数を定義
myfunc<-function(x){
  rnorm(5,x)
}
#各値を引数として、myfunc関数を適用する。戻り値は、各5個の乱数を持った3要素のリスト。
c(1,10,100) %>% map(myfunc)
[[1]]
[1] -0.5811198  0.7399071  0.6681385 -0.4134890  1.5514651

[[2]]
[1] 10.287854 10.007775  8.025936  9.513119 10.354516

[[3]]
[1] 101.20028 100.64113 101.27250  99.90563  99.63272
#以下のように簡略化して書ける。~はfunction、.xは引数を表す。
c(1,10,100) %>% map(~rnorm(5,.x))
#リストの各要素ごとに平均値を計算する場合。戻り値は、身長の平均値、体重の平均値の2要素のリスト
mylist<-list(height=130:140,weight=40:50)
mylist %>% map(mean)
$height
[1] 135

$weight
[1] 45

1.3.3 reduce関数

ベクトルかリストを入力し、二変数関数を順々に適用して1つの値を返す。

c(1,10,100) %>% reduce(sum)
[1] 111
x <- list(c(0, 1), c(2, 3), c(4, 5))
x %>% reduce(c)
[1] 0 1 2 3 4 5

1.3.4 データの変形

横長なデータを作成

data <- tibble(
  Point = c("地点A","地点B"),
  Red = c(0.1,0.3),
  Green = c(0.3,0.2),
  Blue = c(0.5,0.6)
)
print(data)
# A tibble: 2 x 4
  Point   Red Green  Blue
  <chr> <dbl> <dbl> <dbl>
1 地点A   0.1   0.3   0.5
2 地点B   0.3   0.2   0.6

縦長のデータに変換

data2<-data %>% gather(key = Band, value=Refrectance,Red,Green,Blue)

横長のデータに変換

data2 %>% spread(key=Band,value=Refrectance)
# A tibble: 2 x 4
  Point  Blue Green   Red
  <chr> <dbl> <dbl> <dbl>
1 地点A   0.5   0.3   0.1
2 地点B   0.6   0.2   0.3

2 ベクタ基本編

2.1 基本操作

2.1.1 シェープファイルを読み込みたい

landcover<-read_sf("data/landcover.shp",options=c("ENCODING=CP932"))

2.1.2 データの概要を確認したい

head(landcover)
Simple feature collection with 6 features and 2 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 24904.97 ymin: 2710.929 xmax: 24986.65 ymax: 2733.338
epsg (SRID):    NA
proj4string:    +proj=tmerc +lat_0=36 +lon_0=139.8333333333333 +k=0.9999 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs
# A tibble: 6 x 3
  class  subclass                                                 geometry
  <chr>  <chr>                                               <POLYGON [m]>
1 構造物 コンクリート ((24985.42 2713.855, 24986.65 2713.886, 24986.65 27~
2 構造物 コンクリート ((24933.71 2711.357, 24932.91 2712.066, 24933.68 27~
3 植生   芝           ((24921.19 2711.823, 24920.89 2711.823, 24920.48 27~
4 植生   高木         ((24917.46 2711.022, 24917.61 2711.126, 24917.85 27~
5 構造物 コンクリート ((24925.27 2710.929, 24925.06 2710.985, 24924.73 27~
6 植生   高木         ((24917.46 2711.022, 24910.18 2710.995, 24910.08 27~

2.1.3 地図を表示したい

plot(landcover)

2.1.4 データクラスを確認したい

class(landcover)
[1] "sf"         "tbl_df"     "tbl"        "data.frame"

2.1.5 データの構造を確認したい

str(landcover) 
Classes 'sf', 'tbl_df', 'tbl' and 'data.frame': 157 obs. of  3 variables:
 $ class   : chr  "構造物" "構造物" "植生" "植生" ...
 $ subclass: chr  "コンクリート" "コンクリート" "芝" "高木" ...
 $ geometry:sfc_POLYGON of length 157; first list element: List of 1
  ..$ : num [1:5, 1:2] 24985 24987 24987 24986 24985 ...
  ..- attr(*, "class")= chr  "XY" "POLYGON" "sfg"
 - attr(*, "sf_column")= chr "geometry"
 - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA
  ..- attr(*, "names")= chr  "class" "subclass"

2.1.6 属性名を確認したい

names(landcover)
[1] "class"    "subclass" "geometry"

2.1.7 属性値のサマリーを確認したい

summary(landcover)
    class             subclass                  geometry  
 Length:157         Length:157         POLYGON      :157  
 Class :character   Class :character   epsg:NA      :  0  
 Mode  :character   Mode  :character   +proj=tmer...:  0

2.1.8 データの範囲を確認したい

st_bbox(landcover)
     xmin      ymin      xmax      ymax 
24904.970  2710.929 24992.022  2824.005

2.1.9 Rstudioで属性表を確認したい

View(landcover)

2.1.10 インタラクティブな地図を表示したい

mapview(landcover)

2.1.11 日本の行政界を表示したい

library(jpndistrict)
ibaraki<-jpn_pref(admin_name="茨城県")
plot(ibaraki)

2.1.12 シェープファイルを書き出したい

st_write(landcover, "mydata/landcover2.shp",layer_options = "ENCODING=UTF-8",delete_layer = TRUE)
Deleting layer `landcover2' using driver `ESRI Shapefile'
Writing layer `landcover2' to data source `C:\Users\mizutani\Desktop\R_GIS_tutorial\mydata\landcover2.shp' using driver `ESRI Shapefile'
options:        ENCODING=UTF-8 
features:       157
fields:         2
geometry type:  Polygon

2.2 座標系の操作

2.2.1 データの空間参照系を確認したい

st_crs(landcover) 
Coordinate Reference System:
  No EPSG code
  proj4string: "+proj=tmerc +lat_0=36 +lon_0=139.8333333333333 +k=0.9999 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs"

2.2.2 空間参照系をEPSGコードで指定したい

landcover_epsg<-st_set_crs(landcover,2451)

2.2.3 別の座標系に変換したい

landcover_latlon<-st_transform(landcover,4612)

2.2.4 データのジオメトリを確認したい

st_geometry(landcover)
Geometry set for 157 features 
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 24904.97 ymin: 2710.929 xmax: 24992.02 ymax: 2824.005
epsg (SRID):    NA
proj4string:    +proj=tmerc +lat_0=36 +lon_0=139.8333333333333 +k=0.9999 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs
First 5 geometries:

2.3 データの作成・編集

2.3.1 ポイントデータを作成したい

p1<-st_point(c(1,2)) #sfgクラス
p2<-st_point(c(2,1))
p1
p2
mygeom<-st_sfc(p1,p2) #sfcクラス
mygeom
Geometry set for 2 features 
geometry type:  POINT
dimension:      XY
bbox:           xmin: 1 ymin: 1 xmax: 2 ymax: 2
epsg (SRID):    NA
proj4string:    NA
myname <- c("point1","point2") #属性値を入力
point_f<-st_sf(name=myname,geom=mygeom) #sfクラス
point_f
Simple feature collection with 2 features and 1 field
geometry type:  POINT
dimension:      XY
bbox:           xmin: 1 ymin: 1 xmax: 2 ymax: 2
epsg (SRID):    NA
proj4string:    NA
    name        geom
1 point1 POINT (1 2)
2 point2 POINT (2 1)

2.3.2 ラインデータを作成したい

pts1 = rbind(c(-1,-1), c(1,1), c(2,0), c(3,2))
pts2 = rbind(c(1,-1), c(2,2), c(3,1), c(2,4))
ls1 = st_linestring(pts1)
ls2 = st_linestring(pts2)
mygeom<-st_sfc(ls1,ls2) #sfcクラス
myname <- c("line1","line2") #属性値を入力
line_f<-st_sf(myname,mygeom) #sfクラス

2.3.3 ポリゴンデータを作成したい

pts = rbind(c(-1,-1), c(1,-1), c(1,1), c(-1,1), c(-1,-1))
b0 = st_polygon(list(pts)) #sfgクラス
b1 = b0 + c(2,0) #b0のポリゴンのX座標を2ずらしたポリゴンを作成
b2 = b0 * 0.5 + c(0,1.5) #b0のポリゴンを縮小してY座標を1.5ずらしたポリゴンを作成
mygeom<-st_sfc(b0,b1,b2) #sfcクラス
myname<-data.frame(number=c(1,2,3),species=c("dog","cat","cat")) #属性データを作成
poly_f<-st_sf(myname,mygeom) #sfクラス

2.3.4 ランダムポイントを作成したい

pnt<-st_sample(landcover,1000) #landcoverの範囲に1000点
plot(pnt,pch=".",cex=2)

2.3.5 メッシュを作成したい

mesh <- st_make_grid(n=c(10,5), offset = c(24920,2740), cellsize = c(6,6),crs = 2451)
mesh <- st_sf(mesh)
plot(mesh,axes=T)

2.3.6 ポリゴンをラインに変換したい

landcover_l<-landcover %>% st_cast("LINESTRING")

2.3.7 sfクラスからspクラスに変換したい

landcover_sp<-as(landcover,"Spatial")

2.3.8 spクラスからsfクラスに変換したい

landcover_sf<-st_as_sf(landcover_sp)

2.4 地図表示

2.4.1 指定したデータを表示したい

plot(landcover["subclass"])

2.4.2 座標軸を表示したい

plot(landcover["subclass"],axes=T)

2.4.3 凡例の表示幅を指定したい

plot(landcover["subclass"],key.width = lcm(3))

2.4.4 境界だけを表示したい

plot(st_geometry(landcover))

2.4.5 ggplot2で表示したい

ggplot() + 
  geom_sf(data = landcover,aes(fill = subclass))

2.4.6 座標軸の座標系を指定したい

ggplot() + 
  geom_sf(data = landcover,aes(fill = subclass)) +
  coord_sf(datum = 2451)

2.4.7 座標軸と背景を消したい

ggplot() + 
  geom_sf(data = landcover,aes(fill = subclass)) +
  coord_sf(datum = NA) +
  theme_void()

2.4.8 凡例の色を指定したい

mycol <- c("gray", "black", "white", "darkgreen","brown","lightgreen","lightblue", "blue","orange","green","yellow")
ggplot() + 
  geom_sf(data = landcover,aes(fill = subclass)) +
  scale_fill_manual(values = mycol) + 
  coord_sf(datum = NA) +
  theme_void()

2.4.9 データを重ねたい

ggplot() + 
  geom_sf(data = landcover,aes(fill = subclass)) +
  geom_sf(data = pnt,size=0.5) +
  geom_sf(data = mesh,fill = "transparent",color="white",size=1)

2.4.10 ラベルを表示したい

参照:属性にポイントの座標を追加したい

cent_xy<-landcover %>% st_centroid() %>% st_coordinates()
landcover_xy<-cbind(landcover,cent_xy)
ggplot(data=landcover_xy) + 
  geom_sf(aes(fill = subclass))+
  geom_text_repel(aes(x=X,y=Y,label=subclass),size=2)

2.5 データの選択・検索

2.5.1 特定の行のデータを抽出したい

landcover[1,]
Simple feature collection with 1 feature and 2 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 24985.42 ymin: 2712.775 xmax: 24986.65 ymax: 2713.886
epsg (SRID):    NA
proj4string:    +proj=tmerc +lat_0=36 +lon_0=139.8333333333333 +k=0.9999 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs
# A tibble: 1 x 3
  class  subclass                                                 geometry
  <chr>  <chr>                                               <POLYGON [m]>
1 構造物 コンクリート ((24985.42 2713.855, 24986.65 2713.886, 24986.65 27~

2.5.2 特定の列のデータを抽出したい

landcover[,"class"]
Simple feature collection with 157 features and 1 field
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 24904.97 ymin: 2710.929 xmax: 24992.02 ymax: 2824.005
epsg (SRID):    NA
proj4string:    +proj=tmerc +lat_0=36 +lon_0=139.8333333333333 +k=0.9999 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs
# A tibble: 157 x 2
   class                                                          geometry
   <chr>                                                     <POLYGON [m]>
 1 構造物 ((24985.42 2713.855, 24986.65 2713.886, 24986.65 2712.806, 2498~
 2 構造物 ((24933.71 2711.357, 24932.91 2712.066, 24933.68 2712.868, 2493~
 3 植生   ((24921.19 2711.823, 24920.89 2711.823, 24920.48 2711.781, 2492~
 4 植生   ((24917.46 2711.022, 24917.61 2711.126, 24917.85 2711.195, 2491~
 5 構造物 ((24925.27 2710.929, 24925.06 2710.985, 24924.73 2711.118, 2492~
 6 植生   ((24917.46 2711.022, 24910.18 2710.995, 24910.08 2711.091, 2490~
 7 裸地   ((24991.99 2718.085, 24991.99 2717.693, 24990.16 2717.663, 2498~
 8 植生   ((24991.99 2717.693, 24991.99 2710.986, 24931.58 2710.985, 2491~
 9 植生   ((24905.07 2714.537, 24905.01 2720.457, 24905.18 2720.342, 2490~
10 裸地   ((24970.5 2719.856, 24970.72 2719.979, 24971.08 2719.979, 24971~
# ... with 147 more rows

2.5.3 特定の列のデータをベクトルで抽出したい

landcover$class
 [1] "構造物" "構造物" "植生"   "植生"   "構造物" "植生"   "裸地"  
 [8] "植生"   "植生"   "裸地"   "構造物" "構造物" "構造物" "構造物"
[15] "植生"   "構造物" "構造物" "構造物" "構造物" "構造物" "構造物"
[22] "植生"   "構造物" "植生"   "植生"   "構造物" "構造物" "構造物"
[29] "構造物" "構造物" "構造物" "構造物" "裸地"   "構造物" "裸地"  
[36] "構造物" "裸地"   "裸地"   "構造物" "構造物" "構造物" "構造物"
[43] "構造物" "構造物" "構造物" "構造物" "構造物" "構造物" "構造物"
[50] "構造物" "構造物" "植生"   "植生"   "植生"   "植生"   "構造物"
[57] "植生"   "植生"   "構造物" "植生"   "植生"   "植生"   "植生"  
[64] "構造物" "裸地"   "裸地"   "裸地"   "裸地"   "植生"   "植生"  
[71] "その他" "植生"   "植生"   "植生"   "植生"  
 [ reached getOption("max.print") -- omitted 82 entries ]

2.5.4 属性値でデータを抽出したい

landcover_A<-landcover %>% dplyr::filter(class=="構造物")
plot(landcover_A["class"])

2.5.5 複数の属性値でデータを抽出したい

landcover_B<-landcover %>% dplyr::filter(subclass %in% c("高木","川"))
plot(landcover_B["subclass"])

2.5.6 各ポリゴンに重なるポイントのインデックスを知りたい

st_intersects(mesh,pnt)
Sparse geometry binary predicate list of length 50, where the predicate was `intersects'
first 10 elements:
 1: 363, 386, 438, 844
 2: 223, 325, 741, 822
 3: 131, 175, 524, 834
 4: 41, 753, 831
 5: 13, 132, 137
 6: 45, 311, 320, 681
 7: 174, 560, 963
 8: 262
 9: 410, 433, 496, 869
 10: 42, 86, 198, 546, 939

2.5.7 各ポリゴンに重なるポイントの数を知りたい

st_intersects(mesh,pnt) %>% lengths
 [1] 4 4 4 3 3 4 3 1 4 5 3 4 6 3 4 2 5 1 3 8 4 4 3 3 2 3 3 4 2 8 2 0 5 3 5
[36] 8 3 1 3 1 7 4 6 2 7 3 7 2 8 7

2.5.8 交差するオブジェクトを抽出したい

inter<-landcover %>% filter(lengths(st_intersects(.,mesh))>0)
plot(inter)

2.6 空間演算

2.6.1 交差したオブジェクトを切り抜きたい

inter<-st_intersection(landcover,mesh)
plot(inter)

2.6.2 ポリゴンの重心を求めたい

cent<-st_centroid(mesh)
plot(cent)

2.6.3 バッファーを作成したい

st_buffer(cent,2) %>% plot

2.7 データの計測

2.7.1 面積を計測したい

st_area(landcover)
Units: m^2
 [1]    1.2664181    1.1670632   92.4662731    4.7533855  168.4321983
 [6]   20.5295716   88.7457455 1166.3521519   10.6482112   14.3381424
[11]   18.4315028    8.1846281    0.2911716    8.0071366   14.2558465
[16]    6.3188338    4.5953293    5.9931429    8.5189867  501.9557047
[21]    3.0543304    0.8345024  188.8887054  211.9038885  106.0645716
[26]    1.9500176  140.1159740    0.6037353    2.0023096    1.5920548
[31]  137.7566426    0.4645672   34.6844957    0.7316131   43.8001451
[36]    6.7152484   18.2064899   42.4655335    0.5142420   49.6365397
[41]    1.0387100    6.5098833    0.5390834   12.4345711   11.8929091
[46]    5.3195393    3.1879128    3.1124514    3.0177756    3.6326343
[51]    3.8782796    7.4885220   31.9022097   27.0036081   34.7827617
[56]  987.6822113   11.7822736   62.4633140   61.5095772   62.6857929
[61]   63.1763873   62.3910560 1476.8087810   73.5310415   62.1990016
[66]   61.4193748   63.8913620   62.4633140  267.4098246    6.7647361
[71]    0.6137535    6.0617407   62.3549269   62.3891544  255.5535664
 [ reached getOption("max.print") -- omitted 82 entries ]

2.7.2 長さを計測したい

st_length(landcover_l)
Units: m
 [1]   4.509671   4.326394  62.415441  15.704905  74.471046  21.671120
 [7]  77.541357  17.711155   8.133101   8.139008 307.324447   4.509671
[13]   4.326394  13.516729  15.355734  34.242881   7.286367   6.493763
[19]   6.755525  14.739732  13.516729  17.711155  11.448565   2.117576
[25]  26.772089  15.769340  10.323347  10.526287  10.116995  11.918046
[31] 380.063402  15.769340  10.116995  61.600162   8.245661   8.141681
[37]  22.699572  24.822336  29.206505   7.748644   4.963406  60.729771
[43]  53.151124  59.544347   3.264167   3.168179  15.055240   3.179918
[49]   7.313497  61.600162   3.120614   5.678371   5.048299  47.953015
[55]   2.736776  46.131944   7.512688   3.429690  46.129084   3.429690
[61]  11.600841   2.938996  10.039243  12.055736  15.355734  23.410794
[67]   2.951162  34.242881   9.313325   4.082485  11.600841   2.938996
[73]  14.970448  14.777440   9.313325
 [ reached getOption("max.print") -- omitted 155 entries ]

2.7.3 地物の距離を計測したい

st_distance(cent[1,],cent[2:4,])
Units: m
     [,1] [,2] [,3]
[1,]    6   12   18

2.7.4 地物の座標を知りたい

st_coordinates(cent)
       X    Y
1  24923 2743
2  24929 2743
3  24935 2743
4  24941 2743
5  24947 2743
6  24953 2743
7  24959 2743
8  24965 2743
9  24971 2743
10 24977 2743
11 24923 2749
12 24929 2749
13 24935 2749
14 24941 2749
15 24947 2749
16 24953 2749
17 24959 2749
18 24965 2749
19 24971 2749
20 24977 2749
21 24923 2755
22 24929 2755
23 24935 2755
24 24941 2755
25 24947 2755
26 24953 2755
27 24959 2755
28 24965 2755
29 24971 2755
30 24977 2755
31 24923 2761
32 24929 2761
33 24935 2761
34 24941 2761
35 24947 2761
36 24953 2761
37 24959 2761
 [ reached getOption("max.print") -- 13 行を無視しました ]

2.8 属性情報の操作

2.8.1 属性にIDを振りたい

mesh<-mesh %>% mutate(ID=1:nrow(.))
mesh
Simple feature collection with 50 features and 1 field
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 24920 ymin: 2740 xmax: 24980 ymax: 2770
epsg (SRID):    2451
proj4string:    +proj=tmerc +lat_0=36 +lon_0=139.8333333333333 +k=0.9999 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
First 10 features:
   ID                           mesh
1   1 POLYGON ((24920 2740, 24926...
2   2 POLYGON ((24926 2740, 24932...
3   3 POLYGON ((24932 2740, 24938...
4   4 POLYGON ((24938 2740, 24944...
5   5 POLYGON ((24944 2740, 24950...
6   6 POLYGON ((24950 2740, 24956...
7   7 POLYGON ((24956 2740, 24962...
8   8 POLYGON ((24962 2740, 24968...
9   9 POLYGON ((24968 2740, 24974...
10 10 POLYGON ((24974 2740, 24980...

2.8.2 属性にポイントの座標を追加したい

参照:ポリゴンの重心を求めたい

cent_xy<-cbind(cent,st_coordinates(cent))

2.8.3 面積を属性に追加したい

landcover_area<-landcover %>% mutate(AREA=st_area(.))

2.8.4 面積の単位を変更したい

#ヘクタールha、アールa、平方キロ,km^2
landcover_area %>% mutate(AREA=units::set_units(.$AREA,ha))
Simple feature collection with 157 features and 3 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 24904.97 ymin: 2710.929 xmax: 24992.02 ymax: 2824.005
epsg (SRID):    NA
proj4string:    +proj=tmerc +lat_0=36 +lon_0=139.8333333333333 +k=0.9999 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs
# A tibble: 157 x 4
   class  subclass     AREA                                       geometry
   <chr>  <chr>        <S3: units>                           <POLYGON [m]>
 1 構造物 コンクリート 0.000126641806965745 ((24985.42 2713.855, 24986.65~
 2 構造物 コンクリート 0.00011670631985744  ((24933.71 2711.357, 24932.91~
 3 植生   芝           0.00924662730632117  ((24921.19 2711.823, 24920.89~
 4 植生   高木         0.000475338547215705 ((24917.46 2711.022, 24917.61~
 5 構造物 コンクリート 0.0168432198342496   ((24925.27 2710.929, 24925.06~
 6 植生   高木         0.00205295716293094  ((24917.46 2711.022, 24910.18~
 7 裸地   裸地         0.00887457454996586  ((24991.99 2718.085, 24991.99~
 8 植生   芝           0.116635215189833    ((24991.99 2717.693, 24991.99~
 9 植生   高木         0.00106482112367005  ((24905.07 2714.537, 24905.01~
10 裸地   裸地         0.00143381424049762  ((24970.5 2719.856, 24970.72 ~
# ... with 147 more rows

2.8.5 メッシュ内のポイント数を属性に追加したい

cnt<-lengths(st_intersects(mesh,pnt))
mesh2<-mesh %>% mutate(count=cnt)

2.8.6 属性値で計算したい

mesh2<-mesh2 %>% mutate(AREA=st_area(.))
mesh2 %>% mutate(density=count/AREA)
Simple feature collection with 50 features and 4 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 24920 ymin: 2740 xmax: 24980 ymax: 2770
epsg (SRID):    2451
proj4string:    +proj=tmerc +lat_0=36 +lon_0=139.8333333333333 +k=0.9999 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
First 10 features:
   ID count   AREA          density                           mesh
1   1     4 36 m^2 0.11111111 1/m^2 POLYGON ((24920 2740, 24926...
2   2     4 36 m^2 0.11111111 1/m^2 POLYGON ((24926 2740, 24932...
3   3     4 36 m^2 0.11111111 1/m^2 POLYGON ((24932 2740, 24938...
4   4     3 36 m^2 0.08333333 1/m^2 POLYGON ((24938 2740, 24944...
5   5     3 36 m^2 0.08333333 1/m^2 POLYGON ((24944 2740, 24950...
6   6     4 36 m^2 0.11111111 1/m^2 POLYGON ((24950 2740, 24956...
7   7     3 36 m^2 0.08333333 1/m^2 POLYGON ((24956 2740, 24962...
8   8     1 36 m^2 0.02777778 1/m^2 POLYGON ((24962 2740, 24968...
9   9     4 36 m^2 0.11111111 1/m^2 POLYGON ((24968 2740, 24974...
10 10     5 36 m^2 0.13888889 1/m^2 POLYGON ((24974 2740, 24980...

2.8.7 属性列を削除したい

landcover %>% dplyr::select(-subclass)
Simple feature collection with 157 features and 1 field
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 24904.97 ymin: 2710.929 xmax: 24992.02 ymax: 2824.005
epsg (SRID):    NA
proj4string:    +proj=tmerc +lat_0=36 +lon_0=139.8333333333333 +k=0.9999 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs
# A tibble: 157 x 2
   class                                                          geometry
   <chr>                                                     <POLYGON [m]>
 1 構造物 ((24985.42 2713.855, 24986.65 2713.886, 24986.65 2712.806, 2498~
 2 構造物 ((24933.71 2711.357, 24932.91 2712.066, 24933.68 2712.868, 2493~
 3 植生   ((24921.19 2711.823, 24920.89 2711.823, 24920.48 2711.781, 2492~
 4 植生   ((24917.46 2711.022, 24917.61 2711.126, 24917.85 2711.195, 2491~
 5 構造物 ((24925.27 2710.929, 24925.06 2710.985, 24924.73 2711.118, 2492~
 6 植生   ((24917.46 2711.022, 24910.18 2710.995, 24910.08 2711.091, 2490~
 7 裸地   ((24991.99 2718.085, 24991.99 2717.693, 24990.16 2717.663, 2498~
 8 植生   ((24991.99 2717.693, 24991.99 2710.986, 24931.58 2710.985, 2491~
 9 植生   ((24905.07 2714.537, 24905.01 2720.457, 24905.18 2720.342, 2490~
10 裸地   ((24970.5 2719.856, 24970.72 2719.979, 24971.08 2719.979, 24971~
# ... with 147 more rows

2.8.8 属性だけを取り出したい

landcover %>% st_set_geometry(NULL)
# A tibble: 157 x 2
   class  subclass    
 * <chr>  <chr>       
 1 構造物 コンクリート
 2 構造物 コンクリート
 3 植生   芝          
 4 植生   高木        
 5 構造物 コンクリート
 6 植生   高木        
 7 裸地   裸地        
 8 植生   芝          
 9 植生   高木        
10 裸地   裸地        
# ... with 147 more rows

2.8.9 同じ属性のポリゴンをマルチポリゴンにしたい

landcover %>% group_by(class) %>% summarise()
Simple feature collection with 5 features and 1 field
geometry type:  GEOMETRY
dimension:      XY
bbox:           xmin: 24904.97 ymin: 2710.929 xmax: 24992.02 ymax: 2824.005
epsg (SRID):    NA
proj4string:    +proj=tmerc +lat_0=36 +lon_0=139.8333333333333 +k=0.9999 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs
# A tibble: 5 x 2
  class                                                           geometry
  <chr>                                                     <GEOMETRY [m]>
1 その他 MULTIPOLYGON (((24932.53 2751.01, 24934.15 2751.134, 24934.27 27~
2 構造物 MULTIPOLYGON (((24985.42 2713.855, 24986.65 2713.886, 24986.65 2~
3 植生   MULTIPOLYGON (((24964.29 2728.878, 24963.5 2729.035, 24962.74 27~
4 水     POLYGON ((24921.17 2823.986, 24922.01 2823.986, 24921.54 2822.84~
5 裸地   MULTIPOLYGON (((24949.95 2723.714, 24950.26 2723.849, 24950.51 2~

2.8.10 同じ属性の面積を集計したい

landcover_area %>% group_by(class) %>% summarise(TOTALAREA=sum(AREA))
Simple feature collection with 5 features and 2 fields
geometry type:  GEOMETRY
dimension:      XY
bbox:           xmin: 24904.97 ymin: 2710.929 xmax: 24992.02 ymax: 2824.005
epsg (SRID):    NA
proj4string:    +proj=tmerc +lat_0=36 +lon_0=139.8333333333333 +k=0.9999 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs
# A tibble: 5 x 3
  class  TOTALAREA                                                geometry
  <chr>  <S3: units>                                        <GEOMETRY [m]>
1 その他 24.6008204800984 MULTIPOLYGON (((24932.53 2751.01, 24934.15 2751~
2 構造物 3191.21992748005 MULTIPOLYGON (((24985.42 2713.855, 24986.65 271~
3 植生   5891.81899363423 MULTIPOLYGON (((24964.29 2728.878, 24963.5 2729~
4 水     56.9571404562651 POLYGON ((24921.17 2823.986, 24922.01 2823.986,~
5 裸地   661.563330774753 MULTIPOLYGON (((24949.95 2723.714, 24950.26 272~

2.8.11 メッシュ内の同じ属性の面積を集計したい

landcover_mesh<- st_intersection(landcover,mesh) %>% mutate(ID=1:nrow(.),AREA=st_area(.))
landcover_mesh<-landcover_mesh %>% group_by(ID,subclass) %>% summarise(AREA=sum(AREA))

2.8.12 メッシュ内の最大面積の属性をメッシュの属性にしたい

df<-landcover_mesh %>% group_by(ID) %>% filter(AREA==max(AREA)) %>% st_set_geometry(NULL)
mesh %>% left_join(df,by="ID")
Simple feature collection with 50 features and 3 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 24920 ymin: 2740 xmax: 24980 ymax: 2770
epsg (SRID):    2451
proj4string:    +proj=tmerc +lat_0=36 +lon_0=139.8333333333333 +k=0.9999 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
First 10 features:
   ID     subclass          AREA                           mesh
1   1 コンクリート 16.765318 m^2 POLYGON ((24920 2740, 24926...
2   2         建物 10.795827 m^2 POLYGON ((24926 2740, 24932...
3   3         低木  2.600732 m^2 POLYGON ((24932 2740, 24938...
4   4           芝  5.838123 m^2 POLYGON ((24938 2740, 24944...
5   5 コンクリート 13.523171 m^2 POLYGON ((24944 2740, 24950...
6   6         建物 14.234150 m^2 POLYGON ((24950 2740, 24956...
7   7         低木  2.345148 m^2 POLYGON ((24956 2740, 24962...
8   8         低木  1.298345 m^2 POLYGON ((24962 2740, 24968...
9   9         低木  0.110244 m^2 POLYGON ((24968 2740, 24974...
10 10           芝  4.488941 m^2 POLYGON ((24974 2740, 24980...

3 ラスタ基本編

3.1 基本操作

3.1.1 ラスタデータを読み込みたい

ortho<-stack('data/rededge.tif')

3.1.2 ラスタデータをバンドを指定して読み込みたい

ortho1<-raster('data/rededge.tif',layer=1)

3.1.3 データの概要を確認したい

ortho
class       : RasterStack 
dimensions  : 2128, 1703, 3623984, 6  (nrow, ncol, ncell, nlayers)
resolution  : 0.1, 0.1  (x, y)
extent      : 24842.63, 25012.93, 2661.646, 2874.446  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=tmerc +lat_0=36 +lon_0=139.833333333333 +k=0.9999 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
names       : rededge.1, rededge.2, rededge.3, rededge.4, rededge.5, rededge.6 
min values  :    0.0000,    0.0000,    0.0000,    0.0000,  498.1926,    0.0000 
max values  :     65535,     65535,     65535,     65535,     65535,     65535

3.1.4 地図を表示したい

plot(ortho1)

3.1.5 総セル数を確認したい

ncell(ortho)
[1] 3623984

3.1.6 Xセル数、Yセル数、バンド数を確認したい

dim(ortho)
[1] 2128 1703    6

3.1.7 解像度を確認したい

res(ortho)
[1] 0.1 0.1

3.1.8 バンド数を確認したい

nlayers(ortho)
[1] 6

3.1.9 バンド名を確認したい

names(ortho)
[1] "rededge.1" "rededge.2" "rededge.3" "rededge.4" "rededge.5" "rededge.6"

3.1.10 セル値のサマリーを確認したい

summary(ortho1)
          rededge
Min.        0.000
1st Qu.  1751.245
Median   3709.821
3rd Qu. 65535.000
Max.    65535.000
NA's        0.000

3.1.11 データの範囲を確認したい

extent(ortho)
class       : Extent 
xmin        : 24842.63 
xmax        : 25012.93 
ymin        : 2661.646 
ymax        : 2874.446

3.1.12 ラスタデータを書き出したい

writeRaster(ortho1,filename = "mydata/ortho1.tif", format = "GTiff", overwrite = TRUE)

3.2 座標系の操作

3.2.1 データの空間参照系を確認したい

crs(ortho)
CRS arguments:
 +proj=tmerc +lat_0=36 +lon_0=139.833333333333 +k=0.9999 +x_0=0
+y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

3.2.2 空間参照系をEPSGコードで指定したい

crs(ortho)<-CRS("+init=epsg:2451")

3.2.3 別の座標系に変換したい

ortho_utm54<-projectRaster(ortho1, crs=CRS("+init=epsg:3100"))

3.3 ラスタの編集

3.3.1 ラスタースタックからレイヤを取り出したい

#subsetで
ortho1<-subset(ortho,1)
#rasterで
ortho1<-raster(ortho,1)
#リストで
ortho1<-ortho[[1]]

3.3.2 指定範囲で切り抜きたい

e<-extent(24905,24992,2711,2824)
orthocrop<-crop(ortho,e)

3.3.3 マウスで指定範囲を切り抜きたい

plot(ortho1)
orthocrop<-raster::select(ortho)
plot(orthocrop)

3.3.4 バンド名を指定したい

names(ortho) <- c('blue','green','red','nir','rededge','alpha')
ortho
class       : RasterStack 
dimensions  : 2128, 1703, 3623984, 6  (nrow, ncol, ncell, nlayers)
resolution  : 0.1, 0.1  (x, y)
extent      : 24842.63, 25012.93, 2661.646, 2874.446  (xmin, xmax, ymin, ymax)
coord. ref. : +init=epsg:2451 +proj=tmerc +lat_0=36 +lon_0=139.8333333333333 +k=0.9999 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
names       :     blue,    green,      red,      nir,  rededge,    alpha 
min values  :   0.0000,   0.0000,   0.0000,   0.0000, 498.1926,   0.0000 
max values  :    65535,    65535,    65535,    65535,    65535,    65535

3.3.5 解像度を変更したい

ortho_res5 <- ortho[[1]]
res(ortho_res5)<-5
ortho_res5 <- resample(ortho[[1]], ortho_res5, method='bilinear')
plot(ortho_res5)

3.4 値の取得・集計

ヒストグラムを表示したい

hist(ortho1,breaks = 30)

3.4.1 ラスタ値を取得したい

getValues(ortho1)
 [1] 65535 65535 65535 65535 65535 65535 65535 65535 65535 65535 65535
[12] 65535 65535 65535 65535 65535 65535 65535 65535 65535 65535 65535
[23] 65535 65535 65535 65535 65535 65535 65535 65535 65535 65535 65535
[34] 65535 65535 65535 65535 65535 65535 65535 65535 65535 65535 65535
[45] 65535 65535 65535 65535 65535 65535 65535 65535 65535 65535 65535
[56] 65535 65535 65535 65535 65535 65535 65535 65535 65535 65535 65535
[67] 65535 65535 65535 65535 65535 65535 65535 65535 65535
 [ reached getOption("max.print") -- omitted 3623909 entries ]

3.4.2 ラスタの値を座標付きで取得したい

df<-as.data.frame(ortho,xy=TRUE)
head(df)
         x        y  blue green   red   nir rededge alpha
1 24842.68 2874.396 65535 65535 65535 65535   65535     0
2 24842.78 2874.396 65535 65535 65535 65535   65535     0
3 24842.88 2874.396 65535 65535 65535 65535   65535     0
4 24842.98 2874.396 65535 65535 65535 65535   65535     0
5 24843.08 2874.396 65535 65535 65535 65535   65535     0
6 24843.18 2874.396 65535 65535 65535 65535   65535     0

3.4.3 ラスタの値をポイントで抽出したい

参照:ランダムポイントを作成したい

df<-raster::extract(ortho1,as(pnt,"Spatial"),df=TRUE)

3.5 ラスタ演算

3.5.1 ラスタ値を変更したい

ortho2<-ortho1/65536

3.5.2 ラスタ値に関数を適用したい

myfunc<-function(x){
  y=x/65536
  return(y)
}
ortho2 <- calc(ortho1, myfunc)

ラスタを再分類したい

ortho3 <- reclassify(ortho2, c(-Inf,0.25,1, 0.25,0.3,2, 0.3,0.4,3, 0.4,0.5,4, 0.5,Inf, 5))

3.6 地図表示

3.6.1 バンドと色を指定して表示したい

plot(ortho,1,col=gray.colors(30))

3.6.2 バンドをRGBに指定して表示したい

plotRGB(ortho, r = 3, g = 2, b = 1, axes = TRUE, stretch = "hist", main = "True Color Composite")

3.6.3 ズームして表示したい

#ズーム範囲の左上と右下をクリック
zoom(ortho)

3.6.4 ggplot2で表示したい

df<-as.data.frame(ortho1,xy=TRUE)
colnames(df)<-c("x","y","value")
ggplot()+
  geom_raster(data=df, aes(x = x, y = y,fill = value))+
  scale_fill_gradientn(colors=gray.colors(30))

4 応用編

4.1 マルチスペクトル画像から水稲の活性状況を把握する

4.1.1 データを読み込み、情報を確認する

#読み込み
sequoia <- stack('data/sequoia.tif')
#情報確認
crs(sequoia)
CRS arguments:
 +proj=tmerc +lat_0=36 +lon_0=139.833333333333 +k=0.9999 +x_0=0
+y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
ncell(sequoia)
[1] 3918123
dim(sequoia)
[1] 2083 1881    5
res(sequoia)
[1] 0.05 0.05
nlayers(sequoia)
[1] 5
sequoia
class       : RasterStack 
dimensions  : 2083, 1881, 3918123, 5  (nrow, ncol, ncell, nlayers)
resolution  : 0.05, 0.05  (x, y)
extent      : 24732.22, 24826.27, 2702.225, 2806.375  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=tmerc +lat_0=36 +lon_0=139.833333333333 +k=0.9999 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
names       : sequoia.1, sequoia.2, sequoia.3, sequoia.4, sequoia.5 
min values  :  165.1592,    0.0000,  267.5177,  441.6390,    0.0000 
max values  :     65535,     65535,     65535,     65535,     65535 
#表示
plot(sequoia,col=gray.colors(30))

4.1.2 データを前処理する

#処理するバンドを抽出
sequoia <- subset(sequoia,1:4)
#バンド前を付ける
names(sequoia) <- c('green','red','nir','rededge')
#空間参照系を指定
crs(sequoia)<-CRS('+init=EPSG:2451')
#確認
sequoia
class       : RasterStack 
dimensions  : 2083, 1881, 3918123, 4  (nrow, ncol, ncell, nlayers)
resolution  : 0.05, 0.05  (x, y)
extent      : 24732.22, 24826.27, 2702.225, 2806.375  (xmin, xmax, ymin, ymax)
coord. ref. : +init=EPSG:2451 +proj=tmerc +lat_0=36 +lon_0=139.8333333333333 +k=0.9999 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
names       :    green,      red,      nir,  rededge 
min values  : 165.1592,   0.0000, 267.5177, 441.6390 
max values  :    65535,    65535,    65535,    65535 
#処理範囲のshpを読み込む
tanbo<-read_sf("data/tanbo.shp",options=c("ENCODING=CP932"))
plot(tanbo)

#処理範囲でデータを切り抜く
croped<-crop(sequoia,tanbo)
plot(croped)

#処理範囲をマスクする
masked<-mask(croped,tanbo)
#表示して確認
plot(masked)

4.1.3 植生指標(NDVI)を計算する

#NDVIを計算する関数
calc_NDVI <- function(img, n, r) {
  red <- img[[r]]
  nir <- img[[n]]
  return ((nir-red) / (red + nir))
}

#NDVIを計算する。SEQUOIAバンドNIR = 3, red = 2.
ndvi <- calc_NDVI(masked,3, 2)
plot(ndvi, col = rev(terrain.colors(10)), main = 'SEQUOIA-NDVI')

#ズームして確認する
zoom(ndvi)

4.1.4 活性状況を色々な方法で把握する

4.1.4.1 NDVIのヒストグラムを表示してみる

hist(ndvi,
     main = "Distribution of NDVI values",
     xlab = "NDVI",
     ylab="Frequency",
     col = "wheat",
     xlim = c(-0.5, 1),
     breaks = 30,
     xaxt = 'n')
axis(side=1, at = seq(-0.5,1, 0.05), labels = seq(-0.5,1, 0.05))

4.1.4.2 活性状況が良くない場所を抽出してみる

veg <- calc(ndvi, function(x){x[x > 0.6] <- NA; return(x)})
plot(veg, main = 'Veg cover')

4.1.4.3 活性状況を5段階に分類してみる

vegc <- reclassify(ndvi, c(-Inf,0.4,1, 0.4,0.5,2, 0.5,0.6,3, 0.6,0.7,4, 0.7,Inf, 5))
plot(vegc,col = rev(terrain.colors(5)),breaks=0:5, main = 'NDVI based thresholding')

4.1.4.4 メッシュでNDVIを集計してみる

#メッシュを作成して回転させる
#cpを中心にラジアンrだけsfobjを回転する関数
rotate_sf<-function(sfobj,cp,r){
  rot <- function(a) matrix(c(cos(a), sin(a), -sin(a), cos(a)), 2, 2)
  rot_sf<-(sfobj-cp)*rot(-r)+cp
  return(rot_sf)
}

mesh <- st_make_grid(n=c(38,12), offset = c(24787.8,2755.4), cellsize = c(0.5,0.5))
cp<-st_sfc(st_point(c(24787.8,2755.4)))
mesh<-rotate_sf(mesh,cp,pi*50/180)
plot(mesh,axes=T)

#メッシュごとのNDVIの平均をデータフレームで取得
mesh_sp<-as(mesh,"Spatial") #spクラスに変換
df<-raster::extract(ndvi,mesh_sp,fun=mean,df=TRUE)

#メッシュの属性に結合
mesh<-st_sf(mesh) #属性が付けられるようにsfクラスに変換
mesh<-mesh %>% mutate(ndvi=df$layer)
plot(mesh)

#メッシュのndvi値を表示
mesh<-mesh %>% mutate(id=1:nrow(.)) #メッシュにidを追加
mesh_df<-mesh %>% st_set_geometry(NULL) #データフレームに変換
hist(mesh_df$ndvi) #ヒストグラム表示

mesh_df %>% filter(ndvi<0.5) #活性が良くないメッシュIDを表示
        ndvi  id
1  0.4837638   1
2  0.4989155  15
3  0.4506496  17
4  0.4701184  19
5  0.4569584  37
6  0.3599697  38
7  0.4209734 153
8  0.4473416 154
9  0.4399733 155
10 0.4507088 156
11 0.4219699 157
12 0.4665093 158
13 0.4513437 159
14 0.4452101 160
15 0.4568889 161
16 0.4784372 162
17 0.4753804 163
18 0.4800457 164
19 0.4631615 167
20 0.4377891 168
21 0.4620305 169
22 0.4412604 170
23 0.4013176 171
24 0.4142887 172
25 0.4396136 173
26 0.4770755 174
27 0.4804691 175
28 0.4859000 176
29 0.4158148 177
30 0.4467641 178
31 0.4888475 179
32 0.4697119 180
33 0.4632633 181
34 0.4702576 182
35 0.4532415 183
36 0.4476693 184
37 0.4439775 185
 [ reached getOption("max.print") -- 38 行を無視しました ] 
mesh %>% filter(ndvi<0.5) %>% dplyr::select(ndvi) %>% plot()

4.1.5 結果を書き出す

#ggplot2で背景画像とndviのメッシュを重ねる
sequoia_df<-as.data.frame(croped[[4]],xy=TRUE) #ラスタをデータフレームに変換
head(sequoia_df)
         x       y  rededge
1 24782.89 2774.45 25042.36
2 24782.94 2774.45 23915.62
3 24782.99 2774.45 22572.72
4 24783.04 2774.45 25547.49
5 24783.09 2774.45 28664.44
6 24783.14 2774.45 26880.76
#メッシュの中心座標を属性に追加
cent_xy<-mesh %>% st_centroid() %>% st_coordinates()
mesh<-cbind(mesh,cent_xy)
g<-ggplot() +
  geom_raster(data=sequoia_df, aes(x = x, y = y,alpha = rededge))+
  scale_alpha(name = "", range = c(0.6, 0), guide = F)+
  geom_sf(data = mesh %>% filter(ndvi>0.5),aes(fill = ndvi,alpha=0.3))+
  scale_fill_gradientn(colors=c("white","yellow","red"))+
  geom_text_repel(data = mesh %>% filter(id %in% c(1,38,191,343,456)),aes(x=X,y=Y,label = id),nudge_y=10)+
  annotate("text",x=24798,y=2756,label="数字は区画ID")
  
plot(g)

#画像をjpgに保存する
ggsave("mydata/ndvi.jpg",g)
#マスク処理したラスタデータを書き出す
writeRaster(masked,filename = "mydata/masked_sequoia.tif", format = "GTiff", overwrite = TRUE)
#ndviメッシュポリゴンを書き出す
st_write(mesh, "mydata/ndvi_mesh.shp",layer_options = "ENCODING=UTF-8",delete_layer = TRUE)
Deleting layer `ndvi_mesh' using driver `ESRI Shapefile'
Writing layer `ndvi_mesh' to data source `C:\Users\mizutani\Desktop\R_GIS_tutorial\mydata\ndvi_mesh.shp' using driver `ESRI Shapefile'
options:        ENCODING=UTF-8 
features:       456
fields:         4
geometry type:  Polygon

4.2 マルチスペクトル画像から土地被覆を分類する

4.2.1 前処理して、情報を確認する

rededge <- stack('data/rededge.tif')
rededge <- subset(rededge,1:5) #アルファバンドを除外
names(rededge) <- c('blue','green','red','nir','rededge')
crs(rededge)<-CRS("+init=epsg:2451") #CRSを設定

#解像度変更
rededge0.5<- rededge
res(rededge0.5)<-0.5
rededge <- resample(rededge, rededge0.5, method='bilinear')

#切り取り
e<-extent(24900,24995,2710,2825)
rededge<-crop(rededge,e)

#書き出し
writeRaster(rededge,filename = "mydata/cropped-rededge.tif", format = "GTiff", overwrite = TRUE)

#情報確認
crs(rededge)
CRS arguments:
 +init=epsg:2451 +proj=tmerc +lat_0=36 +lon_0=139.8333333333333
+k=0.9999 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0
+units=m +no_defs 
ncell(rededge)
[1] 43700
dim(rededge)
[1] 230 190   5
res(rededge)
[1] 0.5 0.5
nlayers(rededge)
[1] 5

4.2.2 色々な方法で表示する

blueバンド

plot(rededge,1,col=gray.colors(30))

TrueColor

plotRGB(rededge, r = 3, g = 2, b = 1, axes = TRUE, stretch = "lin", main = "RedEdge True Color Composite")

FalseColor

plotRGB(rededge, r = 4, g = 3, b = 2, axes = TRUE, stretch = "lin", main = "RedEdge False Color Composite")

NaturalColor

plotRGB(rededge, r = 3, g = 4, b = 2, axes = TRUE, stretch = "lin", main = "RedEdge Natural Color Composite")

4.2.3 データを分析する

4.2.3.1 バンド間の相関を見てみる

#各バンドを抽出し、16bitの値から0～1の反射率に変換
b_red<-rededge[[3]]/65536
b_nir<-rededge[[4]]/65536
b_rededge<-rededge[[5]]/65536

#相関グラフ
# 赤域と近赤外域
red_nir<-stack(b_red,b_nir)
pairs(red_nir, main = "Scatterplot between Red and NIR")

# 赤域とレッドエッジ域
red_rededge<-stack(b_red,b_rededge)
pairs(red_rededge, main = "Scatterplot between Red and Rededge")

4.2.3.2 バンドの反射率を土地被覆ごとにランダム点で取得してみる

#土地被覆の読み込み
landcover<-read_sf("data/landcover.shp",options=c("ENCODING=CP932"))
plot(landcover)

#サブクラスの一覧をベクトルで取得
subcl<-landcover %>% 
  st_set_geometry(NULL) %>% #sfからデータフレームに変換
  distinct(subclass) #重複を削除
subcl
# A tibble: 11 x 1
   subclass    
   <chr>       
 1 コンクリート
 2 芝          
 3 高木        
 4 裸地        
 5 建物        
 6 水田        
 7 砂地        
 8 その他      
 9 低木        
10 川          
11 草地        
#サブクラスごとにランダム点を500ポイント作成
subclass_sampling<-function(x){
  filter(landcover,subclass==x) %>% 
    st_sample(500)
}
pntlist<-subcl$subclass %>% map(subclass_sampling) 
pnt<-pntlist %>% reduce(c)

#ポイントの属性に土地被覆を追加
pnt<-st_intersection(landcover,pnt)
plot(pnt,pch=".",cex=3)

#各ポイントでバンド値を抽出
band_df <- raster::extract(rededge, as(pnt,"Spatial"),df=TRUE)
#ポイントの属性にバンド値を追加
pnt<-pnt %>% bind_cols(band_df)

4.2.3.3 地物とバンドの反射率の関係を見てみる

#土地被覆(subclass)ごとのバンドの平均値を計算
df_class_mean<-pnt %>% 
  st_set_geometry(NULL) %>% 
  dplyr::select(-ID,-class) %>% 
  group_by(subclass) %>% 
  summarise_all(mean) %>% 
  print
# A tibble: 11 x 6
   subclass      blue green   red   nir rededge
   <chr>        <dbl> <dbl> <dbl> <dbl>   <dbl>
 1 コンクリート 4200. 5487. 5994. 7491.  11450.
 2 その他       9196. 7750. 7135. 9521.  18707.
 3 建物         6568. 7649. 8010. 9054.  12310.
 4 高木         1110. 2867. 1151. 6334.  18198.
 5 砂地         5559. 6804. 7334. 8596.  11913.
 6 芝           2656. 4961. 4024. 8982.  17283.
 7 水田         2655. 3970. 4053. 6012.   9471.
 8 川           1302. 1934. 1317. 2513.   6398.
 9 草地         1027. 3005. 1146. 6910.  19473.
10 低木         1070. 2892. 1230. 6833.  19163.
11 裸地         3290. 5113. 5530. 8390.  14063.
#縦長なデータに変換
df_class_mean<-df_class_mean %>% 
  gather(key=Bands,value=Reflectance,-subclass) %>% 
  print
# A tibble: 55 x 3
   subclass     Bands Reflectance
   <chr>        <chr>       <dbl>
 1 コンクリート blue        4200.
 2 その他       blue        9196.
 3 建物         blue        6568.
 4 高木         blue        1110.
 5 砂地         blue        5559.
 6 芝           blue        2656.
 7 水田         blue        2655.
 8 川           blue        1302.
 9 草地         blue        1027.
10 低木         blue        1070.
# ... with 45 more rows
#反射率に変換
df_class_mean <- df_class_mean %>% mutate(Reflectance=Reflectance/65536)

#Bandsをfactorに変換（グラフ表示の際にBandsの順番を並べるため）
df_class_mean <- df_class_mean %>% mutate(Bands= factor(Bands, levels = c("blue", "green", "red", "rededge", "nir")))

# 各subclassのバンドの平均反射率をグラフ表示
# 結果を見るとレッドエッジの校正がおかしいかも（？）
classcolor <- c("コンクリート"="gray", "その他"="black","建物" = "white","高木"="darkgreen","砂地"="brown","芝"="lightgreen","水田"="lightblue", "川"="blue","草地"="orange","低木"="green","裸地"="yellow")
df_class_mean %>% 
  ggplot(aes(x=Bands ,y=Reflectance,col=subclass))+
  geom_line(aes(group = subclass),size=1)+
  geom_point()+
  labs(title = "Spectral Profile")+
  scale_color_manual(values=classcolor)

4.2.4 教師なし分類をする(5バンド使用)

set.seed(99)
km<-kmeans(values(rededge), centers = 11)
knr<-rededge[[1]]
knr[] <- km$cluster
par(mfrow = c(1,2))
plotRGB(rededge, r = 3, g = 2, b = 1, axes = TRUE, stretch = "lin", main = "RedEdge True Color Composite")
plot(knr, main = 'Unsupervised classification',breaks=0:11,col=rainbow(11))

par(mfrow = c(1,1))

4.2.5 教師あり分類をする

トレーニングデータを作成する

area1 = st_polygon(list(rbind(c(24910,2800), c(24930,2800), c(24930,2820), c(24910,2820), c(24910,2800))))
area2 = st_polygon(list(rbind(c(24950,2800), c(24970,2800), c(24970,2820), c(24950,2820), c(24950,2800))))
area3 = st_polygon(list(rbind(c(24950,2770), c(24970,2770), c(24970,2790), c(24950,2790), c(24950,2770))))
area4 = st_polygon(list(rbind(c(24910,2770), c(24930,2770), c(24930,2790), c(24910,2790), c(24910,2770))))

training_area<-st_sfc(area1,area2,area3,area4,crs=2451)

#トレーニングエリアの確認
plot(landcover %>% st_geometry(),axes=T)
plot(training_area,add=T)

#トレーニングエリアのデータを抽出
training_pnt<-st_intersection(pnt,training_area)
# データフレームに変換
training_df<-training_pnt %>% st_set_geometry(NULL) %>% dplyr::select(-class,-ID)
head(training_df)
# A tibble: 6 x 6
  subclass  blue green   red   nir rededge
  <chr>    <dbl> <dbl> <dbl> <dbl>   <dbl>
1 高木     1094. 3309. 1113. 8023.  22969.
2 高木      948. 2867.  949. 7083.  22065.
3 高木      674. 2093.  641. 4870.  15273.
4 高木      723. 2291.  701. 5412.  17460.
5 高木      865. 2710.  899. 6478.  19231.
6 高木      757. 2271.  788. 5483.  18185.

4.2.5.1 決定木で分類する

#モデルの作成
model <- rpart(as.factor(subclass)~., data=training_df, method="class")
plot(model, uniform=TRUE, main="Classification Tree",margin=0.05)
text(model, cex = 0.7)

#モデルから全域を予測
pr <- predict(rededge, model, type="class")

#予測結果をプロット
df<-as.data.frame(pr,xy=TRUE)

classcolor <- c("コンクリート"="gray", "その他"="black","建物" = "white","高木"="darkgreen","砂地"="brown","芝"="lightgreen","水田"="lightblue", "川"="blue","草地"="orange","低木"="green","裸地"="yellow")
ggplot()+
  geom_raster(data=df, aes(x = x, y = y,fill = layer_value))+
  scale_fill_manual(values=classcolor)+ coord_equal(ratio = 1) 

#予測結果を書き出す
writeRaster(pr,filename = "mydata/RP_predicted_landcover.tif", format = "GTiff", overwrite = TRUE)

4.2.5.2 ランダムフォレストで分類する

#モデルの作成
model <- randomForest(as.factor(subclass)~.,data=training_df)

#モデルから全域を予測
pr <- predict(rededge, model)

#予測結果をプロット
df<-as.data.frame(pr,xy=TRUE)
classcolor <- c("コンクリート"="gray", "その他"="black","建物" = "white","高木"="darkgreen","砂地"="brown","芝"="lightgreen","水田"="lightblue", "川"="blue","草地"="orange","低木"="green","裸地"="yellow")
ggplot()+
  geom_raster(data=df, aes(x = x, y = y,fill = layer_value))+
  scale_fill_manual(values=classcolor)+ coord_equal(ratio = 1) 

#予測結果を書き出す
writeRaster(pr,filename = "mydata/RF_predicted_landcover.tif", format = "GTiff", overwrite = TRUE)

4.2.5.3 SVMで分類する

#モデルの作成
model <- ksvm(as.factor(subclass)~.,data=training_df)

#モデルから全域を予測
pr <- predict(rededge, model)

#予測結果をプロット
classcolor <- c("コンクリート"="gray", "その他"="black","建物" = "white","高木"="darkgreen","砂地"="brown","芝"="lightgreen","水田"="lightblue", "川"="blue","草地"="orange","低木"="green","裸地"="yellow")
df<-as.data.frame(pr,xy=TRUE)
ggplot()+
  geom_raster(data=df, aes(x = x, y = y,fill = layer_value))+
  scale_fill_manual(values=classcolor)+ coord_equal(ratio = 1) 

#予測結果を書き出す
writeRaster(pr,filename = "mydata/SVM_predicted_landcover.tif", format = "GTiff", overwrite = TRUE)