Automating GIS Processes 2024 写経 Lesson 7 Processing and Analysis of Raster Data（１）

In the first part of our lesson, we reviewed the basics of raster files and explored how to read and work with both single-band and multiband raster files, as well as retrieve raster images via WMS. Now, in the second part of the lesson, we will focus on how to process and analyze raster data. Specifically, we will cover:

• Merging multiple raster files to create a raster mosaic
• Clipping rasters using a polygon
• Reclassifying raster data
• Performing slope analysis

このレッスン（レッスン７）の前半のパートにて、私たちはラスタファイルの基礎を振り返り、シングルバンドとマルチバンドのラスタファイルをどのように読み込み、扱う方法を探求しました。WMS経由で、ラスタのイメージを取得する方法も。
さて、この後半のパートでは、ラスタデータの処理と分析の仕方にフォーカスします。特に、以下を網羅します。

• 複数のラスタファイルを統合し、ひとつのラスタモザイクを作る
• ポリゴンで、ラスタデータをクリッピングする
• ラスタデータを再分類する
• 傾斜分析を行う

For this part of the lesson, we will continue using the elevation model data provided by the National Land Survey of Finland. This time, we will work with all four raster tiles (L4133A, L4133B, L4133C, L4133D), which cover parts of Helsinki’s city center. These tiles have already been downloaded to the data directory and will be used in this lesson.

Now let’s read all four raster tiles and visualize them using subplots.

このパートでは、継続して、National Land Survey of Finlandが提供している、標高モデルのデータを使います。
今回は、すべてのラスタタイル（L4133A, L4133B, L4133C, L4133D）を取り扱います。これらは、ヘルシンキ市街地の一部を網羅します。これらのタイルは、すでにdataディレクトリにダウンロードされていて、このレッスンで使われます。

さて、４つのラスタタイルを読み込み、サブプロットで可視化しましょう。

（ファイルはこちらにあります。）
https://github.com/Automating-GIS-processes/notebooks/tree/main/lesson-7/data

まずは、rioxarrayライブラリをインストールします。

!pip install rioxarray

サブプロットで可視化します。

import rioxarray
import matplotlib.pyplot as plt

# Load the rasters and reproject them to EPSG:3067
raster1 = rioxarray.open_rasterio('data/L4133A.tif').rio.reproject("EPSG:3067")
raster2 = rioxarray.open_rasterio('data/L4133B.tif').rio.reproject("EPSG:3067")
raster3 = rioxarray.open_rasterio('data/L4133C.tif').rio.reproject("EPSG:3067")
raster4 = rioxarray.open_rasterio('data/L4133D.tif').rio.reproject("EPSG:3067")

# Create a subplot with 2x2 grid
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# Plot each raster in a separate subplot
raster1.plot(ax=axes[0, 0], cmap='viridis')
axes[0, 0].set_title('Raster 1 (L4133A)')

raster2.plot(ax=axes[0, 1], cmap='viridis')
axes[0, 1].set_title('Raster 2 (L4133B)')

raster3.plot(ax=axes[1, 0], cmap='viridis')
axes[1, 0].set_title('Raster 3 (L4133C)')

raster4.plot(ax=axes[1, 1], cmap='viridis')
axes[1, 1].set_title('Raster 4 (L4133D)')

# Adjust layout
plt.tight_layout()
plt.show()

Note: We are again using matplotlib for visualization of our raster datasets. Everything we’ve previously learned about visualization with matplotlib and styling still applies here. We have extensively used matplotlib during the GeoPython Lesson 7 and throughout this course, especially in Lesson 5 where we focused on creating static maps.

メモ：私たちは、ラスタのデータセットの可視化のために、また、matplotlibを使っています。私たちが以前に学んだmatplotlibによる可視化とスタイルの整え方は、ここでも適用されます。 GeoPythonのレッスン７やこのコースで、特に、静的な地図を作ることに注力したレッスン５で、私たちは手広くmatplotlibを使っています。

As you can see, the four raster tiles cover adjacent areas in the Helsinki city center. However, they are still separate files! In the next step, we will merge these tiles into a single raster mosaic, creating a seamless representation that covers the entire central area of Helsinki.

ご覧の通り、４つのラスタのタイルは、ヘルシンキの市街地の隣接するエリアを網羅しています。
しかし、それらはまだ、別々のファイルです！次のステップとして、この４つのタイルをひとつのラスタモザイクに統合し、ヘルシンキの市街地エリア全体を網羅する継ぎ目のないものを作りましょう。

