はじめに

下記記事の続き
https://qiita.com/iwasaki_kenichi/items/3144e923b0a573eac6e1
https://qiita.com/iwasaki_kenichi/items/e2f0c6777263a93c7230

sfとラスタデータについて学び、レイヤーを分析する準備ができました。ようやく空間分析を勉強していきます。 sfとラスターデータの両方に、バッファリングやハルの計算などの単層の種類の分析や、交差、オーバーラップ、マスキング、クリッピングなどの多層操作を実行できる一連の機能を勉強していきます。

###バッファ層
バッファーの計算は重要な空間分析スキルです。バッファーは下記のようなものがあります。

学校から1キロメートル以内の道路の数を特定する
敏感な自然地域の近くにある有害廃棄物サイトの数を計算する

技術的には、投影されていない座標参照システムでデータをバッファリングできます。しかし、バッファ距離は重要な要素のため、投影されていないデータであれば、予め何らかのCRSに投影してバッファリングすることが推奨されています。

使用する関数名:`st_as_sf`,`st_transform()`,`st_buffer()`,`plotRGB()`

1.r

# Review df
df

# Convert the data frame to an sf object             
df_sf <- st_as_sf(df, coords = c("longitude", "latitude"), crs = 4326)

# Transform the points to match the manhattan CRS
df_crs <- st_transform(df_sf, crs = crs(manhattan, asText = TRUE))

# Buffer the points
df_buf <- st_buffer(df_crs, dist = 1000)

# Plot the manhattan image (it is multi-band)
plotRGB(manhattan)
plot(st_geometry(df_buf), col = "firebrick", add = TRUE)
plot(st_geometry(df_crs), pch = 16, add = TRUE)

以下、結果です。

> # Review df
> df
                      place longitude latitude
1     Empire State Building -73.98566 40.74844
2 Museum of Natural History -73.97398 40.78132
> 
> # Convert the data frame to an sf object
> df_sf <- st_as_sf(df, coords = c("longitude", "latitude"), crs = 4326)
> 
> # Transform the points to match the manhattan CRS
> df_crs <- st_transform(df_sf, crs = crs(manhattan, asText = TRUE))
> 
> # Buffer the points
> df_buf <- st_buffer(df_crs, dist = 1000)
> 
> # Plot the manhattan image (it is multi-band)
> plotRGB(manhattan)
> plot(st_geometry(df_buf), col = "firebrick", add = TRUE)
> plot(st_geometry(df_crs), pch = 16, add = TRUE)

バッファリングは、st_buffer()で実現することができました。それよりも、データをRに取り込み、適切な座標参照系を確保するほうが難しそうです。

###ポリゴンの重心計算
バッファリングと同様に、ポリゴンの重心の計算は、値の割り当てに使用され、マップのラベリングを支援するための基盤となる重要なジオプロセッシングタスクです。 sfでのこのための関数はst_centroid()です。

使用する関数名:`st_read`,`st_transform()`,`st_centroid()`,`st_geometry()`

2.r

# Read in the neighborhods shapefile
neighborhoods <- st_read("neighborhoods.shp")

# Project neighborhoods to match manhattan
neighborhoods_tf <- st_transform(neighborhoods, crs = 32618)

# Compute the neighborhood centroids
centroids <- st_centroid(neighborhoods_tf)

# Plot the neighborhood geometry
plot(st_geometry(neighborhoods_tf), col = "grey", border = "white")
plot(centroids, pch = 16, col = "firebrick", add = TRUE)

結果は下記の通り

多角形の重心の計算方法を知っていると、役立ちそうです。

##ベクトルデータの周囲に境界ボックスを作成する
sfを使用して、ベクトルデータの周囲の境界ボックスを計算できます。これらは、分析のために共通領域にレイヤーをクリップするためのポリゴンを作成したり、影響のある領域を特定したりするのに役立ちます。
sfには、バウンディングボックスの座標を抽出するための関数があります。st_bbox()です。これらの座標から新しいsfオブジェクト（ポリゴン）を作成します。これを行うには、sfのst_make_grid()を使用します。

使用する関数名:`st_geometry`,`st_bbox()`,`st_make_grid()`,`st_geometry()`

3.r

# Plot the neighborhoods and beech trees
plot(st_geometry(neighborhoods), col = "grey", border = "white")
plot(beech, add = TRUE, pch = 16, col = "forestgreen")

