GeoPandasの公式チュートリアル実践

GeoPandasとは

GeoPandasは、地理空間データ（点・線・面などのジオメトリを持つデータ）をPythonで扱うためのライブラリで、pandasのDataFrameに地理機能を拡張したものです。表形式データの操作性（フィルタ、結合、集計など）と、Shapelyによるジオメトリ演算、PyProjによる座標参照系（CRS）変換、Fiona/GeoJSON I/O、Matplotlibによるプロットを一体化して扱えます。

主な特徴

  • pandasと同じ感覚で、ジオメトリ列（GeoSeries）を含むGeoDataFrameを操作できる
  • 交差・包含・バッファ・ユニオンなどの幾何演算を簡潔に記述可能
  • 座標参照系（CRS）の保持と変換が容易（to_crs）
  • Shapefile、GeoJSON、GeoPackageなど主要フォーマットの読み書きに対応
  • すばやい可視化（.plot()）が可能で、背景地図（contextily等）とも連携しやすい

よく使われる用途

  • 行政界や道路網などの境界データの読み込みと可視化
  • ポイントデータ（店舗、センサー、事故地点など）の空間結合・近接解析
  • 集計結果のコロプレス図作成
  • 異なる投影法への座標変換と距離・面積計算

基本的な使い方

import geopandas as gpd

# データの読み込み（ShapefileやGeoJSONなどを自動判別）
gdf = gpd.read_file("areas.geojson")

# CRSの確認と変換（例: WGS84 -> Webメルカトル）
print(gdf.crs)                 # 例: EPSG:4326
gdf = gdf.to_crs(epsg=3857)

# フィルタや列操作（pandasと同様）
subset = gdf[gdf["population"] > 50000][["name", "population", "geometry"]]

# 空間結合（ポイントをポリゴンに割り当てるなど）
points = gpd.read_file("shops.geojson").to_crs(gdf.crs)
joined = gpd.sjoin(points, gdf, how="left", predicate="within")

# 幾何演算（バッファ、交差）
buffered = points.buffer(500)  # 500mの緩衝帯（投影座標系で）
intersections = gdf.intersection(buffered.unary_union)

# 可視化
ax = gdf.plot(edgecolor="black")
points.plot(ax=ax, markersize=5)

周辺ツールとの関係

  • Shapely: ジオメトリ型と幾何演算のエンジン
  • PyProj: 座標参照系と変換
  • Fiona: GISファイルのI/O
  • Rasterio/xarray: ラスター（画像・標高）系データとの連携に有用
  • contextily: タイル背景地図を追加して地図表現を強化

押さえておくポイント

  • 距離や面積の計算は、投影座標系（メートルなど）へto_crsしてから行う
  • 大規模データは処理が重くなるため、Dask-GeoPandasやPostGISの活用を検討
  • predicate引数（"intersects" "within" "contains" など）で空間結合の条件を明示する

GeoPandasは、pandasの操作性をそのまま地理空間データに拡張することで、調査・可視化・分析を一気通貫で行えるのが最大の利点です

公式チュートリアル

公式：https://geopandas.org/en/stable/getting_started.html

インストールするもの

pip install geopandas

チュートリアルでは「geodatasets」というパッケージのデータを使う

公式：https://geodatasets.readthedocs.io/en/latest/

pip install geodatasets

あとマップを描画するためにmatplotlibが必要

pip install matplotlib

ファイルの読み込み：gdf.read_file()

nybb = ニューヨーク行政地区のデータ

import geopandas
from geodatasets import get_path

path_to_data = get_path("nybb")
gdf = geopandas.read_file(path_to_data)

読み込み結果

ファイルの書き出し：gdf.to_file()

拡張子はファイル名から推測する（Shapefileなら.shp、GeoJSONなら.geojson）が、driverで直接指定もできる

gdf.to_file("my_file.geojson", driver="GeoJSON")

もっと詳しい情報：Reading and writing files

面積を計算する

各ポリゴン（この場合はMultiPolygon）の面積を計測するには、GeoDataFrame.area属性にアクセスします。この属性はpandas.Seriesを返します。GeoDataFrame.areaは、アクティブなジオメトリ列にGeoSeries.areaを適用したものに過ぎないことに注意してください。

