geopandas

GeoPandas - 地理空間データ分析

GeoPandas を使うと、使い慣れた Pandas API で空間結合、幾何学的操作、座標変換を行えます。「geometry」を DataFrame のほかの列と同じように扱いながら、面積、距離、交差を計算できます。

使用する場面

• 空間形式（Shapefile、GeoJSON、GeoPackage、KML）の読み込みと書き込み。
• 空間結合の実行（例：「このポリゴン内にあるポイントはどれか？」）。
• 座標参照系の変換（緯度・経度からメートル単位への再投影）。
• 幾何学的分析（バッファー、重心、凸包の計算）。
• 主題図の作成（コロプレス図）。
• 空間的な関係（contains、overlaps、touches、within）の計算。
• OpenStreetMap データや衛星画像由来のベクターデータの操作。

参考ドキュメント

公式ドキュメント：https://geopandas.org/
インタラクティブチュートリアル：https://geopandas.org/en/stable/gallery/index.html
検索パターン：gpd.read_file、gdf.to_crs、gpd.sjoin、gdf.buffer、gdf.explore

基本原則

GeoDataFrame

GeoDataFrame は、少なくとも 1 つの GeoSeries 列（通常は geometry という名前）を持つ pandas.DataFrame です。各行はフィーチャー（ポイント、ライン、ポリゴン）を表します。

座標参照系（CRS）

CRS のないデータは、グリッド上の単なる数値です。実世界での計算（km² 単位の面積など）を行うには、CRS（WGS84 - EPSG:4326 や UTM など）を定義する必要があります。

述語と集合演算

空間分析では、二項述語（intersects、within、contains）と集合論的演算（union、intersection、difference）を使用します。

クイックリファレンス

インストール

pip install geopandas pyarrow pyproj fiona shapely

標準的なインポート

import geopandas as gpd
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from shapely.geometry import Point, LineString, Polygon

基本パターン - 読み込みと描画

import geopandas as gpd

# Load built-in dataset
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))

# Filter and Project
europe = world[world.continent == 'Europe']
europe = europe.to_crs(epsg=3035) # Equal Area projection for Europe

# Plot
europe.plot(column='pop_est', legend=True, cmap='viridis')

重要なルール

✅ 推奨事項

• 必ず CRS を確認する - 空間操作の前に gdf.crs を確認してください。
• 計測には投影座標系を使用する - 面積や距離を計算する前に、UTM のような投影座標系（メートル／フィート単位）を使用してください。
• 空間インデックスを使用する - 大規模なデータセットでは、gdf.sindex を使用するか、クエリを高速化するために sjoin が使用されていることを確認してください。
• ジオメトリを検証する - gdf.is_valid を使って、壊れたポリゴン（自己交差）を見つけてください。
• 可視化のために簡略化する - gdf.simplify() を使って、複雑な境界線の描画を高速化してください。
• .explore() を使用する - Jupyter で手軽にインタラクティブマップを作成できます（Leaflet／Folium を使用）。

❌ 禁止事項

• 度単位で面積を計測する - 緯度・経度の CRS（EPSG:4326）で .area を計算しないでください。結果は「平方度」（意味のない単位）になります。
• ループで反復処理する - 行ごとのループを避け、ベクトル化された空間操作を使用してください。
• トポロジーを無視する - 「touches」と「intersects」は異なること（境界と内部）に注意してください。
• アクティブジオメトリの設定を忘れる - GeoDataFrame に複数のジオメトリ列がある場合は、gdf.set_geometry() で使用する列を指定してください。

アンチパターン（絶対に避ける）

# ❌ BAD: Manual distance calculation on Lat/Lon (ignores Earth's curvature)
def dist(p1, p2):
    return np.sqrt((p1.x - p2.x)**2 + (p1.y - p2.y)**2)

# ✅ GOOD: Reproject and use vectorized distance
gdf = gdf.to_crs(epsg=3857) # Web Mercator (meters)
distances = gdf.distance(other_point)

# ❌ BAD: Manually checking points in polygons
for i, poly in countries.iterrows():
    for j, pt in cities.iterrows():
        if poly.geometry.contains(pt.geometry):
            print("Found")

# ✅ GOOD: Spatial Join (Optimized with spatial index)
cities_with_country = gpd.sjoin(cities, countries, predicate='within')

ジオメトリの作成と操作

座標からの作成

# From a Pandas DataFrame with Lat/Lon
df = pd.DataFrame({'City': ['NY', 'London'], 'Lat': [40.7, 51.5], 'Lon': [-74.0, -0.1]})
gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.Lon, df.Lat))
gdf.set_crs(epsg=4326, inplace=True)

幾何学的操作

# Centroids and Envelopes
gdf['centroid'] = gdf.centroid
gdf['envelope'] = gdf.envelope # Bounding box

