Python備忘録 国土数値情報「京都府」を題材にGeopandasの基本（１）（https://shikuuk.blogspot.com/2023/01/pythongeopandas.html）で未解決のままの面積計算を学習しました。
面積計算にあたっては「Map projection（地図投影）」と「Coording Reference System（座標系）」を理解する必要があるようです。
難しいことは脇に置き、まずは（１）でも使用した「国土数値情報（https://nlftp.mlit.go.jp/ksj/index.html）」にある「政策区域－行政区域」情報（https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N03-v3_1.html）から「京都府」のShapeFileデータを読み込みます。
以下、出力です。
| N03_001 | N03_002 | N03_003 | N03_004 | N03_007 | geometry | |
|---|---|---|---|---|---|---|
| 0 | 京都府 | None | None | None | None | MULTIPOLYGON (((135.31410 35.45566, 135.31255 ... |
| 1 | 京都府 | None | 京都市 | 北区 | 26101 | POLYGON ((135.72539 35.17080, 135.72552 35.170... |
| 2 | 京都府 | None | 京都市 | 上京区 | 26102 | POLYGON ((135.75062 35.03822, 135.75079 35.038... |
| 3 | 京都府 | None | 京都市 | 左京区 | 26103 | POLYGON ((135.80481 35.31708, 135.80586 35.316... |
| 4 | 京都府 | None | 京都市 | 中京区 | 26104 | POLYGON ((135.73195 35.02251, 135.73195 35.022... |
以下のコードで座標系を確認します。
以下、出力です。
<Geographic 2D CRS: EPSG:6668> Name: JGD2011 Axis Info [ellipsoidal]: - Lat[north]: Geodetic latitude (degree) - Lon[east]: Geodetic longitude (degree) Area of Use: - name: Japan - onshore and offshore. - bounds: (122.38, 17.09, 157.65, 46.05) Datum: Japanese Geodetic Datum 2011 - Ellipsoid: GRS 1980 - Prime Meridian: Greenwich
難しいことは脇に置きますが、座標系の名前が「Name: JGD2011」と表示されています。JGD2011はネットで検索すると、世界測地系に準拠した日本の座標系といった意味のことが書かれていると思います。
あと大切そうなのは「EPSG:6668」です。EPSGは「European Petroleum Survey Group」の団体名の略語で、GISで使用される様々な座標系に数値コードが割り振られており、ここから先の座標変換の主要なパラメータのようです。
このGeoDataFrame「kyoto」の座標系をHelsinkiUniv.のAutomating GIS Processingに従って「Lambert Azimuthal Equal-Area projection」に変換してみます。
「Lambert Azimuthal Equal-Area projection」のEPSGコードは3035みたいです。
変換コードは以下のようになります。おそろしく簡単です。こんなに簡単に座標系の変換ができるというのは、昭和のおじさんにとってはちょー驚きです。昔は座標系の定義から２つの座標系の変換式を作成して、それをコーディングしないといけなかったのではないかと愚考します。
以下、出力です。
| N03_001 | N03_002 | N03_003 | N03_004 | N03_007 | geometry | |
|---|---|---|---|---|---|---|
| 0 | 京都府 | NaN | NaN | NaN | NaN | MULTIPOLYGON (((9882965.159 9287745.942, 98830... |
| 1 | 京都府 | NaN | 京都市 | 北区 | 26101 | POLYGON ((9891118.447 9327971.717, 9891111.857... |
| 2 | 京都府 | NaN | 京都市 | 上京区 | 26102 | POLYGON ((9903549.650 9330896.818, 9903541.942... |
| 3 | 京都府 | NaN | 京都市 | 左京区 | 26103 | POLYGON ((9871748.858 9334926.407, 9871721.335... |
| 4 | 京都府 | NaN | 京都市 | 中京区 | 26104 | POLYGON ((9906166.009 9329171.029, 9906171.630... |
geometryの列が緯度経度から変換されているのがわかりました。
このGeoDataFrame'kyoto_EPSG3035'を図化してみます。わかりやすくするため、JGD2011座標系の図化を左側に置いてみます。
何じゃこりゃ！左と右で全然違う絵になりました。よくよく見ると、左側のJGD2011座標系の絵が約180°回転しているようにも見えます。フィンランドと日本は地球のほぼ反対側にあるので、このように見えるのでしょうか？
EPSG3035の座標系を利用して、行(row)ごとの面積を計算してみます。
結果は以下の通りです。
0 4.612204e+09
1 9.488360e+07
2 7.028689e+06
3 2.467674e+08
4 7.409009e+06
...
995 3.359364e+01
996 5.900918e+01
997 6.089839e+07
998 6.489418e+01
999 1.083829e+08
Length: 1000, dtype: float64index=0は京都府全体のデータを表しており、4,612,204,000m2、すなわち、4,612km2となっています。京都府のHP(https://www.pref.kyoto.jp/intro/)によれば、京都府の面積は4,613.2km2で面積は当たらずとも遠からずのようです。index=1は京都市の北区です。京都市のHP(https://www.city.kyoto.lg.jp/sogo/page/0000015581.html)によると北区の面積は94.88km2となっており、94,883,600m2=94.88km2とこちらもほぼ同等です。
フィンランドの座標系で日本国内の面積がほぼほぼ計算できるというのは少し驚きました。たまたまなのかもしれません。
次は日本でよく用いられる座標系を使ってみました。
まずはUTM座標系(EPSG6690)を使ってみます。UTM座標については例えば国土地理院のこちらのページ(https://www.gsi.go.jp/chubu/minichishiki10.html)などが参考になります。
| N03_001 | N03_002 | N03_003 | N03_004 | N03_007 | geometry | |
|---|---|---|---|---|---|---|
| 0 | 京都府 | NaN | NaN | NaN | NaN | MULTIPOLYGON (((528502.344 3923621.755, 528362... |
| 1 | 京都府 | NaN | 京都市 | 北区 | 26101 | POLYGON ((566055.565 3892225.761, 566067.559 3... |
| 2 | 京都府 | NaN | 京都市 | 上京区 | 26102 | POLYGON ((568464.281 3877538.646, 568479.746 3... |
| 3 | 京都府 | NaN | 京都市 | 左京区 | 26103 | POLYGON ((573157.049 3908504.365, 573252.560 3... |
| 4 | 京都府 | NaN | 京都市 | 中京区 | 26104 | POLYGON ((566773.853 3875783.757, 566774.152 3... |
このGeodataFrameを図化します。こちらもわかりやすいようにJGD2011(EPSG6668)と並べて表示します。
結果はこちらでございます。さすが、日本でよく利用される座標系です。JGD2011とほぼ相似形の地図になっています。
0 4.608775e+09
1 9.481760e+07
2 7.023885e+06
3 2.466021e+08
4 7.403931e+06
...
995 3.356714e+01
996 5.896263e+01
997 6.085044e+07
998 6.484298e+01
999 1.082965e+08
Length: 1000, dtype: float64index=0（京都府全体）は4,609km2[京都府HP(https://www.pref.kyoto.jp/intro/)によれば、京都府の面積は4,613.2km2]、index=1（京都市北区）は94.82km2[京都市HP(https://www.city.kyoto.lg.jp/sogo/page/0000015581.html)は94.88km2]で微妙な誤差。
次に平面直角座標系（詳しくは、国土地理院HP(https://www.gsi.go.jp/sokuchikijun/datum-main.html)などを参照してください）を使用します。全国を19の平面直角座標系に分割して管理しており、京都府は第6系に位置します。平面直角座標系日本第6系(EPSG:6674)に変換します。
| N03_001 | N03_002 | N03_003 | N03_004 | N03_007 | geometry | |
|---|---|---|---|---|---|---|
| 0 | 京都府 | NaN | NaN | NaN | NaN | MULTIPOLYGON (((-62259.781 -60173.959, -62400.... |
| 1 | 京都府 | NaN | 京都市 | 北区 | 26101 | POLYGON ((-25014.199 -91956.629, -25002.206 -9... |
| 2 | 京都府 | NaN | 京都市 | 上京区 | 26102 | POLYGON ((-22752.451 -106670.878, -22736.995 -... |
| 3 | 京都府 | NaN | 京都市 | 左京区 | 26103 | POLYGON ((-17747.335 -75746.387, -17652.098 -7... |
| 4 | 京都府 | NaN | 京都市 | 中京区 | 26104 | POLYGON ((-24460.803 -108409.179, -24460.566 -... |
面積を計算します。
0 4.611648e+09
1 9.486623e+07
2 7.027372e+06
3 2.467202e+08
4 7.407626e+06
...
995 3.359090e+01
996 5.900436e+01
997 6.089317e+07
998 6.488888e+01
999 1.083790e+08
Length: 1000, dtype: float64index=0（京都府全体）は4,612km2[京都府HP(https://www.pref.kyoto.jp/intro/)によれば、京都府の面積は4,613.2km2]、index=1（京都市北区）は94.87km2[京都市HP(https://www.city.kyoto.lg.jp/sogo/page/0000015581.html)は94.88km2]でこちらも先ほどよりはマシです微妙な誤差。
結局、地図から求める面積は座標系によって微妙に変わるから、真値というのはないのかなぁと思った今日この頃です。
最後に平面直角座標系（日本第6系）に変換した結果を図化しておきます。当然のことながら、JGD2011、UTM座標系と相似形の図になります。
コメント
コメントを投稿