先のデータセットには「dtype: geometry」の列が1つしかないので、この列が自動的に「アクティブ」となっている。よって空間メソッドを実行した際もこの列に適用される

gdf = gdf.set_index("BoroName") # 見やすいように行政区名をindexにする
gdf["area"] = gdf.area

ポリゴンの境界と重心の取得

boundary=境界

centroid＝重心

各ポリゴン（LineString）の境界を取得するには、GeoDataFrame.boundaryにアクセスします。

gdf["boundary"] = gdf.boundary

これで「dtype: geometry」の列が2つになった

さらに重心の列を作成。これで「dtype: geometry」の列が3つになった

gdf["centroid"] = gdf.centroid

boundaryとcentroidはどんな意味がある？

「boundary」と「centroid」は、ポリゴンの幾何学的な性質を表す派生ジオメトリです。

boundary（境界）

  • 意味: ポリゴンの外周（および穴があればその内周）だけを取り出した「線」の集合
  • 形: Polygon → LineString（外周）と複数あれば MultiLineString（外周＋内周）
  • 性質: 面積は持たず、座標は元のCRSのまま。ライン解析（長さ、交差判定など）に使う
  • 注意: LineString/Polygon以外では定義が異なる（例: LineStringのboundaryは両端点、Pointのboundaryは空）

centroid（重心）

  • 意味: ポリゴンの面積に基づく幾何学的重心（“ど真ん中”の代表点）
  • 形: 1点のPoint。MultiPolygonなら全体の面積加重で1点
  • 性質: 代表位置の取得、ラベリング、空間結合のキーなどに使う
  • 注意1: くびれの強い形やドーナツ状では、centroidがポリゴンの外に出ることがある（地図表示用には point_on_surface/representative_point を検討）
  • 注意2: 緯度経度（EPSG:4326）でのcentroidは“球面上の見かけの中心”で、測地学的な意味は弱い。距離や面積、重心を厳密に扱うときは投影座標系（メートル系）にto_crsしてから計算する

要するに、boundaryは「外周の線」を取り出す操作、centroidは「面の中心点」を求める操作で、どちらも元のポリゴンから導かれる派生ジオメトリです。

距離の計算

各重心が最初の重心の位置からどれだけ離れているかを測定する

first_point = gdf["centroid"].iloc[0]
gdf["distance"] = gdf["centroid"].distance(first_point)

マップを描画する

GeoPandasは地図をプロットすることもできるので、空間上でジオメトリがどのように表示されるかを確認できます。アクティブなジオメトリをプロットするには、GeoDataFrame.plot()を呼び出します。別の列で色分けするには、その列を最初の引数として渡します。以下の例では、アクティブなジオメトリ列をプロットし、「area」列で色分けしています。また、凡例も表示します（legend=True）

gdf.plot("area", legend=True)

「geometry」で形を描写し、area（面積）が大きいほど黄色になるという設定

重心をアクティブにして、同じようにマップを描画した場合

gdf = gdf.set_geometry("centroid")
gdf.plot("area", legend=True)

さらに両方を重ねて表示することもできる

ax = gdf["geometry"].plot()
gdf["centroid"].plot(ax=ax, color="black")

もっと詳しいガイド：Mapping and plotting tools

凸包の作成

geometry列をアクティブに戻す

gdf = gdf.set_geometry("geometry")

凸包列を作成

gdf["convex_hull"] = gdf.convex_hull

凸包を半透明にして描画

ax = gdf["convex_hull"].plot(alpha=0.5)
gdf["boundary"].plot(ax=ax, color="white", linewidth=0.5)

バッファの作成

バッファ列の作成

gdf["buffered"] = gdf.buffer(10000)
gdf["buffered_centroid"] = gdf["centroid"].buffer(10000)

マップの描画

ax = gdf["buffered"].plot(alpha=0.5)
gdf["buffered_centroid"].plot(ax=ax, color="red", alpha=0.5)
gdf["boundary"].plot(ax=ax, color="white", linewidth=0.5)

凸包とバッファはどんな意味を持つ？