Create raster mosaic

Raster files are usualy large in size, therefor it is quite common for providers to publish them in smaller pieces or tiles. While this makes it easier to transfer the data, it may not be so practical when it comes to the actual analysis. For example, as seen above, the elevation model from center of Helsinki is divided into 4 separate raster files.

Now we want to merge multiple raster files together and create a raster mosaic. This can be done with the merge_arrays() -function in rioxarray.

ラスタファイルは一般的にはサイズが大きいです。
それゆえ、提供側はそれをより小さい一部やファイルで配信するのがかなり一般的です。
これはデータを転送することをより簡単にしますが、分析をする際にはたらそれほど実践的でないかもしれません。
例えば、上記のように、ヘルシンキの中心部の標高モデルは４つの別々のラスタファイルに分割されています。

さて、私たちは複数のラスタファイルを一つに統合し、１つのラスタモザイクを作りたいと思っています。
これは、rioxarrayライブラリのmerge_arrays()関数によって実現できます。

Aligning raster tiles In our case, the raster files are coming from the same source, and are actually pieces of a larger data. If that was not the case, it would have been a good idea to start by aligning our raster tiles, ensuring that they can be combined without alignment problems. The Align function helps align your raster files in terms of resolution, projection, and spatial extent (docs.xarray).

ラスタファイルを揃える
今回の場合は、ラスタファイルは同じソースから作成され、実際にはより大きなデータの一部です。
そうでないのならば、位置合わせの問題なしに結合されることが保証されているラスタファイルを揃えることから始めるのが良いかもしれません。
Align関数は、解像・投影・そして空間拡張について、ラスタファイルを揃えることに役立ちます。

import xarray as xr
# Ensure all rasters have the same CRS
assert raster1.rio.crs == raster2.rio.crs == raster3.rio.crs == raster4.rio.crs, "Rasters have different CRS"

# Align the rasters to the same resolution and spatial extent
raster1, raster2 = xr.align(raster1, raster2)
raster3, raster4 = xr.align(raster3, raster4)

Now let’s merge our data:

では、データを統合しましょう。（比較しやすいように、cmap（カラーマップ）は、統合前と同じviridisにしています。）

from rioxarray.merge import merge_arrays

# Merge the four raster data into one
mosaic_merged = merge_arrays([raster1, raster2, raster3, raster4])

# Save the merged raster (optional)
mosaic_merged.rio.to_raster('merged_raster.tif')

# Plot the mosaic raster
mosaic_merged.plot(cmap='viridis')

Clipping raster using a polygon

Clipping a raster using vector data (i.e., a polygon) is another common operation with raster data. In this part of lesson, we use wfs to get the administrative boundaries in Helsinki area.

ベクタデータ（例えばポリゴン）を使ってラスタデータをクリッピングすることもまた、ラスタデータの一般的な操作です。
ここでは、ヘルシンキの行政の境界を取得するために、Web Feature Service（wfs）を使います。

まずは、wfsを使うために必要な、owslibライブラリをインストールします。

! pip install owslib

wfsでヘルシンキの行政の境界を取得し、先頭の５件を表示します。

from owslib.wfs import WebFeatureService
import geopandas as gpd

# Define the bounding box for Helsinki (EPSG:3879)
bbox_helsinki = "25400000,6670000,25500000,6680000"

# Fetch the Paavo data limited to the Helsinki area using WFS
paavo_data_helsinki = gpd.read_file(
    "https://geo.stat.fi/geoserver/wfs"
    "?service=wfs"
    "&version=2.0.0"
    "&request=GetFeature"
    "&typeName=postialue:pno"  # Adjust to the correct layer if needed
    "&srsName=EPSG:3879"
    f"&bbox={bbox_helsinki},EPSG:3879"
)

# Display the first few rows of the data
paavo_data_helsinki.head()

# Save to a file if needed
# paavo_data_helsinki.to_file("paavo_data_helsinki.gpkg", driver="GPKG")

Now we can use the contents of the column posti_alue to select a ppostal area of our choice. Let’s try Kamppi; we know that the postal code for Kamppi is 00100.

これで、郵便番号ごとのエリア（ppostal area）を選択するために、posti_alue列の内容を使うことができます。
では、Kamppiエリアで試してみましょう。Kamppiエリアの郵便番号は、00100です。

# Select Kamppi polygon based on postal code
kamppi = paavo_data_helsinki[paavo_data_helsinki['posti_alue'] == '00100']

