日本政府統計の「地図で見る統計（統計GIS）で使われているGISデータで、同じデータを異なるフォーマットでダウンロードできるようになっている。どこがどう違うのか調べてみた。
shpは昔からあるバイナリフォーマットで、GMLは最近整備されてきたXMLフォーマット。python geopandas で読み込む（内部的には fiona）場合は、shp を使う。fiona が内部で使っている GDAL で、必要な XML 仕様を完全にサポートしていないかららしい。
後は projection の差で、proj4 な定義を見るとこうなっていた（geopandasの解説）。神戸市で取り出していて、平面直角座標系は地域によって少し違うはず。
| 呼称 | proj4文字列 |
|---|---|
| 日本測地系緯度経度 | +ellps=bessel +no_defs +proj=longlat |
| 世界測地系緯度経度 | +ellps=GRS80 +no_defs +proj=longlat |
| 日本測地系平面直角座標系 | +lon_0=134.3333333333333 +y_0=0 +lat_0=36 +ellps=bessel +no_defs +units=m +proj=tmerc +x_0=0 +k=0.9999 |
| 世界測地系平面直角座標系 | +lon_0=134.3333333333333 +y_0=0 +lat_0=36 +ellps=GRS80 +no_defs +units=m +proj=tmerc +x_0=0 +k=0.9999 |
平面直角座標系は、地面と同じように平面に見える数字になる。unitsがmになってますね。平面直角座標系で、日本測地系と世界測地系は楕円球体が違うということですかね。日本測定系と呼ばれているものはベッセル楕円体ということで、つまり旧日本測地系。いまどき WGS84 なデータがたくさんあるので、日本測地系 2000 な「世界測地系」の GRS80 を使うのが良さそう。
いろいろ計算したい場合は、「世界測地系緯度経度」をダウンロードして計算しておいて、最後に図にするときに平面直角座標系に変換するといいかな。to_crs() で一発でできるし。
data.to_crs("+lon_0=134.3333333333333 +y_0=0 +lat_0=36 +ellps=GRS80 +no_defs +units=m +proj=tmerc +x_0=0 +k=0.9999").plot()
面積を使って値を割り振るとかする時は、平面直角座標系で計算しないとダメか。