# Compute the coordinates of the bounding box
st_bbox(beech)

# Create a bounding box polygon
beech_box <- st_make_grid(beech, n = 1)

# Plot the neighborhoods, add the beech trees and add the new box
plot(st_geometry(neighborhoods), col = "grey", border = "white")
plot(beech, add = TRUE, pch = 16, col = "forestgreen")
plot(beech_box, add = TRUE)

結果は、下記の通りです。

> st_bbox(beech)
     xmin      ymin      xmax      ymax 
-74.25460  40.49894 -73.70078  40.91165

バウンディボックスのようなものは、GISソフト的なもので画像の縮尺が自由にできるものがないと、少し不便だなと感じました。少しわかりにくいですね。こういうときは、shinyとかを組み合わせると良いのでしょうか。それともshpを書き出しておいて、GISソフト等にて読み込むようにするのでしょうか。。。（後者はあまり利点がなさそうな。。。）

##分割された領域を1つの領域にまとめる
いわゆるディゾルブです。指定した属性に基づいて一つのシェープファイルに集約することをいいます。

使用する関数名:`st_buffer()`,`st_geometry()`,`st_union()`

4.r

# Buffer the beech trees by 3000
beech_buffer <- st_buffer(beech, 3000)

# Limit the object to just geometry
beech_buffers <- st_geometry(beech_buffer)

# Compute the number of features in beech_buffer
length(beech_buffers)

# Plot the tree buffers
plot(beech_buffers)

# Dissolve the buffers
beech_buf_union <- st_union(beech_buffers)

# Compute the number of features in beech_buf_union
length(beech_buf_union)

# Plot the dissolved buffers
plot(beech_buf_union)

結果は下記です。

ディゾルブすると、下記のようになります。

重なっているバッファが、一つの面に集約されていることが確認できました。

##分割された領域を1つの領域にまとめる
いわゆる凸包です。与えられた点群を包む最小の凸多角形のことです。

使用する関数名:`st_union()`,`st_convex_hull()`

5.r

# Look at the data frame to see the type of geometry
head(beech)

# Convert the points to a single multi-point
beech1 <- st_union(beech)

# Look at the data frame to see the type of geometry
head(beech1)

# Confirm that we went from 17 features to 1 feature
length(beech)
length(beech1)

# Compute the tight bounding box
beech_hull <- st_convex_hull(beech1)

# Plot the points together with the hull
plot(beech_hull, col = "red")
plot(beech1, add = TRUE)

以下は結果です。

plotの引数でadd = TRUEは、上書きで描写するオプションです。
凸包を計算する場合は、まずここのshpをユニオンして、集約する必要があることに注意が必要です。
headすると、MULTIPOINTとなることが確認できます。

##空間結合
多くの分析は、地理情報と地理情報を結合する必要があります。たとえば、各近傍にある樹木の数を知りたいが、樹木データに近傍属性がない場合、st_join（）を使用した空間結合を行います。
st_join（）関数は入力としてsfデータフレームを必要とし、sfジオメトリのみのオブジェクトを受け入れません。 st_sf（）関数を使用して、sfジオメトリオブジェクトをsfデータフレームに変換します。

使用する関数名:`st_join()`,`st_sf()`

6.r

# Plot the beech on top of the neighborhoods
plot(st_geometry(neighborhoods))
plot(beech, add = TRUE, pch = 16, col = "red")

# Determine whether beech has class data.frame
class(beech)

# Convert the beech geometry to a sf data frame
beech_df <- st_sf(beech)

# Confirm that beech now has the data.frame class
class(beech_df)

# Join the beech trees with the neighborhoods
beech_neigh <- st_join(beech_df, neighborhoods)

# Confirm that beech_neigh has the neighborhood information
head(beech_neigh)

結果は下記です。

最初のplotはneightborhoodsを表示し、2回目のplotでは、beechを表示させています。

class(beech)

では、"sfc_POINT" "sfc"と出力されました。
その後、st_sf()を適用したbeech_dfのclass()は"sf" "data.frame"でした。data.frame型に変換されたことがわかります。

その後、st_join()をすると

> head(beech_neigh)
Simple feature collection with 6 features and 5 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 569653 ymin: 4488294 xmax: 608444.9 ymax: 4514770
epsg (SRID):    32618
proj4string:    +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
  county_fip     boro_name ntacode                                    ntaname
