Pythonを用いた緯度経度座標系ポリゴンデータの面積計算

はじめに

以前に投稿した記事において，緯度経度座標系のグリッド面積を求めるコードを紹介した．その際，pyprojの更新によって緯度と経度の定義順序が逆転していることについて触れた．この際は手作業で緯度と経度とを逆転させた．しかし，より複雑なポリゴンデータに対する適用のためにはこの操作を自動化させる必要がある．
本記事では，geopandasパッケージで公開している世界地図ポリゴンデータを題材として，各地域の面積を算出するコードを作成した．

使用モジュール

今回はpython3.7系でコードを作成した．pyprojのversionは2.2である．

import matplotlib.pyplot as plt
import geopandas as gpd
from functools import partial
import pyproj
from shapely.ops import transform

世界地図データの読込

デフォルトで様々な項目（人口，GDPなど）が入っているが，本記事ではgeometryのみを用いる．

world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
print(world.head())

ポリゴンデータから面積を算出するコード

（失敗例）緯度経度座標ポリゴンにareaメソッド適用

geometry列に並ぶshapely形式データに対し，geom.areaを適用することで，緯度経度座標系を用いて算出された面積を計算できる．しかし，これはkm2の単位での出力ではなく，高緯度ほど相対的に面積が過大評価となる．

world['geom.area'] = world['geometry'].area
print(world.head())

正積図法へ投影後に面積を計算

同じコードが以前紹介した記事にも記載してある．正積図法をみたすEPSGコードを指定している．

def geom_to_area(geom, epsg=3410): # epsg=3410：正積図法
    project = partial( # 関数の入力を部分的に与える．
        pyproj.transform, 
        pyproj.CRS.from_epsg(4326), # WGS84      # ver 2.2.x
        pyproj.CRS.from_epsg(epsg))              # ver 2.2.x
    trans = transform(project, geom) # 座標系変換
    return trans.area # 面積を返す

（失敗例）緯度経度を入れ替えずに面積を計算

試しに，テーブルの2行目にあるタンザニアの面積を算出してみる．テーブル内にあるgeometryをそのまま面積計算に用いると，下記のような結果になる．

tanz = world['geometry'][1] # Polygon
print(geom_to_area(tanz)

数字は計算で求められた面積を示しており，約77.4万km2となった．wikiによると，実際の値は約94.5万km2であった．境界線を簡略化しているとしても小さすぎる値である．
ここで，面積計算コード内で座標変換されたshapely形式のデータも取得してみると，上図に示すように，xyが反転してしまっていることが分かる．これが誤った面積を算出してしまう要因である．
解決策として，（やや奇妙ではあるが）事前にxy座標を反転させてから，正積図法への投影および面積計算を行うことがあげられる．

xy座標の反転

この記事を参照して，PolygonもしくはMultiPolygonのxy座標を反転するコードを作成した．

def Swap_xy(geom):
    # (x, y) -> (y, x)
    def swap_xy_coords(coords):
        for x, y in coords:
            yield (y, x)

    # if geom.type == 'Polygon':
    def swap_polygon(geom):
        ring = geom.exterior
        shell = type(ring)(list(swap_xy_coords(ring.coords)))
        holes = list(geom.interiors)
        for pos, ring in enumerate(holes):
            holes[pos] = type(ring)(list(swap_xy_coords(ring.coords)))
        return type(geom)(shell, holes)

    # if geom.type == 'MultiPolygon':
    def swap_multipolygon(geom):
        return type(geom)([swap_polygon(part) for part in geom.geoms])

    # Main
    if geom.type == 'Polygon':
        return swap_polygon(geom)
    elif geom.type == 'MultiPolygon':
        return swap_multipolygon(geom)
    else:
        raise TypeError('Unexpected geom.type:', geom.type)

MultiPolygonの例として，テーブルからカナダも取得して検証を行った．wikiの数字と比較して数%程度の誤差はあるものの，ポリゴンデータから面積の計算は適切に行えていることが分かる．

tanz = world['geometry'][1] # Polygon
cana = world['geometry'][3] # MultiPolygon

tanz_swap = Swap_xy(tanz)
cana_swap = Swap_xy(cana)

tanz_area = geom_to_area(tanz_swap) #  936857756916.421 (真値：94.5万km2)
cana_area = geom_to_area(cana_swap) # 9986763834110.908 (真値： 998万km2)

geopandasテーブルに面積計算関数の適用

上記で定義した関数を組み合わせて，面積計算関数を作成した．これを用いて，新たな列をテーブルに追加した．

def CalcArea(geom):
    geom_swap = Swap_xy(geom)
    area = geom_to_area(geom_swap)
    return area * 1e-6 # m2 -> km2 

world['area_km'] = world['geometry'].map(CalcArea)
print(world.head())

描画

新たに作成した列を用いて世界地図の描画．グリーンランドや南極大陸が，大きい国々と比較してそこまで大きくないことがわかる．

world.plot(figsize=(14, 7), column='area_km', legend=True)

重ねて，失敗例で示した，緯度経度座標系のまま計算してしまった面積も描画した．いわゆる，地図上での見た目の面積である．正しい地図とは大違いの結果となっていることが分かる．

おわりに

地域ごとの面積を簡単に比較できるサイトを紹介しておしまい．
http://thetruesize.com/