alterlab-geomaster

GeoMaster

GIS、リモートセンシング、空間解析、地球観測向け機械学習を網羅する地理空間科学スキルです。70以上のトピックと、8つのプログラミング言語による500以上のコード例を収録しています。

このスキルを使用する場面

次のような要望がある場合に、このスキルを使用します。

• 衛星画像や航空画像（Sentinel、Landsat、MODIS、SAR、ハイパースペクトル）を処理し、スペクトル指数を計算する
• ラスターやDEMを扱う — 地形指標、モザイク、再投影、ゾーン統計
• クラウドネイティブデータ（STAC、COG、Planetary Computer、Earth Engine、Copernicus Data Space）を検索して読み込む
• 地球観測データに機械学習を学習または適用する、あるいは水文学・海洋・大気の空間ワークフローを実行する

該当しない場合

シナリオ	代わりに使用
ラスター／地球観測要素を含まない純粋なベクター処理 — sjoin、バッファー、オーバーレイ、ディゾルブ、コロプレス図	alterlab-geopandas
天体座標、FITS画像、天空WCS	alterlab-astropy
流体力学そのもののシミュレーション（Navier-Stokes、浅水方程式ソルバー）	alterlab-fluidsim
人文科学研究向けに史料をマッピングする、または歴史地図にジオリファレンスを付与する	alterlab-digital-humanities

インストール

GDALを基盤とするスタックでは、パッケージマネージャーを1つだけ選択してください。rasterio／fiona／pyproj／shapelyの最新のpipホイールにはGDAL／GEOS／PROJが同梱されているため、pipのみの（uv）環境でも、システムにGDALをインストールせずにコアライブラリを利用できます。同じ環境でconda-GDALとpip版rasterio／fionaを混在させないでください。両者には異なるGDALバイナリが含まれており、ABIの不一致によってセグメンテーション違反が発生します。単独のgdal Pythonバインディングにはバイナリが同梱されていないため、一致するシステム／condaのlibgdalが必要です。PDALとrsgislibにも信頼できるpipホイールがないため、condaでインストールしてください（オプションB）。

# Option A — uv/pip env (recommended here): wheels bundle GDAL/GEOS/PROJ.
# rasterio/fiona cover most GDAL needs; standalone `gdal` -> use Option B.
uv pip install rasterio fiona shapely pyproj geopandas
uv pip install torchgeo earthengine-api
uv pip install scikit-learn xgboost torch-geometric
uv pip install osmnx networkx folium keplergl
uv pip install cartopy contextily mapclassify
uv pip install xarray rioxarray dask-geopandas
uv pip install pystac-client planetary-computer odc-stac rio-cogeo
uv pip install "laspy[lazrs]" open3d      # PDAL: use conda (no pip wheel)

# Option B — conda env (best for PDAL / rsgislib / system GDAL tooling)
conda install -c conda-forge gdal rasterio fiona shapely pyproj geopandas \
    rsgislib pdal python-pdal postgis libspatialite

クイックスタート

Sentinel-2からのNDVI

import rasterio
import numpy as np

with rasterio.open('sentinel2.tif') as src:
    red = src.read(4).astype(float)   # B04
    nir = src.read(8).astype(float)   # B08
    ndvi = (nir - red) / (nir + red + 1e-8)
    ndvi = np.nan_to_num(ndvi, nan=0)

    profile = src.profile
    profile.update(count=1, dtype=rasterio.float32)

    with rasterio.open('ndvi.tif', 'w', **profile) as dst:
        dst.write(ndvi.astype(rasterio.float32), 1)

GeoPandasによる空間解析

import geopandas as gpd

# Load and ensure same CRS
zones = gpd.read_file('zones.geojson')
points = gpd.read_file('points.geojson')

if zones.crs != points.crs:
    points = points.to_crs(zones.crs)

# Spatial join and statistics
joined = gpd.sjoin(points, zones, how='inner', predicate='within')
stats = joined.groupby('zone_id').agg({
    'value': ['count', 'mean', 'std', 'min', 'max']
}).round(2)