1        081        Queens    QN06                 Jamaica Estates-Holliswood
2        085 Staten Island    SI25                      Oakwood-Oakwood Beach
3        081        Queens    QN43                                  Bellerose
4        061     Manhattan    MN28                            Lower East Side
5        081        Queens    QN72                                   Steinway
6        085 Staten Island    SI01 Annadale-Huguenot-Prince's Bay-Eltingville
  boro_code                 geometry
1         4 POINT (602156.8 4508959)
2         5 POINT (573777.4 4491199)
3         4 POINT (608444.9 4510249)
4         1 POINT (586276.7 4507258)
5         4   POINT (592365 4514770)
6         5   POINT (569653 4488294)

という結果になりました。
結合していることがわかります。

ちなみにhead(neightborhoods)を行うと、

> head(beech_neigh)
> head(neighborhoods)
Simple feature collection with 6 features and 5 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 583943.4 ymin: 4496263 xmax: 603361.8 ymax: 4521077
epsg (SRID):    32618
proj4string:    +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
  county_fip boro_name ntacode            ntaname boro_code
1        047  Brooklyn    BK88       Borough Park         3
2        081    Queens    QN52      East Flushing         4
3        081    Queens    QN48         Auburndale         4
4        081    Queens    QN51        Murray Hill         4
5        081    Queens    QN27      East Elmhurst         4
6        005     Bronx    BX35 Morrisania-Melrose         2
                        geometry
