Automating GIS Processes 2024 写経 Lesson 7 Exploring raster data in Python（２）

前回の続きです。

Loading raster through WMS

Let’s take a quick look at an example of how we can retrieve raster files or map images via Web Map Service (WMS). WMS is a standard protocol developed by the Open Geospatial Consortium (OGC) that allows users to request and retrieve georeferenced map images over the internet, typically in common formats such as PNG, JPEG, and GeoTIFF. It’s widely used in GIS applications, serving as a common method for organizations, like the City of Helsinki in our case, to share raster data.

Web Map Service (WMS)からどのようにしてラスタファイルもしくは地図画像を取得するのか、事例を見てみましょう。
WMSは、Open Geospatial Consortium (OGC)によって定められた標準の規約で、これにより、使用者は、インターネットで、ジオリファレンスされた（地理情報が割り当てられた）PNG・JPEG・GeoTIFFといった一般的なフォーマットで、地図イメージを要求・取得ができます。
これはGISアプリケーションでは広く使われていて、例えば、ヘルシンキ市のような組織が、ラスタデータを共有すために、一般的な方法として提供しています。

In the example below, we fetch a base map of Helsinki from the Helsingin kantakartta dataset, which is provided by the Helsinki City Geographic Information Service through their WMS. This dataset includes detailed maps of the Helsinki area, containing information about land use, urban planning, and public infrastructure. We use OWSLib to retrieve the data and then read and visualize it using the rioxarray and rasterio libraries. The layer we work with, “avoindata:Kantakartta”, is in the WGS84 (EPSG:4326) coordinate system and covers part of the city.

以下、私たちは、Helsingin kantakartta dataset から、ヘルシンキのベースマップを取得します。これは、WMSを通じたヘルシンキ市の地理情報サービスによって提供されています。
このデータセットはヘルシンキの詳細な地図群で、土地使用、都市計画、公共インフラを含みます。
私たちは、データを取得するためにOWSLib を使い、rioxarray・rasterioライブラリを使って、読み込み、可視化します。
私たちが扱うレイヤーavoindata:Kantakartta（フィンランド語でオープンデータ：ベースマップという意味らしい）のCRSはWGS84 (EPSG:4326) で、ヘルシンキ市の一部を網羅します。

では実践です。
まずは、owslibをインポートします。これをしないと、from owslib.wms import WebMapServiceでエラーになります。

! pip install owslib

次に、WMS（URLは"https://kartta.hel.fi/ws/geoserver/avoindata/wms" ）からデータを取得します。
前半は、複数レイヤーを取得（名前を表示しているのは最初の３レイヤーだけ）しています。
後半は、avoindata:Kantakarttaレイヤーのみ、取得する範囲の座標を絞って取得します。取得が成功した場合は、プロットします。

import requests
from owslib.wms import WebMapService
from io import BytesIO
import rasterio
from rasterio.plot import show

# WMS URL from Helsingin kantakartta dataset
wms_url = "https://kartta.hel.fi/ws/geoserver/avoindata/wms"

# Connect to the WMS service
wms = WebMapService(wms_url, version="1.3.0")

# List available layers (Only print the first 3 to avoide clutter)
for i, (layer_name, layer) in enumerate(wms.contents.items()):
    if i < 3:
        print(f"Layer: {layer_name} - {layer.title}")
    else:
        print("More layers available...")
        break

# Define the layer and bounding box (for Helsinki region)
layer_name = "avoindata:Kantakartta"
bbox = (24.93, 60.16, 24.97, 60.18)  # Example bounding box in EPSG:4326 (longitude, latitude)

# Define the GetMap request parameters
wms_request_url = wms.getmap(
    layers=[layer_name],
    srs='EPSG:4326',  # Coordinate system (WGS84)
    bbox=bbox,        # Bounding box for Helsinki area
    size=(800, 800),  # Image size in pixels
    format='image/png',  # Request PNG image format
    transparent=True
)

# Extract the actual URL from the ResponseWrapper
url = wms_request_url.geturl()  # This extracts the URL string

# Fetch the map from the WMS as a PNG image
response = requests.get(url)

# Check if the request was successful
if response.status_code == 200:
    # Load the image into rasterio
    with rasterio.open(BytesIO(response.content)) as src:
        # Plot the map image
        show(src)
else:
    print(f"Failed to retrieve WMS data. Status code: {response.status_code}")

Loading a Multiband Raster File