# Buffering (Creating a zone around features)
# Warning: Do this in a projected CRS (meters)
stations_buffered = metro_stations.to_crs(epsg=32633).buffer(500) # 500 meters

# Unifying overlapping polygons
total_area = gdf.union_all()

空間クエリ

空間結合（sjoin）

# Find which district each school belongs to
schools_in_districts = gpd.sjoin(schools, districts, how="inner", predicate="within")

# sjoin types:
# 'intersects' (default), 'contains', 'within', 'touches', 'overlaps', 'crosses'

オーバーレイ（集合演算）

# Intersection of two layers (e.g., protected area vs. forest)
forest_in_park = gpd.overlay(forests, parks, how='intersection')

# 'union', 'intersection', 'difference', 'symmetric_difference'

座標参照系（CRS）

再投影

# Checking the current CRS
print(gdf.crs)

# Reprojecting to a specific EPSG code
# EPSG:4326 -> WGS84 (Degrees, used for GPS)
# EPSG:3857 -> Web Mercator (Meters, used for web maps)
gdf_meters = gdf.to_crs(epsg=3857)

# Match CRS of another layer
gdf_2 = gdf_2.to_crs(gdf_1.crs)

可視化

静的マップとインタラクティブマップ

# Layered plotting
fig, ax = plt.subplots(figsize=(10, 10))
base = countries.plot(ax=ax, color='white', edgecolor='black')
cities.plot(ax=base, marker='o', color='red', markersize=5)

# Interactive exploration (requires folium)
cities.explore(column='population', cmap='magma', m=None)

実践的なワークフロー

1. 近接分析（ポイントがバッファー内にあるか）

def find_entities_near_road(roads, entities, distance_m=1000):
    """Find all entities within 1km of any road."""
    # 1. Project to a metric CRS (e.g., UTM)
    roads_m = roads.to_crs(epsg=3857)
    entities_m = entities.to_crs(epsg=3857)
    
    # 2. Create buffer
    road_buffer = roads_m.buffer(distance_m)
    
    # 3. Create a GeoDataFrame from buffer to use in sjoin
    buffer_gdf = gpd.GeoDataFrame(geometry=road_buffer, crs=roads_m.crs)
    
    # 4. Spatial Join
    nearby = gpd.sjoin(entities_m, buffer_gdf, predicate='within')
    return nearby

2. 境界に合わせたデータのクリッピング

def clip_data(data, boundary):
    """Clip a large vector dataset to a specific boundary polygon."""
    return gpd.clip(data, boundary)

# Usage: city_parks = clip_data(national_parks, city_limits)

3. 面積の被覆率の計算

def calculate_land_use_pct(region, land_use_layer):
    """Calculate what % of 'region' is covered by each land use type."""
    # Ensure CRS matches and is projected
    land_use_layer = land_use_layer.to_crs(region.crs)
    
    # Intersect region with land use
    intersections = gpd.overlay(land_use_layer, region, how='intersection')
    
    # Calculate area
    intersections['area'] = intersections.area
    total_area = region.area.sum()
    
    return intersections.groupby('class')['area'].sum() / total_area * 100

パフォーマンスの最適化

空間インデックス（sindex）の使用

# Check if a point is within any polygon in a large GDF efficiently
spatial_index = countries.sindex

# Find possible matches using bounding boxes first
possible_matches_index = list(spatial_index.intersection(target_point.bounds))
possible_matches = countries.iloc[possible_matches_index]

# Precise check only on candidates
precise_match = possible_matches[possible_matches.intersects(target_point)]

大容量ファイルの読み込み（Parquet）

# GeoJSON is slow to read/write. GeoParquet is significantly faster and smaller.
gdf.to_parquet("large_data.parquet")
gdf = gpd.read_parquet("large_data.parquet")

よくある落とし穴と解決策

CRS の不一致

# ❌ Problem: sjoin returns 0 results even if data looks overlapping
# ✅ Solution: Align CRS
if cities.crs != districts.crs:
    cities = cities.to_crs(districts.crs)

無効なジオメトリ（自己交差）

# ❌ Problem: Overlay or Area calculation fails/gives weird results
# ✅ Solution: Fix with buffer(0) or check validity
invalid = gdf[~gdf.is_valid]
gdf['geometry'] = gdf['geometry'].buffer(0) # Common trick to fix minor topology errors

バッファーによるメモリ不足

# ❌ Problem: Buffering millions of points with high resolution
# ✅ Solution: Use low resolution or simplify
gdf.buffer(100, resolution=4) # Default is 16, 4 is often enough for analysis

GeoPandas は、従来の GIS ソフトウェアと Python のデータサイエンススタックをつなぎます。Pandas のコードを 1 行書くような手軽さで、空間分析を行えます。