# Transfrom the polygon to the same CRS as our raster data
kamppi = kamppi.to_crs(mosaic_merged.rio.crs)

kamppi.plot()

Note: Similar to other map overlay analyses, it is essential to ensure that the CRS (Coordinate Reference System) of the raster files and the clipping feature are matching.

メモ 地図のオーバーレイ分析などと同じように、ラスタデータのCRSとクリッピングするフィーチャーのCRSが一致していることが肝要です。

Once again, let’s make sure that the CRS are matching and then let’s proceed with the clipping.

ではもう一度、CRSがあっていることを確認しましょう。そして、クリッピングの作業を進めましょう。

# Double check that the CRS match
assert raster1.rio.crs == kamppi.crs , "CRS Mismatch"

# Clip the raster using the Kamppi polygon
clipped_mosaic = mosaic_merged.rio.clip(kamppi.geometry, kamppi.crs)

# Save the clipped raster (optional)
clipped_mosaic.rio.to_raster("clipped_mosaic_kamppi.tif")

clipped_mosaic.plot()

これだけだと、なんだかよくわからいないですね。

clip()はクリップされなかったエリアの値をrio.nodataにします。そのため、相対的にrio.nodataが強調される結果となっていまう。

ひとまず、rio.nodataだけを見てみましょう。
なお、clipped_mosaic.rio.nodataが非数です。where()については、後述します。

clipped_mosaic.where(clipped_mosaic == clipped_mosaic.rio.nodata).plot()

緑色のエリアがrio.nodataです。

plot()は非数（NaN）を描画しないので、プロットの前に、clipped_mosaicのrio.nodataを非数（NaN）にマスクし、その後、プロットします。

nodata_value = clipped_mosaic.rio.nodata

# Mask the NoData values
clipped_raster = clipped_mosaic.where(clipped_mosaic != nodata_value)
clipped_raster.plot()

（クリップされなかったエリアは非数（NaN）で描画されないため）クリップした箇所の色の違いが際立つようになりました。

The function np.where in NumPy

The np.where function in NumPy is a conditional function that allows you to select elements from an array based on a condition. It can be thought of as an element-wise if-else for arrays.

NumPyライブラリのnp.where関数は条件関数です。条件の配列を使って、配列内の要素を選ぶことができます。
各要素のif-elseが配列になっている、と考えることもできます。

Basic Syntax

np.where(condition, x, y)

• condition: A boolean array or condition that’s checked for each element.
• x: The value (or array) to use when the condition is True.
• y: The value (or array) to use when the condition is False.

• 条件：真偽値の配列、もしくは、各要素に対してチェックされる条件
• X：条件が成立する場合の値（もしくは配列）
• y：条件が成立しない場合の値（もしくは配列）

For each element, if the condition is True, np.where returns the corresponding value from x; otherwise, it returns the corresponding value from y.

各々の要素について、条件が成立するのならば、np.where関数はxから対応する値を返します。条件が成立しないのならば、np.where関数はyから対応する値を返します。

Example
Here’s a simple example to clarify how it works:

以下はどのようにこの関数が機能するかのサンプルです。

import numpy as np

arr = np.array([1, -2, 3, -4, 5])
result = np.where(arr > 0, arr, 0)
print(result)

Output:

[1, 0, 3, 0, 5]

In this case:

• The condition arr > 0 checks if each element is positive.
• Where the condition is True (for positive numbers), np.where returns the original arr value.
• Where the condition is False (for negative numbers), it returns 0.

この場合では、

• 条件arr > 0は、配列の各々の要素がプラス値かをチェックします
• 各々の要素で、条件に該当する（値がプラス値）の場合、np.whereはarrの値を要素の値としてそのまま返します
• 各々の要素で、条件に該当しない（値がマイナス値）の場合（正確には0またはマイナス値の場合）、np.whereは0を返します

Practical Uses
np.where is often used for:

• Masking values in an array (e.g., setting invalid values to NaN).
• Conditional selection in data processing and filtering.
• Replacing values in an array based on certain conditions.

np.whereはよく以下の用途で使われます

• 配列内の値のマスキング（例：不正な値を非数（NaN）にする）
  上記のケース。where(clipped_mosaic != nodata_value)のように、第二・第三引数を省略した場合、「条件に当てはまらない要素（「nodata_valueでない」でない、つまり、nodata_valueの要素）は、すべて NaN (非数) に置き換わる。
• データ処理やフィルタリングにおける条件つきの選択
• 特定の条件に基づく、配列内の値の置き換え

続きます。