第02回座標変換と投影変換¶

科目：ビッグデータ応用Ⅱ〔特徴量エンジニアリング〕
担当：小野原彩香
対象：3年生（後期）

この授業の進め方¶

この授業は個人作業・自己ペース型です。以下を前提として設計しています。

ハンズオン手順は授業時間内に全部終わらせる必要はありません。続きは持ち帰って自分のペースで進めてください。
演習課題の締め切りは2週間後です。
グループワークや口頭発表はありません。全作業は個人で行います。
質問はメールで受け付けています（返信は数日以内）。
コマンドは一行ずつコピー&ペーストして実行するのが確実です。
各手順の後に「うまくいかないとき」の対処法を記載しています。

学習目標¶

座標参照系（CRS）の概念を理解し、主要なEPSGコードを説明できる
ogrinfo でベクターデータのCRSと属性を確認できる
ogr2ogr で異なるCRS間の座標変換を実行できる
gdalinfo でラスターデータのCRSとメタデータを確認できる
gdalwarp でラスターデータの投影変換を実行できる
発展：緯度・経度の「1度」が場所によって何メートルになるかを、自分で導ける

1. 座標参照系（CRS）の基礎¶

1.1 座標参照系とは¶

地球上の位置を数値で表すためのルールを座標参照系（Coordinate Reference System; CRS）といいます。

「北緯43度、東経141度が札幌」というとき、その「度」の基準になっているのがCRSです。CRSが違うと、同じ数値でも指す場所が微妙にずれます。

CRSには大きく2種類があります。

種類	説明	単位	用途
地理座標系（Geographic CRS）	緯度・経度で位置を表す。地球を楕円体として扱う	度（°）	GPSデータ、広域マップ
投影座標系（Projected CRS）	地球の曲面を平面に投影して表す	メートル（m）	距離・面積の計算、地形図

大事なポイント: 距離を「メートル」で正確に計算したいときは、投影座標系を使う必要があります。緯度・経度の「度」のままでは、距離の計算が正確にできません。

導出①：緯度・経度の「1度」は、場所によって何メートルか変わる

なぜ「度」のままでは距離が計算できないのか。地球を半径 \(R\) の球とみなして考えます。

緯度方向（南北）：子午線は地球を縦にひと回りする大円です。その1周が \(2\pi R\) で、それを360度に分けるので、緯度1度ぶんの長さは

\[\frac{2\pi R}{360} = \frac{\pi R}{180}\]

となり、どの緯度でもほぼ同じです。

経度方向（東西）：一方、緯度 \(\varphi\) の緯線（同じ緯度を結んだ円）は、赤道の円より小さい円です。その半径は \(R\cos\varphi\) になります。したがって経度1度ぶんの長さは

\[\frac{2\pi R\cos\varphi}{360} = \frac{\pi R}{180}\cos\varphi\]

で、高緯度ほど \(\cos\varphi\) の分だけ縮みます。極（\(\varphi = 90°\)）では \(\cos 90° = 0\) なので、経度が何度違っても距離は0です。

数値で確かめる：WGS84の楕円体で計算すると、この資料に出てくる札幌（\(\varphi = 43.062096°\)）では次のようになります。

向き	1度あたりの長さ
緯度（南北）	111,093.9 m
経度（東西）	81,458.8 m

比は \(81{,}458.8 / 111{,}093.9 = 0.7332\) で、\(\cos(43.062096°) = 0.7306\) とほぼ一致します。わずかな差は、地球が完全な球ではなく扁平だからです。

赤道と比べると違いがはっきりします。

緯度	緯度1度	経度1度
0°（赤道）	110,574.3 m	111,319.5 m
43.062096°（札幌）	111,093.9 m	81,458.8 m
60°	111,412.3 m	55,800.0 m

札幌で「0.01度」動かすと、北へなら 1,110.9 m ですが、東へなら 814.6 m です。同じ「0.01度」で 296.4 m も違います。これが「度のままでは距離の計算が正確にできない」の中身です。

楕円体の正確な半径の式（子午線曲率半径・卯酉線曲率半径）や、平面直角座標系が実際にどういう式で平面へ写しているかは、この授業で扱う範囲外です。ここでは「経度1度は \(\cos\varphi\) で縮む」という性質までを押さえてください。

使った数学：円周の長さ \(2\pi R\) と、弧を角度で割ること。三角比 \(\cos\)。

1.2 主要なCRSとEPSGコード¶

CRSは EPSG（European Petroleum Survey Group）コードという番号で管理されています。番号を覚えれば、どのCRSかすぐに指定できます。