In GIS and remote sensing, multiband raster files are commonly used to store imagery, such as satellite images, where each band represents a different wavelength of the electromagnetic spectrum (e.g., red, green, blue, near-infrared). In such files, each pixel has multiple values corresponding to the reflectance in different bands.

GISやリモートセンシングでは、各バンド（帯域）が異なる電磁スペクトルの波長（例：赤、緑、青、近赤外線）を表す衛星画像のような画像を保存するために、マルチバンドのラスタファイルがよく使われます。
そのようなファイルでは、各々のピクセルは、異なるバンドの反射率に対応する、複数の値を持っています。

For this section, we will explore how to load a multiband raster file, examine the individual bands, and visualize the data. We’ll be working with a multiband Sentinel-2 satellite image from Nuuksio National Park, located near Helsinki (sentinel2_Nuuksio_small.tif).

このセクションでは、マルチバンドのラスタファイルの読み方、個々のバンドの検証の仕方、データの可視化の仕方を探求します。

ヘルシンキの近くのヌークシオ国立公園（Nuuksio National Park）からのマルチバンドのセンチネル2衛星画像（sentinel2_Nuuksio_small.tif）を使って作業します。

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

まずは、rioxarrayをインポートします。

!pip install rioxarray

次に、ラスタファイルを読み込みます。

import rioxarray
import matplotlib.pyplot as plt
import numpy as np

Nuuksio_sentinel_path = "data/sentinel2_Nuuksio_small.tif"

# Load the multiband raster file
Nuuksio_sentinel = rioxarray.open_rasterio(Nuuksio_sentinel_path)

Let’s plot our raster:

では、ラスタデータをプロットしてみましょう。

Nuuksio_sentinel.plot()

WHat just happened? The output is probabely not what you expected!

When you try to plot a multiband raster directly using the plot() function, xarray interprets the data as a multi-dimensional array and creates a summary, such as a bar chart, to represent the entire dataset. This happens because xarray doesn’t know you’re working with geospatial raster data and tries to plot all bands together, leading to an unexpected result like a bar chart.

To correctly visualize the raster data, you need to plot individual bands or create an RGB composite. But before getting there, let’s have a look at the metadata:

何が起こったのでしょうか？結果は期待したものとはおそらく違います！

マルチバンドのラスタデータをplot()でプロットしようとすると、xarray は、データセット全体を表現するために、データを多次元の配列として解釈し、例えば、棒グラフのようなサマリーを作ります。
これは、xarrayが地理空間のラスタデータを処理しているということを理解できず、すべてのバンドをまとめてプロットしようとして、棒グラフのような意図しない結果になっています。

このラスタデータを正しく可視化するために、個々のバンドをプロットするか、RGB合成画像を作る必要があります。
しかし、その前に、このラスタデータのメタデータを見てみましょう。

# Print metadata and band information
print(Nuuksio_sentinel)

From this metadata, we learn that the raster file contains 12 bands, each with dimensions of 723x910 pixels, and the pixel values are stored as 32-bit floating point numbers. The coordinates provide the geospatial reference with x and y representing the spatial extent, while each band corresponds to a different spectral band. Additionally, the file has attributes such as a scale factor of 1.0 and no offset applied.

このメタデータから、ラスタファイルが12個のバンドを含んでいること、各バンドが723x910ピクセルの二次元であること、ピクセルは32ビットの浮動小数点の値を保持していることがわかります。
座標は、空間的な範囲を示すxとyで、地理空間の参照値を提供しています。 一方、各々のバンドは異なるスペクトルバンドに対応しています。
加えて、このファイルは、尺度係数（scale factor）：1.0、オフセット；なし、のような属性を持っています。

However, this brief metadata does not tell us which spectral band does each raster band refer to. However, this is the kind of information that you should get from the data source. In our case, Sentinel-2 User Guide is a good resource to find this information.

For Sentinel-2 imagery, the bands are predefined and consistent across products. Here is a quick reference for the Sentinel-2 bands:

しかしながら、この簡素なメタデータからは、各々のラスタバンドがどのスペクトルバンドを参照しているのかはわかりません。ただ、これがデータソースから得られる情報です。
今回の場合は、Sentinel-2 User Guide　Sentinel-2 Bandsが、この情報（どのスペクトルバンドを参照しているのか）を見つけるためのよい情報源になります。

センチネル2衛星画像では、そのバンドは前もって定義されていて、一貫性が保たれています。以下はセンチネル2のバンドの早見表です。

These are the bands in our raster data, in order from 1 to 12. Depending on our data processing needs, we can use this reference to select the relevant bands for tasks such as true color composites (using Bands 4, 3, and 2) or vegetation analysis (using Bands 8 and 4 for NDVI).

