地理座標系（投影法）のまとめと，緯度経度格子面積マップ作成

[2019/08/22] pyprojの新旧versionについての記述を追加

座標系まとめ

Wikiの一覧が非常に見やすいほか，ESRIの解説がわかりやすい．その他，この記事やこのまとめなども大変参考になる．

地理座標系

「緯度経度（つまり角度）で表現される座標系」のことを地理座標系と呼ぶ。地球をどのような立体で近似しているか，もしくは座標系の原点（地球の重心の位置）をどこにとるのか（＝測地系）によって同じ緯度経度でも位置がずれたりする．測地系も世界系・日本系それぞれでいくつか存在するが，世界系ではWGS84楕円体という1984年に策定された基準が有名である（WGS84楕円体の地理座標系はEPSG:4326）．その一方で，以後に示す座標系はすべて，投影座標系と分類される．

正距円筒（World Equidistant Cylindrical）図法

普段よく目にする地図の投影法で，「緯度経度が正方格子をなすように作られたもの」である．南北方向には正距．東西方向には緯度が大きくなるほど引き延ばされる．WGS84楕円体に対して投影したもののEPSGは32663．

球面メルカトル（Spherical Mercator）図法

「正距円筒図法の南北方向を引き延ばして，正角※としたもの」である．もしくは，「球体の地球を南北方向にのびる円筒で包み，地球中心に置いた点光源が円筒上に作る影で投影したもの」である．緯度が高くなればなるほどグリッド幅が広くなり，極ではその長さが無限大となる（80°Sから84°Nまでの使用が推奨されている）．Google Mapをはじめとするスクロール地図はWGS 84/Pseudo-Mercator（EPSG:3857）と呼ばれる球面メルカトル図法の一つであり，WGS84楕円体に対するもの（EPSG:3395）とは少し異なるよう．いずれにせよ，南北に85.06°までの緯度での使用が推奨される点に注意が必要.狭いスケールの描画であれば，正角地図としての使用価値が高い．
※正角：等角写像．緯線経線がどこでも直交，どこに座標をとってもそこから等距離の点を結んだ軌跡が円で表示される←参考：テイソーの指示楕円

ランベルト正積円筒（Lambert Cylindrical Equal-Area）図法

EPSG:9834/9835など（後述）で定義され「正距円筒図法の南北方向を圧縮して，正積としたもの」であり，アフリカの面積が想像以上に大きいことを再確認できる図法．
※正積：どこに座標をとってもそこから一定距離の範囲を塗りつぶした領域が等面積で表示される←参考：テイソーの指示楕円

ランベルト正積方位（Lambert azimuthal equal-area）図法

EPSG:1027/9820などで定義され「（地図中央からの方位と距離が正確な，全体が円形の地図をつくる）正距方位図法の半径方向を圧縮して，正積としたもの」である．

正弦曲線（Sinusoidal）図法

一見ぐにゃっとして見栄えがよろしくないが，「緯線を線分で表示し，その長さを緯線円周（高緯度ほど短い）と対応させることで，経線が正弦曲線の山と相似になるようにしたもの」であるときけば理にかなった図法とわかる．緯度方向に正距であり，正距円筒図法の経度方向引き延ばしバージョンともいえる．また，正積図法の1つでもある．NASAの人工衛星センサで非常に有名なMODISの画像の配布にはこの投影法が用いられている．ESRI:54008などで定義される．

ユニバーサル横メルカトル（UTM）図法

横メルカトル図法は，「球体の地球を横方向にのびる円筒で包み，地球中心に置いた点光源が円筒上に作る影で投影したもの」である（球面メルカトル図法の円筒の向きを変えただけ）．円筒軸と直交する面を通る（包んだ円筒との接線となる）経線のうち，地図の中心になる方（半円周分）を中央経線と呼ぶ．球面メルカトル図法では真横真縦に経緯線が投影された．それに対し，横メルカトル図法で真横真縦に引かれるのはそれぞれ赤道と中央経線のみ．他は歪んでいる（参考リンク）．横メルカトル図法における東西方向，南北方向は，真横方向，真縦方向を意味しない．ただし，球面メルカトル図法と同様に正角図法であるため，歪んだ経緯線どうしは直交する．
UTM直交座標系は，中央経線を変えて描画した地図を細長い長方形に切断したものの集合である．縦は，中央経線上における南緯80°から北緯84°までの間．横は，赤道上における中心経線±3°の間．中央経線を6°ごとに変化させ，60枚の長方形地図が作成される．
先述の通り，赤道と中央経線を除く経緯線は長方形の縦横と平行とならない（ほぼほぼ平行となるように6°ごとに切り取っている）．そのため，長方形上に新たな直交座標系（単位m）を定義し，使いやすくしたものがUTM直交座標系である．
UTMのEPSGは北半球と南半球における60ゾーンそれぞれに対して指定される（北半球:326XX，南半球:327XX (XX:01~60)）．

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

正方格子で区切った全球の画像データを長方形の地図に描画すると，緯度経度で囲まれた面積は同じようにみえる．しかし，特に高緯度地域では東西方向を引き延ばして描画される（つまり，高緯度地域のグリッド実面積は小さい）．この辺を適切に扱わないと，領域平均の計算などで不都合が生じる．

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

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

コード

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

（追記）pyprojのversion（pyproj.__version__で確認可能）によって文法が異なる．version 1.9.xおよび最新の2.2.xについてまとめた．

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(geolist): # shape()を用いて緯度経度のリストからポリゴンを定義
    geom = shape({'type': 'Polygon',
        'coordinates': [geolist]})
    return geom