WGS84（EPSG:4326）¶

世界測地系1984。GPSが採用している世界標準の座標系
単位：十進法の緯度・経度（例: 43.0618° N, 141.3545° E が札幌）
Google Maps、OpenStreetMap、GeoJSONのデフォルト
距離・面積の計算には不向き（単位が「度」なので）

例: 札幌市の座標（WGS84）
  経度（Longitude）: 141.354370°
  緯度（Latitude）:   43.062096°

JGD2011（EPSG:6668）¶

日本測地系2011（Japanese Geodetic Datum 2011）
東日本大震災後（2011年）に改定された日本の標準測地系
WGS84とほぼ同じ精度だが、日本国内の精度が高い
国土地理院のデータはこの座標系を使用することが多い

平面直角座標系（EPSG:6669〜6687）¶

日本全国を19のゾーンに分け、各ゾーンで平面投影した座標系
単位：メートル。距離・面積の計算に適している
北海道の大部分は第XII系（EPSG:6680）

系番号	EPSGコード	主な対象地域
第I系	EPSG:6669	長崎県、鹿児島県の一部
第II系	EPSG:6670	福岡県、熊本県など
……	……	……
第XII系	EPSG:6680	北海道（大部分）
第XIII系	EPSG:6681	北海道（東部）

1.3 CRSの変換が必要な場面¶

GeoJSON（WGS84）と国土地理院データ（JGD2011）を重ね合わせる → CRSを統一する
距離・面積を正確に計算したい → 平面直角座標系に変換する
Web地図サービスに表示したい → WGS84に戻す

2. GDALとOGRの概要¶

2.1 GDALとは¶

GDAL（Geospatial Data Abstraction Library）は、地理空間データを読み書きするためのオープンソースライブラリです。

GDAL：ラスターデータ（画像、標高データなど）を扱う
OGR：ベクターデータ（点・線・面のGIS データ）を扱う

実際にはどちらも同じパッケージに含まれており、コマンドラインからは以下のツールを使います。

コマンド	対象	主な用途
`ogrinfo`	ベクター	データの確認（CRS、属性、件数）
`ogr2ogr`	ベクター	フォーマット変換、座標変換、空間演算
`gdalinfo`	ラスター	メタデータ・CRS確認
`gdalwarp`	ラスター	投影変換、リサンプリング

2.2 インストールの確認¶

まず、GDALが使える状態か確認してください。

# GDALのバージョンを確認（これが動けばOK）
ogr2ogr --version

期待される出力例：

GDAL 3.7.3, released 2023/11/01

うまくいかないとき（command not found）: macOS では brew install gdal でインストールできます。Anacondaを使っている場合は conda install -c conda-forge gdal を試してください。それでも解決しない場合はメールで連絡してください。

3. ogr2ogrによる座標変換¶

3.1 ogrinfo：データの確認¶

変換の前に、元のデータのCRSと内容を確認します。-al はすべてのレイヤー、-so はサマリーのみ（全データは表示しない）の意味です。

# ファイルの概要を確認
ogrinfo -al -so data/raw/llm_log.geojson

出力例（このような情報が表示されます）：

INFO: Open of `data/raw/llm_log.geojson'
      using driver `GeoJSON' successful.