これらのバンドが今回のラスタデータの中にあり、1から12の順番になっています。
データを処理する際、必要に応じて、関連するバンドを選択するため、このリファレンスを使うことができます。

バンドの選択としては、例えば、以下があります。

• バンド4・3・2（赤色・緑色・青色）によるトゥルーカラー合成（True Color Composites）
• バンド8（近赤外線）・4（赤色）よる植生分析（Vegetation analysis）（正規化植生指数（NDVI：Normalized Difference Vegetation Index）用）

Let’s start by accessing the red band from our raster file. In (rioxarray), an individual band can be accessed using the sel() method, where you specify the band number. According to the reference above, band=4 retrieves the red band:

では、ラスタデータの赤色のバンドにアクセスすることから始めましょう。rioxarrayでは、個々のバンドはsel()メソッドでアクセスでき、その際、バンドの番号を指定します。
前述のリファレンスに基づき、band=4で赤色のバンドを取得します。

# Access red band
red_band = Nuuksio_sentinel.sel(band=4)

Let’s see what happen now if we try to plot this specific band of the raster file:

ラスタファイルの特定のバンド（band=4）をプロットして、どうなったか見てみましょう。

red_band.plot()

Now our raster file is visualized (the red band). We could also select the band and plot on the fly. let’s try for the blue band (band=2).

ラスタファイル（の赤色のバンド）が可視化されています。
すぐに別のバンドを選択して、プロットできます。では、青色のバンド（band=2）で試してみましょう。

Nuuksio_sentinel.sel(band=2).plot()

Creating an RGB Composite

In remote sensing and GIS, creating an RGB composite is a common method for visualizing multiband satellite imagery in natural or false colors. By combining the red, green, and blue bands of a satellite image, we can generate a true-color image that resembles what the human eye would see, or we can create false-color composites for enhanced analysis (e.g., vegetation health, water bodies, etc.).

リモートセンシングやGISでは、RGB合成は、マルチバンドの衛星画像を自然色（natural colors）または疑似色（false colors）で可視化するための、一般的な手段です。

衛星画像の赤色・緑色・青色のバンドを結合することで、人の目が見ているものと似たトゥルーカラーの画像を生成することができたり、拡張的な解析（例えば、植生の健康状態や水域（Water bodies）等）のために、疑似色を作成することができます。

For Sentinel-2 imagery, the red, green, and blue bands correspond to:

• Red: Band 4
• Green: Band 3
• Blue: Band 2

センチネル2衛星画像での赤色・緑色・青色のバンドはこちら：

• 赤色： Band 4
• 緑色： Band 3
• 青色： Band 2

We accomplish this by loading the individual red, green, and blue bands, stacking them together into an RGB image, and normalizing the values to be suitable for display. The Numpy’s np.dstack() function stacks arrays along a new depth axis (3rd dimension), which is useful for combining multiple 2D arrays (e.g., raster bands) into a 3D RGB image (NumPy dstack documentation).

私たちは、個々の赤色・緑色・青色のバンドを読み込み、それらをまとめてRGB画像にまとめて積み込み、表示できるように値を正規化することで、上記の処理を実現します。

Numpyのnp.dstack()関数は、配列を奥行き方向（depth方向、第3軸）に積み重ねます。
（つまり、関数の先頭の"d"はdepthの"d"）
なので、この関数は、複数の二次元配列（例えば、ラスタデータのバンド）を三次元のRGB画像として結合する際に使いやすいです。

（わかりやすくするために、np.dstack関数の引数の書き方を少し変えています）

# Extract the red, green, and blue bands
red_band = Nuuksio_sentinel.sel(band=4)
green_band = Nuuksio_sentinel.sel(band=3)
blue_band = Nuuksio_sentinel.sel(band=2)

# Stack the bands together to create an RGB composite
rgb_image = np.dstack(
    (
        red_band.values,
        green_band.values,
        blue_band.values
    )
)

# Normalize the image values between 0 and 1 by dividing by the max value
rgb_image = rgb_image / np.max(rgb_image)

# Plot the RGB composite image
plt.imshow(rgb_image)
plt.title("RGB Composite of Nuuksio National Park (Sentinel-2)")
plt.axis('off')  # Hide the axis for better visualization
plt.show()

Great! Now we can visualize the image in true color, just as it would appear to the human eye.

すばらしい！人が目で見るような、トゥルーカラーの画像が可視化できました。