Google Earth Engine時系列

import ee
import pandas as pd

ee.Initialize(project='your-project')
roi = ee.Geometry.Point([-122.4, 37.7]).buffer(10000)

s2 = (ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
      .filterBounds(roi)
      .filterDate('2020-01-01', '2023-12-31')
      .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20)))

def add_ndvi(img):
    return img.addBands(img.normalizedDifference(['B8', 'B4']).rename('NDVI'))

s2_ndvi = s2.map(add_ndvi)

def extract_series(image):
    stats = image.reduceRegion(ee.Reducer.mean(), roi.centroid(), scale=10, maxPixels=1e9)
    return ee.Feature(None, {'date': image.date().format('YYYY-MM-dd'), 'ndvi': stats.get('NDVI')})

series = s2_ndvi.map(extract_series).getInfo()
df = pd.DataFrame([f['properties'] for f in series['features']])
df['date'] = pd.to_datetime(df['date'])

コア概念

データ型

型	例	ライブラリ
ベクター	Shapefile、GeoJSON、GeoPackage	GeoPandas、Fiona、GDAL
ラスター	GeoTIFF、NetCDF、COG	Rasterio、Xarray、GDAL
点群	LAS、LAZ	Laspy、PDAL、Open3D

座標系

• EPSG:4326（WGS 84）- 地理座標、緯度／経度。保存に使用します
• EPSG:3857（Web Mercator）- ウェブマップ専用（面積／距離の計算には使用しないでください！）
• EPSG:326xx/327xx（UTM）- メートル単位の計算に使用します。ゾーンあたりの歪みは1%未満です
• UTMを自動検出するには、gdf.estimate_utm_crs()を使用します

# Always check CRS before operations
assert gdf1.crs == gdf2.crs, "CRS mismatch!"

# For area/distance calculations, use projected CRS
gdf_metric = gdf.to_crs(gdf.estimate_utm_crs())
area_sqm = gdf_metric.geometry.area

OGC標準

• WMS：Web Map Service — ラスターマップ
• WFS：Web Feature Service — ベクターデータ
• WCS：Web Coverage Service — ラスターカバレッジ
• STAC：Spatiotemporal Asset Catalog — 最新のメタデータ

一般的な操作

スペクトル指数

def calculate_indices(image_path, bands=(2, 3, 4, 8, 11)):
    """NDVI, EVI, SAVI, NDWI from Sentinel-2.

    `bands` maps (B02, B03, B04, B08, B11) to the 1-based band indices in
    your file. EVI's constants (6, 7.5, +1) assume surface reflectance in
    [0, 1]; raw L2A DN are scaled by 10000, so divide first or EVI is wrong.
    """
    with rasterio.open(image_path) as src:
        B02, B03, B04, B08, B11 = (src.read(b).astype(float) for b in bands)

    ndvi = (B08 - B04) / (B08 + B04 + 1e-8)
    evi = 2.5 * (B08 - B04) / (B08 + 6*B04 - 7.5*B02 + 1 + 1e-8)
    savi = ((B08 - B04) / (B08 + B04 + 0.5)) * 1.5
    ndwi = (B03 - B08) / (B03 + B08 + 1e-8)

    return {'NDVI': ndvi, 'EVI': evi, 'SAVI': savi, 'NDWI': ndwi}

ベクター操作

# Buffer (use projected CRS! 1000 = metres in UTM).
# Keep the result in the projected frame; don't bolt projected geometry
# back onto a geographic gdf or you silently mix CRS in one column.
gdf_proj = gdf.to_crs(gdf.estimate_utm_crs())
gdf_proj['buffer_1km'] = gdf_proj.geometry.buffer(1000)

# Spatial relationships
intersects = gdf[gdf.geometry.intersects(other_geometry)]
contains = gdf[gdf.geometry.contains(point_geometry)]