Layer name: llm_log
Geometry: Point
Feature Count: 4832
Extent: (140.985, 42.711) - (141.792, 43.387)
Layer SRS WKT:
GEOGCRS["WGS 84",
    ...
EPSG:4326

id: String (0.0)
timestamp: String (0.0)
category: String (0.0)
district: String (0.0)

確認すべき項目：

Geometry：データの種類（Point / Polygon / Line）
Feature Count：データの件数
EPSG:XXXX：現在のCRS
Extent：データが存在する範囲（この時点では「度」単位）

うまくいかないとき（ファイルが見つからない）: ERROR 4: No such file or directory が出た場合、pwd コマンドで現在地を確認し、cd ~/bigdata2/pipeline_practice で移動してから再実行してください。

3.2 ogr2ogr：フォーマット変換¶

# GeoJSON → GeoPackage（.gpkg）に変換（推奨形式）
ogr2ogr \
  -f "GPKG" \
  data/processed/llm_log.gpkg \
  data/raw/llm_log.geojson

# 変換できたか確認
ogrinfo -al -so data/processed/llm_log.gpkg

3.3 ogr2ogr：座標変換¶

ベクターデータのCRS変換は -t_srs（target SRS = 変換先のCRS）オプションで行います。

# WGS84（EPSG:4326）→ 平面直角座標系第XII系（EPSG:6680）に変換
ogr2ogr \
  -f "GeoJSON" \
  -t_srs EPSG:6680 \
  data/processed/llm_log_6680.geojson \
  data/raw/llm_log.geojson

変換できたか確認します。Extent（範囲）の単位が「度」から「メートル」に変わっているはずです。

ogrinfo -al -so data/processed/llm_log_6680.geojson
# Extent: (-120543.2, -27384.5) - (47231.8, 44892.1) のような値（メートル単位）

# WGS84（EPSG:4326）→ JGD2011（EPSG:6668）に変換
ogr2ogr \
  -f "GeoJSON" \
  -t_srs EPSG:6668 \
  data/processed/llm_log_jgd2011.geojson \
  data/raw/llm_log.geojson

元のCRSを明示する場合（ファイルにCRS情報がない場合）：

# -s_srs（source SRS）で元のCRSを指定する
ogr2ogr \
  -f "GeoJSON" \
  -s_srs EPSG:4326 \
  -t_srs EPSG:6680 \
  data/processed/llm_log_6680.geojson \
  data/raw/llm_log_no_crs.csv

うまくいかないとき（変換後のデータがおかしい）: Extentが極端な値になっている場合、-s_srs を明示していない可能性があります。入力ファイルにCRS情報が正しく入っているか ogrinfo -al -so で確認してください。

4. gdalinfo / gdalwarp によるラスターデータの操作¶

4.1 gdalinfo：ラスターデータの確認¶

# DEM（数値標高モデル）などのラスターデータの情報を確認
gdalinfo data/raw/dem_hokkaido.tif

出力例（抜粋）：

Driver: GTiff/GeoTIFF
Size is 3600, 2400
Coordinate System is:
PROJCRS["JGD2011 / Japan Plane Rectangular CS XII", ...]
  EPSG:6680
Origin = (141.000000000000000, 43.500000000000000)
Pixel Size = (0.000277777777778,-0.000277777777778)
Band 1 Block=3600x1 Type=Float32, ColorInterp=Gray
  Minimum=0.000, Maximum=1492.100, Mean=89.234, StdDev=142.531

確認すべき項目：

Coordinate System：CRS（EPSGコード）
Size：画像のピクセルサイズ（横×縦）
Pixel Size：1ピクセルが何メートル（または何度）に対応するか
Bandの統計値：標高データなら最小・最大・平均

4.2 gdalwarp：ラスターデータの投影変換¶

# EPSG:6680 → WGS84（EPSG:4326）に変換
gdalwarp \
  -s_srs EPSG:6680 \
  -t_srs EPSG:4326 \
  -r bilinear \
  data/raw/dem_hokkaido.tif \
  data/processed/dem_hokkaido_wgs84.tif

オプションの説明：

オプション	内容
`-s_srs EPSG:6680`	入力ファイルのCRS
`-t_srs EPSG:4326`	出力ファイルのCRS
`-r bilinear`	リサンプリング方法（bilinear=バイリニア補間。標高データに推奨）

# 変換結果を確認
gdalinfo data/processed/dem_hokkaido_wgs84.tif | grep "Coordinate System" -A 5

5. ハンズオン手順¶

所要時間の目安: 全手順で45〜75分程度です。授業時間内に終わらなければ持ち帰ってください。どこまで終わったかをメモしておくと再開しやすいです。

手順1: データのCRSを確認する¶

# 作業ディレクトリに移動
cd ~/bigdata2/pipeline_practice

# 出力ディレクトリを作成
mkdir -p data/processed

# 元データのCRS確認
ogrinfo -al -so data/raw/llm_log.geojson

以下をメモしておきましょう（演習課題で使います）：

Geometry の種類：_____
Feature Count：_____
EPSGコード：_____
Extent（範囲）：_____

手順2: WGS84から平面直角座標系に変換する¶

# WGS84 → 平面直角座標系第XII系（EPSG:6680）に変換
ogr2ogr \
  -f "GeoJSON" \
  -t_srs EPSG:6680 \
  data/processed/llm_log_6680.geojson \
  data/raw/llm_log.geojson

# 変換結果を確認（Extentの単位がメートルに変わっているか確認）
ogrinfo -al -so data/processed/llm_log_6680.geojson

手順3: 変換前後の座標を比較する¶

# 変換前の座標（度単位）
ogrinfo -al -q data/raw/llm_log.geojson | grep "POINT" | head -3

# 変換後の座標（メートル単位）
ogrinfo -al -q data/processed/llm_log_6680.geojson | grep "POINT" | head -3

手順4: Pythonで変換結果を読み込んで確認する¶

以下のコードを src/check_crs.py として保存し、python3 src/check_crs.py で実行してください。

# src/check_crs.py
import geopandas as gpd

# 変換前（WGS84）
gdf_wgs84 = gpd.read_file("data/raw/llm_log.geojson")
print("CRS（変換前）:", gdf_wgs84.crs)
print("座標の例（変換前）:", gdf_wgs84.geometry.iloc[0])
print()

# 変換後（平面直角座標系）
gdf_6680 = gpd.read_file("data/processed/llm_log_6680.geojson")
print("CRS（変換後）:", gdf_6680.crs)
print("座標の例（変換後）:", gdf_6680.geometry.iloc[0])

python3 src/check_crs.py

うまくいかないとき（ModuleNotFoundError: geopandas）: pip install geopandas または conda install -c conda-forge geopandas を実行してください。

手順5: 変換後のデータで距離を計算する（余裕があれば）¶

# src/calc_distance.py
import geopandas as gpd

# 平面直角座標系（単位: メートル）で距離を計算
gdf = gpd.read_file("data/processed/llm_log_6680.geojson")

p1 = gdf.geometry.iloc[0]
p2 = gdf.geometry.iloc[1]
distance_m = p1.distance(p2)
print(f"2地点間の距離: {distance_m:.1f} m（= {distance_m/1000:.2f} km）")

# WGS84での距離（参考：「度」単位なので実際の距離ではない）
gdf_wgs = gpd.read_file("data/raw/llm_log.geojson")
p1_wgs = gdf_wgs.geometry.iloc[0]
p2_wgs = gdf_wgs.geometry.iloc[1]
print(f"WGS84での'距離': {p1_wgs.distance(p2_wgs):.6f}°（これは距離ではない）")

6. 演習課題¶

締め切り：2週間後（次の次の授業の前日23:59）
提出方法：作成したファイルをzipにまとめてメールで送付
再提出：締め切り後も1週間以内であれば再提出可能

課題1（必須）¶

学内配布の llm_log.geojson について、ogrinfo で以下を確認し、結果をテキストファイル（results/ogrinfo_result.txt）に保存してください。

現在のCRS（EPSGコード）
データの件数
Extent（データの範囲）
属性列の名前と型

# 結果をファイルに保存する場合はリダイレクトを使う
ogrinfo -al -so data/raw/llm_log.geojson > results/ogrinfo_result.txt

課題2（必須）¶

ogr2ogr を使って以下の2種類の変換を実行し、変換後のCRSを ogrinfo で確認した結果もファイルに保存してください。

llm_log.geojson（WGS84）→ llm_log_jgd2011.geojson（JGD2011, EPSG:6668）
llm_log.geojson（WGS84）→ llm_log_6680.geojson（平面直角座標系XII系, EPSG:6680）

課題3（必須）¶

変換後のデータ（EPSG:6680）をPythonで読み込み、全データ点の重心（centroid）の座標を計算して表示してください。

import geopandas as gpd

gdf = gpd.read_file("data/processed/llm_log_6680.geojson")
centroid = gdf.geometry.unary_union.centroid
print(f"全データの重心: X={centroid.x:.1f}m, Y={centroid.y:.1f}m")

課題4（任意・加点対象）¶

重心座標をWGS84に逆変換して、緯度・経度で表示してください。

まとめ¶

項目	内容
地理座標系	緯度・経度で表す（単位：度）。距離計算に不向き
投影座標系	平面に投影して表す（単位：m）。距離・面積計算に適する
WGS84（EPSG:4326）	GPS・GeoJSONの標準。世界規模のデータに使用
JGD2011（EPSG:6668）	日本の標準測地系。国土地理院データに多い
平面直角座標系XII系（EPSG:6680）	北海道用。距離・面積計算に使用
ogrinfo	ベクターデータのCRS・属性・件数を確認するコマンド
ogr2ogr	ベクターデータのフォーマット変換・CRS変換コマンド
gdalinfo	ラスターデータのメタデータ・CRS確認コマンド
gdalwarp	ラスターデータの投影変換コマンド

次回（第03回）は、変換後の地点データとポリゴンデータを空間結合し、区ごとの集計を行います。

ビッグデータ応用Ⅱ〔特徴量エンジニアリング〕 / 担当：小野原彩香

第02回 座標変換と投影変換¶