1 MULTIPOLYGON (((586594.6 44...
2 MULTIPOLYGON (((601719.2 45...
3 MULTIPOLYGON (((603361.8 45...
4 MULTIPOLYGON (((600944.6 45...
5 MULTIPOLYGON (((596125.7 45...
6 MULTIPOLYGON (((592999.6 45...

という結果になります。
また、head(beech)をすると

> head(beech)
POINT (586276.7 4507258)
Geometry set for 6 features 
geometry type:  POINT
dimension:      XY
bbox:           xmin: 569653 ymin: 4488294 xmax: 608444.9 ymax: 4514770
epsg (SRID):    32618
proj4string:    +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
First 5 geometries:
POINT (592365 4514770)
POINT (602156.8 4508959)
POINT (573777.4 4491199)
POINT (608444.9 4510249)
POINT (586276.7 4507258)
POINT (592365 4514770)

となります。

##バッファ内の物体を判定
st_intersects()です。複数のポリゴンを入力し、ポリゴンが交差するもののみ、物体があると判定します。

使用する関数名:`st_intersects()`,`st_contains()`,`st_intersection()`

7.r

# Identify neighborhoods that intersect with the buffer
neighborhoods_int <- st_intersects(buf, neighborhoods)

# Identify neighborhoods contained by the buffer
neighborhoods_cont <- st_contains(buf, neighborhoods)

# Get the indexes of which neighborhoods intersect
# and are contained by the buffer
int <- neighborhoods_int[[1]]
cont <- neighborhoods_cont[[1]]

# Get the names of the names of neighborhoods in buffer
neighborhoods$ntaname[int]

# Clip the neighborhood layer by the buffer (ignore the warning)
neighborhoods_clip <- st_intersection(buf, neighborhoods)

# Plot the geometry of the clipped neighborhoods
plot(st_geometry(neighborhoods_clip), col = "red")
plot(neighborhoods[cont,], add = TRUE, col = "yellow")

結果は下記の通り。

今回使用しているのは、st_intersects(),st_contains(),st_intersection()の3つです。どれも領域をクリップするために使う空間関数です。
このあたりは、本サイトがとても詳しいです。

##物体間の距離を測定
ポリゴン間の距離を測定することは、空間分析スキルとして重要です。距離を測定する関数はrgeoやgeosphereパッケージにもありますが、sfにもあります。st_distance()です。

エンパイアステートビルからすべての公園までの距離を測定し、最も近い公園を特定します。

使用する関数名:`st_intersects()`,`st_contains()`,`st_intersection()`

7.r

# Read in the parks object
parks <- st_read("parks.shp")
plot(st_geometry(parks))

# Test whether the CRS match
st_crs(empire_state) == st_crs(parks)

# Project parks to match empire state
parks_es <- st_transform(parks, crs = st_crs(empire_state))

st_crs(empire_state) == st_crs(parks_es)


# Compute the distance between empire_state and parks_es
d <- st_distance(empire_state, parks_es)

# Take a quick look at the result
head(d)

# Find the index of the nearest park
nearest <- which.min(d)

# Identify the park that is nearest
parks_es[nearest, ]

st_crs(empire_state) == st_crs(parks)で、CRSが同じかどうか判定していますが、ここでFALSEが出力されます。
そのため、st_transformによってCRSを同じにしています。
念の為、st_crs(empire_state) == st_crs(parks_es)と、再度検証していますが、これはTRUEと出力されます。

次に、

d <- st_distance(empire_state, parks_es)
head(d)

としています。
ここでの出力は、

Units: [m]
[1]  3055.791 19969.336 22737.903 13846.358  4604.069 18541.779

です。

最後の

# Find the index of the nearest park
nearest <- which.min(d)

# Identify the park that is nearest
parks_es[nearest, ]

で、一番近い公園を抽出しています。
結果は下記です。

Simple feature collection with 1 feature and 14 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 585392 ymin: 4511320 xmax: 585418.2 ymax: 4511379
epsg (SRID):    32618
proj4string:    +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
                                                      location communityb
188 Broadway, Av of the Americas, bet. W. 32 St. and W. 33 St.        105
    nys_senate            signname zipcode us_congres gispropnum borough
188         27 Greeley Square Park   10001         12       M032       M
    waterfront nys_assemb councildis acres     typecatego      address
188         No         75          3 0.144 Triangle/Plaza 894 6 AVENUE
                          geometry
188 MULTIPOLYGON (((585411 4511...

エンパイヤステートビルから一番近い公園はGreeley Square Parkであることがわかりました。
公園ごとに距離がでてくるため、最後の行を行い、抽出する必要はあります。

##ラスタデータをマスクする
マスクとクロップは同じような処理です。ラスタを特定の関心領域に切り抜くことができます。mask()を使うと、基本的には関心領域をラスタの上に置き、境界外のラスタセルにはNA値が割り当てられます。
ただし、現在rasterパッケージはsfオブジェクトをサポートシていないです。そのため、
as(input,"Spatial")により、Spatialオブジェクトに変換する必要があります。

使用する関数名:`st_transform()`,`st_area()`,`unclass()`,`mask()`

9.r

# Project parks to match canopy
parks_cp <- st_transform(parks, crs = st_crs(canopy, asText = TRUE))

# Compute the area of the parks
areas <- st_area(parks_cp)

# Filter to parks with areas > 30000
parks_big <- filter(parks_cp, unclass(areas) > 30000)

# Plot the canopy raster
plot(canopy)

# Plot the geometry of parks_big
plot(st_geometry(parks_big))

# Convert parks to a Spatial object
parks_sp <- as(parks_big, "Spatial")

# Mask the canopy layer with parks_sp and save as canopy_mask
canopy_mask <- mask(canopy, mask = parks_sp)

# Plot canopy_mask -- this is a raster!
plot(canopy_mask)

結果は下記です。

1枚目のplotはcanopyのラスタデータです。
2枚目のplotはparks_bigのポリゴンデータです。
3枚目のplotはcanopy_maskは、2枚目を使って、1枚目をマスクしたデータです。

parks_spは、parks_sp <- as(parks_big, "Spatial")によりSpatialに変換しています。

正直このあたりの処理はGISソフトを見ながら処理したほうがいいのではないかと感じました。処理を繰り返すようならこの方法が適していますが。

##ラスタデータをクロップする
mask()を使った場合、ラスタのエクステントは変更されません。入力ラスタのエクステントと、マスク自体のエクステントが異なっていてもmask()を実行した後には異なるエクステントになります。
しかし、多くの場合はラスタのエクステントを他のレイヤーと共有したい場合があります。crop() を使うと、ラスタのエクステント (外接枠) が入力されたクロップレイヤーのエクステントと一致するようにラスタをクロップします。しかし、外接枠内ではマスキングは行われません (ラスタセルは NA に設定されていません)。

大規模な公園に基づいて NYC のキャノピーレイヤーをマスクして切り抜き、比較してみましょう。マスクされたラスターには多くの NA 値が含まれており (空白があります)、元のキャノピーレイヤーと同じ範囲になっていることに気づくはずです。切り抜いたレイヤーでは、切り抜いたキャノピーレイヤーの範囲が大規模公園の範囲と一致していることに気づくはずです (基本的にはズームインされています)。

使用する関数名:`crop()`,`mask()`

9.r

# Convert the parks_big to a Spatial object
parks_sp <- as(parks_big, "Spatial")

# Mask the canopy with the large parks 
canopy_mask <- mask(canopy, mask = parks_sp)

# Plot the mask
plot(canopy_mask)

# Crop canopy with parks_sp
canopy_crop <- crop(canopy, parks_sp)

# Plot the cropped version and compare
plot(canopy_crop)

1枚目のplotはmask()したデータです。
2枚目のplotはcrop()したデータです。

マスクはエクステントが変更されません。マスク外のラスタ値はすべてNAになります。
クロップは、他のレイヤーに合わせて切り取られることがわかりました。

##ラスタデータのラスタ値
マスキングやトリミングを行うだけではなく、関心領域のセル値を知りたい場合もあります。例えば、ランドマークや大きな公園内のキャノピーの割合を知りたい場合などです。そこで便利なのが extract() 関数です。

便利なことに、このことは後の分析で説明しますが、抽出されたセルに適用される関数を extract() に与えることができます。例えば、extract()を使用して近隣別にラスタ値を抽出し、fun = mean引数を指定すると、近隣別の平均セル値を返します。

他のラスタ関数と同様に、まだ sf オブジェクトを受け入れるように設定されていないので、Spatial オブジェクトに変換する必要があります。

使用する関数名:`st_transform()`,`extract()`

10.r

# Project the landmarks to match canopy
landmarks_cp <- st_transform(landmarks, crs = st_crs(canopy, asText = TRUE))

# Convert the landmarks to a Spatial object
landmarks_sp <- as(landmarks_cp, "Spatial")

# Extract the canopy values at the landmarks
landmarks_ex <- extract(canopy, landmarks_sp)

# Look at the landmarks and extraction results
landmarks_cp
landmarks_ex

結果は下記です。

> landmarks_cp
Simple feature collection with 3 features and 1 field
geometry type:  POINT
dimension:      XY
bbox:           xmin: 1827046 ymin: 2183388 xmax: 1827608 ymax: 2187197
epsg (SRID):    NA
proj4string:    +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
                      place                    geom
1     Empire State Building POINT (1827046 2183388)
2 Museum of Natural History POINT (1827152 2187197)
3 Central Park (The Ramble) POINT (1827608 2186846)
> landmarks_ex
[1]  0.00 18.22 63.96

extract()は非常に便利なツールで、ポイントだけではなく、ポリゴンでも使用できることがわかります。

##オーバーレイを使ったラスター計算
同じ情報源である米国地質調査所のキャノピー層と「不浸透性」層を使用します。不浸透性は、水が表面を通過できるかどうかを測定します。つまり、不浸透率の高い表面は水を通さない道路であり、不浸透率の低い表面は芝生のようなものかもしれません。

この演習で行うことは、基本的には、樹木の樹冠の割合が低い（[除去]80％）地域を見つけることで、最も都市部の場所を特定することです。これを行うために、ラスター計算を行う関数 f を定義しました。

使用する関数名:`raster()`,`overlay()`

11.r

# Read in the canopy and impervious layer
canopy <- raster("canopy.tif")
impervious <- raster("impervious.tif")

# Function f with 2 arguments and the raster math code
f <- function(rast1, rast2) {
  rast1 < 20 & rast2 > 80
}

# Do the overlay using f as fun
canopy_imperv_overlay <- overlay(canopy, impervious, fun = f)

# Plot the result (low tree canopy and high impervious areas)
plot(canopy_imperv_overlay)

結果は以下です。

ラスター関数 overlay() を使ったラスター計算ができるようになりました。
20%未満の樹冠と80%以上の不浸透性のある地域に限定し、これらの地域はマンハッタンとブルックリンの一部を含む都市部の最も都市的な地域になります。

#おわりに
ベクターとラスターに対する演算方法を学ぶことができました。
ベクターでは、バッファリング、重心、凸包、ディゾルブ等
ラスターでは、マスク、クロップ、セルの値抽出、オーバーレイ等