# Geometric operations
gdf['centroid'] = gdf.geometry.centroid
gdf['simplified'] = gdf.geometry.simplify(tolerance=0.001)

# Overlay operations
intersection = gpd.overlay(gdf1, gdf2, how='intersection')
union = gpd.overlay(gdf1, gdf2, how='union')

地形解析

def terrain_metrics(dem_path, azimuth=315.0, altitude=45.0):
    """Slope (deg), aspect (deg clockwise from north, downslope), hillshade [0, 1].

    Assumes a projected, north-up DEM whose horizontal units match its elevation
    units (reproject geographic DEMs first, or the gradients are meaningless).
    """
    with rasterio.open(dem_path) as src:
        dem = src.read(1).astype(float)
        xres, yres = src.res                      # pixel size — gradients need it

    dz_drow, dz_dx = np.gradient(dem, yres, xres)  # rows run north -> south
    dz_dy = -dz_drow                               # northward gradient
    slope = np.degrees(np.arctan(np.hypot(dz_dx, dz_dy)))
    aspect = np.degrees(np.arctan2(-dz_dx, -dz_dy)) % 360

    # Hillshade: cos(zenith)cos(slope) + sin(zenith)sin(slope)cos(azimuth - aspect)
    zenith, az = np.radians(90.0 - altitude), np.radians(azimuth)
    s, a = np.radians(slope), np.radians(aspect)
    hillshade = np.clip(np.cos(zenith) * np.cos(s) +
                        np.sin(zenith) * np.sin(s) * np.cos(az - a), 0, 1)

    return slope, aspect, hillshade

ネットワーク解析

import osmnx as ox
import networkx as nx

# Download and analyze street network
G = ox.graph_from_place('San Francisco, CA', network_type='drive')
G = ox.add_edge_travel_times(ox.add_edge_speeds(G))  # both return the graph

# Shortest path
orig = ox.distance.nearest_nodes(G, -122.4, 37.7)
dest = ox.distance.nearest_nodes(G, -122.3, 37.8)
route = nx.shortest_path(G, orig, dest, weight='travel_time')

画像分類

from sklearn.ensemble import RandomForestClassifier
import rasterio
from rasterio.features import rasterize

def classify_imagery(raster_path, training_gdf, output_path):
    """Train RF and classify imagery."""
    with rasterio.open(raster_path) as src:
        image = src.read()
        profile = src.profile
        transform = src.transform

    # Extract training data
    X_train, y_train = [], []
    for _, row in training_gdf.iterrows():
        mask = rasterize([(row.geometry, 1)],
                        out_shape=(profile['height'], profile['width']),
                        transform=transform, fill=0, dtype=np.uint8)
        pixels = image[:, mask > 0].T
        X_train.extend(pixels)
        y_train.extend([row['class_id']] * len(pixels))

    # Train and predict
    rf = RandomForestClassifier(n_estimators=100, max_depth=20, n_jobs=-1)
    rf.fit(X_train, y_train)

    prediction = rf.predict(image.reshape(image.shape[0], -1).T)
    prediction = prediction.reshape(profile['height'], profile['width'])

    profile.update(dtype=rasterio.uint8, count=1)
    with rasterio.open(output_path, 'w', **profile) as dst:
        dst.write(prediction.astype(rasterio.uint8), 1)

    return rf

最新のクラウドネイティブワークフロー

STAC + Planetary Computer

import pystac_client
import planetary_computer
import odc.stac

# Search Sentinel-2 via STAC
catalog = pystac_client.Client.open(
    "https://planetarycomputer.microsoft.com/api/stac/v1",
    modifier=planetary_computer.sign_inplace,
)

search = catalog.search(
    collections=["sentinel-2-l2a"],
    bbox=[-122.5, 37.7, -122.3, 37.9],
    datetime="2023-01-01/2023-12-31",
    query={"eo:cloud_cover": {"lt": 20}},
)