Convex hull（凸包）とBuffer（バッファ）は、元のジオメトリから新しいジオメトリを「構成的」に作る代表的な操作です。GeoPandasではそれぞれ gdf.convex_hull、gdf.buffer(distance, ...) で得られます。

Convex hull（凸包）

  • 意味: 対象の点群や図形をすべて包み込む「最小の凸多角形」。2点なら線分、1点なら点に退化します。
  • 使いどころ: 点群の全体形状の概略化、外接ポリゴンの作成、範囲推定や衝突判定の前処理など。
  • 出力: Polygon（場合によりLineString/Pointに退化）。
  • 備考: GeoSeries.convex_hull は各ジオメトリの凸包を返します。より“くびれ”のある形がほしい場合は concave hull（凹包）もあります

Buffer（バッファ）

  • 意味: ジオメトリの周囲「指定距離以内のすべての点」の集合。点や線の周りに“帯（面）”を作ります。負の距離で内側に縮めることも可能です。
  • 使いどころ: 近接圏の作成（例: 駅から500m圏）、線路の保護帯、ポリゴンの膨張・収縮、距離ベースの選別や空間結合の前処理。
  • 出力: ふつうはPolygon（線は帯状、点は円状）。cap_style（端の形）、join_style（角の形）、single_sided などで形状を制御できます。
  • 注意: 距離の単位はCRSに依存します。緯度経度（EPSG:4326）のままでは度単位になるため、メートル系の投影（UTMや等距離・等積系）へ to_crs してから距離をメートルで指定するのが実務的です

ブルックリン周辺の確認

まずはブルックリンのポリゴンを取得する

brooklyn = gdf.loc["Brooklyn", "geometry"]

そしてバッファリングしたデータとどこが交差するか確認

gdf["buffered"].intersects(brooklyn)

バッファは1万フィートなので、Bronxだけ1万フィート以上離れてることがわかる

あるいは、バッファリングされた重心が元の行政区ポリゴン内に完全に含まれているかどうかを確認することもできます。この場合、両方のGeoSeriesが整列され、各行ごとにチェックが実行されます。

gdf["within"] = gdf["buffered_centroid"].within(gdf)
gdf["within"]

そしてこれを描画して確認

gdf = gdf.set_geometry("buffered_centroid")
ax = gdf.plot(
    "within", legend=True, categorical=True, legend_kwds={"loc": "upper left"}
)
gdf["boundary"].plot(ax=ax, color="black", linewidth=0.5)

Projections（地図投影・投影法）とは

Projections（地図投影・投影法）は、地球のような曲面（楕円体／球）の位置を、平面の座標（x, y）へ写像する数学的な方法です。地図やGISで平面上に描くために不可欠で、どの投影を選ぶかで「何を正しく保つか（面積・形・距離・方位）」が変わります。

基本ポイント

  • 目的: 曲面上の経緯度を平面座標へ変換して、距離・面積・長さを扱いやすくする
  • CRSとの関係: CRS（座標参照系）は「地理座標系（経度・緯度）」または「投影座標系（平面のx, y）」の定義。投影法は投影座標系の一部要素
  • 歪み: 地球を平面に写す以上、面積・形状・距離・方位のすべてを同時に厳密には保てない

代表的な投影の性質

  • 等積（Equal-area）: 面積を正しく保つ（Albers, Mollweide など）。統計地図の面積比較に適する
  • 正角（Conformal）: 小領域の形と角度を保つ（Mercator, Transverse Mercator, Lambert Conformal Conic）。角度・形重視、航海図や大縮尺図に適する
  • 等距離（Equidistant）: ある基準からの距離を正しく保つ（Azimuthal Equidistant など）。距離測定の特定用途
  • 妥協（Compromise）: 歪みを平均的に抑える（Robinson, Winkel Tripel など）。世界地図の見栄え重視

実務での選び方（例）

  • 面積や統計のコロプレス図: 等積投影（地域別の面積比較が必要なとき）
  • 角度・形状重視、施設図・道路網解析: 正角投影（UTMやローカルのTM）
  • 距離や到達圏（バッファ）重視: 等距離系またはローカルのメートル系投影
  • 広域かつ見た目重視: 妥協投影

