Python
GIS
shapely
pyproj

緯線経線格子面積マップ作成

座標系まとめ

緯度経度座標系(世界測地系WGS84)

EPSG:4326で定義され,「球体上の任意の座標点を緯度経度と対応付けるもの」である(測地座標系と分類される).その一方で,以後に示す座標系はすべて,投影座標系と分類される.

正距円筒図法

普段よく目にする地図の投影法で,「緯度経度が正方格子をなすように作られたもの」である.南北方向には正距.東西方向には緯度が大きくなるほど引き延ばされる.

球面メルカトル図法

EPSG:3857で定義され,「正距円筒図法の南北方向を引き延ばして,正角(等角写像)としたもの」である.緯度が高くなればなるほどグリッド幅が広くなり,極ではその長さが無限大となる.Google Mapをはじめとするスクロール地図で用いられ,その場合(Web Mercator),南北に85度までの緯度でしか表示せず,それよりも緯度の大きい地点は捨象する.

ユニバーサル横メルカトル(UTM)図法

横メルカトル図法は,「球面メルカトル図法における円筒を横方向に適用したもの」である.南北方向に正距であり,東西方向には,円周に該当する経線から離れるほど引き延ばしが生じる.UTM図法は,「地球全体を経度6度ごと(60ゾーン)に分けて投影したもの」である.南緯80度から北緯84度までの間を投影の対象とする.UTMのEPSGは60ゾーンそれぞれに対して指定される.

以上のような詳細は,この記事このまとめなど.

緯度経度ラスタ(画像)データについて注意すべきこと

長方形の地図に書き起こせば緯度経度で囲まれた面積は同じようにみえるが,特に高緯度地域では面積を主に縦(緯度)方向に縮小して描画される(つまり,高緯度地域のグリッド実面積は大きい).この辺を適切に扱わないと,領域平均の計算などで不都合が生じる.

今回は,グリッドごとの面積をpythonコードで出力できるようにした.
本記事では,その一例として1度×1度の全球マップの各グリッドの面積を算出する.

コード

パッケージにはpyprojとshapelyを用いた.ベースとなる変換部分は次のようなコードになる(shapely公式の記述など参照).

list_to_geom_to_area
from functools import partial
import pyproj
from shapely.geometry import shape
from shapely.ops import transform

def list_to_geom(l): # shape()を用いて緯度経度のリストからポリゴンを定義
    geom = shape({'type': 'Polygon',
        'coordinates': [l]})
    return geom

def geom_to_area(geom): # ポリゴンの座標系を変換し,面積の出力
    project = partial( # 関数の入力を部分的に与える.
        pyproj.transform, 
        pyproj.Proj(init='epsg:4326'), # WGS84
        pyproj.Proj(init='epsg:3857')) # Mercator
    trans = transform(project, geom) # 座標系変換
    return trans.area # 面積を返す

lは,ポリゴンを形成する頂点座標を順番に並べたリストである.今回はそれが緯度経度方向の長方形となる.面積の算出のために,この長方形を再現する投影座標系である球面メルカトル図法(EPSG=3857)を用いている,これを指定しないと,おそらく頂点同士を結んだ実最小面積を算出する.この場合,緯度方向の最短距離は緯線に従わない.

上で定義した関数を用いて,1度格子のマップの面積を次のように計算させる(単位はm).

import numpy as np

deg=1.

def area_map_deg(deg):
    outputfile = '%02ideg_area.npy'%(deg*10)
    # 最南西点の座標と解像度の定義(アフィン変換係数の表記に合わせている).
    ll_x,h_res,_,ll_y,_,v_res = [
    -180,  deg,0, -90,0,  deg]
    print 'll_x:',ll_x
    print 'll_y:',ll_y
    print 'h_res:',h_res
    print 'v_res:',v_res

    nhpx = int(360/deg)
    nvpx = int(180/deg)

    # 結果格納用
    outar = np.zeros((nvpx,nhpx))

    for v in range(nvpx):
        btm = ll_y + v_res*v       # 下
        top = btm  + v_res         # 上
        for h in range(nhpx):
            lft = ll_x + h_res*h   # 左
            rgt = lft  + h_res     # 右
            # 座標のリスト(長方形,5点指定して必ず最後にパスを閉じる.)
            geolist = [
             [lft,btm],
             [lft,top],
             [rgt,top],
             [rgt,btm],
             [lft,btm]]
            geom = list_to_geom(geolist) # ポリゴンを定義
            area = geom_to_area(geom)    # 面積の計算
            outar[v,h] = area # 面積の代入
    # 結果の保存
    np.save(outputfile, outar)

結果を簡単に描画すると次のようになる(100km×100kmから1000km×1000kmまでのログスケール).赤道付近のグリッドはおおよそ111km×111km程度となるため,おそらくうまくいっている.
1deg.png