# Load as xarray (cloud-native!)
# pystac-client >= 0.7: use search.items() (get_items() is deprecated).
items = list(search.items())[:5]
data = odc.stac.load(
    items,
    bands=["B02", "B03", "B04", "B08"],
    crs="EPSG:32610",
    resolution=10,
)

# Calculate NDVI on xarray
ndvi = (data.B08 - data.B04) / (data.B08 + data.B04)

クラウド最適化GeoTIFF（COG）

import rasterio
from rasterio.session import AWSSession
from rasterio.windows import Window

# Read a window of a cloud-hosted COG (HTTP range requests fetch only that window).
# Credentials go on a rasterio.Env, not on rasterio.open; for public buckets use
# AWSSession(aws_unsigned=True).
with rasterio.Env(AWSSession(aws_access_key_id=..., aws_secret_access_key=...)):
    with rasterio.open('s3://bucket/path.tif') as src:
        subset = src.read(1, window=Window(col_off=1000, row_off=1000,
                                           width=1000, height=1000))

# Write a valid COG (tiled + internal overviews + COG layout) with rio-cogeo;
# a plain tiled GeoTIFF without overviews is not a COG
from rio_cogeo.cogeo import cog_translate, cog_validate
from rio_cogeo.profiles import cog_profiles

cog_translate('input.tif', 'output_cog.tif', cog_profiles.get('deflate'))
is_valid, errors, warnings = cog_validate('output_cog.tif')

パフォーマンスのヒント

# 1. Spatial indexing (10-100x faster queries)
gdf.sindex  # Auto-created by GeoPandas

# 2. Chunk large rasters
with rasterio.open('large.tif') as src:
    for i, window in src.block_windows(1):
        block = src.read(1, window=window)

# 3. Dask for big data (lazy, chunked, dask-backed DataArray)
import rioxarray
da_raster = rioxarray.open_rasterio('large.tif', chunks=(1, 1024, 1024))
# .data is the underlying dask.array if you need the raw chunked array:
# dask_array = da_raster.data

# 4. Use Arrow for I/O
gdf.to_file('output.gpkg', use_arrow=True)

# 5. GDAL caching
from osgeo import gdal
gdal.SetCacheMax(2**30)  # 1GB cache

# 6. Parallel processing
rf = RandomForestClassifier(n_jobs=-1)  # All cores

ベストプラクティス

1. 空間操作の前に、必ずCRSを確認します
2. 面積／距離の計算には、投影CRSを使用します
3. ジオメトリを検証します：gdf = gdf[gdf.is_valid]（gdf.geometry.make_valid()で修復します）
4. 欠損ジオメトリを削除します：gdf = gdf[gdf.geometry.notna() & ~gdf.geometry.is_empty]（fillna(None)はジオメトリ列では機能しません）
5. 効率的な形式を使用します：GeoPackage > Shapefile、大規模データにはParquet
6. 光学画像にクラウドマスキングを適用します
7. 再現可能な研究のために、データの系譜を保持します
8. 解析スケールに適した解像度を使用します

詳細ドキュメント

• 座標系 - CRS の基礎、UTM、変換
• コアライブラリ - GDAL、Rasterio、GeoPandas、Shapely
• リモートセンシング - 衛星ミッション、スペクトル指数、SAR
• 機械学習 - RS 向けのディープラーニング、CNN、GNN
• GIS ソフトウェア - QGIS、ArcGIS、GRASS との統合
• 科学分野 - 海洋、水文学、農業、林業
• 高度な GIS - 3D GIS、時空間、トポロジー
• ビッグデータ - 分散処理、GPU アクセラレーション
• 業界での応用 - 都市計画、災害管理
• プログラミング言語 - Python、R、Julia、JS、C++、Java、Go、Rust
• データソース - 衛星カタログ、API
• トラブルシューティング - よくある問題、デバッグ、エラーリファレンス
• コード例 - 500 以上の例

---

GeoMaster は、基本的な GIS 操作から高度なリモートセンシング、機械学習まで、あらゆる内容を網羅しています。

AlterLab Academic Skills スイートの一部です。