日本でよく使う投影座標系の例

  • UTM（ユニバーサル横メルカトル）ゾーン52–54（地域により使い分け）
  • 平面直角座標系（JGD2011 / EPSG:6668–6687 など、計19系）
  • いずれもメートル単位で扱えるため、距離・面積計算に適する

GeoPandasでの扱い

  • 現在のCRSの確認: gdf.crs（例: EPSG:4326 は経緯度、単位は度）
  • 投影の変更（再投影）: gdf.to_crs(epsg=<投影EPSG>)
  • 面積・距離・バッファは投影座標系で計算するのが基本（単位は投影の座標単位に従う）

注意点

  • Webメルカトル（EPSG:3857）は面積保存ではないため、正確な面積用途には不向き
  • 「距離」「面積」「バッファ」は意図した単位の投影で実施する
  • 投影の選択は対象地域の広がり・緯度・用途（面積重視か距離重視か）で決める

ニューヨーク行政区の投影法

各GeoSeriesには、GeoSeries.crsでアクセスできる座標参照系（CRS）があります。CRSは、GeoPandasにジオメトリの座標が地球上のどこに位置しているかを伝えます。

場合によっては、CRSは地理座標系（緯度と経度）であり、その場合のCRSはWGS84で、出典コードはEPSG:4326です。では、ニューヨーク行政区のGeoDataFrameの投影法を見てみましょう。

gdf.crs

ジオメトリはEPSG:2263形式で表現され、座標はフィート単位です。GeoSeries.to_crs() を使えば、GeoSeriesをEPSG:4326などの別のCRSに簡単に再投影できます。

gdf = gdf.set_geometry("geometry")
boroughs_4326 = gdf.to_crs("EPSG:4326")
boroughs_4326.plot()
boroughs_4326.crs

プロットの軸に沿った座標の違いに注目してください。以前は120,000～280,000（フィート）でしたが、現在は40.5～40.9（度）になっています。この場合、boroughs_4326にはWGS84の「geometry」列がありますが、その他の情報（重心など）は元のCRSのままです。

距離や面積に依存する演算では、地理座標系（度）ではなく、投影座標系（メートル、フィート、キロメートルなど）を使用する必要があります。GeoPandasの演算は平面座標系ですが、度は球面上の位置を反映します。そのため、度を使用した空間演算では正しい結果が得られない場合があります。例えば、gdf.area.sum()（投影座標系）の結果は8,429,911,572平方フィートですが、boroughs_4326.area.sum()（地理座標系）の結果は0.083です。

重要なのは「小さいからダメ」というより「単位が度²・度で、平面の距離や面積の単位（m、ft、km²など）ではない」ため、そもそも比較不能・不正確という点です。

ポイントの整理

  • 地理座標系（EPSG:4326など）での area は「度²」。度は角度なので、面積の物理単位ではない
  • 緯度によって“1度の長さ”が変わるため、度・度²は場所ごとに実長が異なる
  • GeoPandas／Shapelyの幾何演算は平面（planar）前提なので、距離・面積を正しく扱うには投影座標系（メートル・フィート系）に変換してから計算する必要がある
  • 引用の例で 8,429,911,572 ft² と 0.083 の差が出ているのは、前者が投影CRSでの面積（平方フィート）、後者が地理CRSでの度²だから

実務での対処

  • 面積・距離・バッファを行う前に、適切な投影へ再投影する
    等積投影（面積用途）や、対象地域に合ったUTM・ローカル投影（距離用途）がおすすめ
# 例：メートルで扱いたい場合（地域に合うEPSGを選ぶ）
gdf_m = gdf.to_crs(<meters-based EPSG>)
gdf_m["area_m2"] = gdf_m.area
  • フィートが必要ならフィート系の投影CRSへ
gdf_ft = gdf.to_crs(<feet-based EPSG>)
gdf_ft["area_ft2"] = gdf_ft.area

補足

  • 正確な“地球楕円体”上の面積が必要なら、pyproj.Geod の測地線面積を使う手もある
  • バッファや距離計測でも同じで、度のまま距離を入れると「度」単位のバッファになってしまうため、必ず投影してから行うこと

-GeoPandas, Python