def geom_to_area(geom, epsg):
    '''
    ポリゴンの座標系を変換し，面積の出力
    '''
    project = partial( # 関数の入力を部分的に与える．
        pyproj.transform, 
        pyproj.Proj(init='epsg:4326'), # WGS84     # ver 1.9.x
        pyproj.Proj(init='epsg:%i' %epsg))         # ver 1.9.x
        # 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

geolistは，ポリゴンを形成する頂点座標を順番に並べたリストである．今回はそれが緯度経度方向の長方形となる．例えば，次のコードによって赤道真上の正方格子の面積を求めることができる．

test_area

geolist = [[0.,0.],[0.,1.],[1.,1.],[1.,0.],[0.,0.]]
geom = list_to_geom(geolist)

print(geom.area)

結果は1.0と出力される．1°格子の正方形の面積を単位°で計算しただけである．
メートル単位での面積の算出のために，geom_to_area(geom,epsg)で座標系を指定した座標系に変換する．

上で定義した関数を用いて，1°格子のマップの面積を次のように計算させる（単位はm）．

import numpy as np

def area_map_deg(lat_deg, lon_deg, epsg):
    # 最南西点の座標と解像度の定義（アフィン変換係数の表記に合わせている）．
    ll_x, h_res, _, ll_y, _, v_res = [
    -180, lon_deg, 0, -90, 0, lat_deg]
    print('ll_x:',ll_x)
    print('ll_y:',ll_y)
    print('h_res:',h_res)
    print('v_res:',v_res)

    nhpx = int(360/lon_deg)
    nvpx = int(180/lat_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]]   # ver 1.9.x
            # geolist = [
            #  [btm, lft],
            #  [top, lft],
            #  [top, rgt],
            #  [btm, rgt],
            #  [btm, lft]] # ver 2.2.x
            geom = list_to_geom(geolist)      # ポリゴンを定義
            area = geom_to_area(geom, epsg)   # 面積の計算
            outar[v, h] = area # 面積の代入
    return outar

# 関数の実行（numpy配列を取得）
lat_deg  = 1.
lon_deg  = 1.
epsg = 3857 # <- 結果的に3410がいいという結論になる。
result = area_map_deg(lat_deg, lon_deg, epsg)

# 結果の保存
outputfile = '%02ix%02ideg_area.npy'%(lat_deg*10, lon_deg*10)
np.save(outputfile, result)

結果　～面積計算における座標系選択の重要性～

まず，このサイトを参考に，座標系を球面メルカトル図法（epsg=3857）に変換する．これによって，緯度経度直交座標系かつ単位がメートルの情報に緯度経度情報を変換して面積を計算することができる．結果を簡単に描画すると次のようになる（km2単位）．

赤道付近のグリッドはおおよそ111km×111km程度となるため，おそらくうまくいっている．しかし，高緯度地域ほど格子面積が大きくなっているのは明らかにおかしい（予想と真逆である）．

座標系を正距円筒図法（epsg=32662）に変更してみる．こちらも緯度経度直交座標系かつ単位がメートルの投影法である．結果を下に描画する．

全グリッドの面積がほぼ全く同じ値となってしまった．

何が起きている

この辺りで理解が進んでいく．まず前提として，座標変換はうまくいっている（コーディングミスはない）．しかし，「単位メートルというのはあくまでその投影法内で定義されている単位である」という点に注意しなければならない．

正距円筒図法では，緯度の増加に伴って東西方向を引き延ばし，東西方向には「経度が1°変化したときの距離（メートル）が，どの緯度においても赤道上での変化と同じ値となる」ようにメートルを定義している．その結果，各格子の面積が等しくなる．

対して，球面メルカトル図法では，緯度の増加に伴って東西方向を引き延ばし，正角となるよう同時に南北方向も引き延ばしている．その結果，南北方向には「東西方向に引き延ばした比率に合わせて南北方向にも引き延ばし，東西方向で定義した（正距円筒図法の定義に同じ）メートルに合わせて縦横方向の長さの比率がそろう」ようにメートルを定義している．その結果，高緯度になるほど格子面積が南北長さの増加の分だけ大きくなる．

つまり．．．

正積の図法を使わないと正しい面積の計算ができないということが分かる．
そこで，ランベルト正積円筒（epsg=9835）を試みる．すると，次のエラーが返る．

# ver 1.9.x
  File "_proj.pyx", line 84, in _proj.Proj.__cinit__ (_proj.c:1170)
RuntimeError: no options found in 'init' file

上のコードでpyproj.Proj(init='XXXX')としていたところで，initに使用できる（登録されている）座標系じゃないよ，というメッセージ．一覧はここから見れる（注：旧versionで用いられているものであり，新しいpyprojでは該当するファイルが別形式っぽい）．確かに9835は存在しない．initにはesri（リンク）なども使用できる（注：同上）．

より調べてみると，どうやらepsg=3410で定義されるNSIDC EASE-Grid Globalという投影法がepsgで使用できる唯一の全球正積地図のようである．というわけでこれで実行．結果は次のようになった（km2単位）．

これで成功っぽい．
ちなみに，正弦曲線（Sinusoidal）図法についても同様に計算を行い，ほぼ全く同じ結果が得られた．とにかく，面積の計算には，正積な投影法を用いなくてはいけないということが分かった．

        pyproj.Proj(init='esri:54008')  # 1.9.x
        pyproj.CRS("ESRI:54008")        # 2.2.x (←geolistは1.9.xのものを使用)

どうも，新旧versionでほとんどの投影法に関して緯度経度の定義順序を逆転させたらしい．ESRIのものもそのうち置き換わると考えられる．