GIS × Python Tutorial 3.0 ~ pyproj ~
はじめに
この記事では python での pyproj の使用方法についてを解説していきます。
pyproj とは
pyproj は地理空間の変換に使用されるライブラリで、このライブラリを扱う事で異なる座標系間での座標変換を行う事ができます。
座標系とは
地球は完全な球体ではありません。地球の形は「楕円体」であり、これを2次元上に投影する際に、経緯線網とよばれる経度と緯度の線を用いて平面に投影しています。
地理座標系
地理座標系の原点は地球の重心にあり、その重心からの角度で座標を表現するものです。通常よく目にする経緯度（経度: 140.23567, 緯度: 39.00345）などがこれにあたり、経度は東西方向、緯度は南北方向を示します。地球上のどこにいても座標を特定する事が可能な、非常に広域な座標系です。 しかしデメリットとして距離や面積を正確に求めるのが難しくなります。
有名な地理座標系には、GPSに使用されているWGS84（World Geodetic System 1984 : 世界測地系1984）があります。
投影座標系
投影座標系は地球の3次元形状を2次元平面に投影して、XY座標で表現する座標系です。投影座標系の原点は様々な場所にあり、データの範囲などによってその原点を選択する事が出来ます。これは原点からの距離（メートルなど）で座標を表現し、距離や面積を正確に求める事が出来ます。地球は平面ではありませんので、目的地に近い原点を持つ投影座標系を選択する事で、より正確な距離や面積を求める事が出来ます。
日本国内で使用する場合主に2つの座標系から選択する事が多いです。
UTM座標系
UTM座標系（Universal Transverse Mercator）は、地図投影法の一種で、南北方向の精度を向上させた座標系です。地球は60の南北のゾーンに分割され、各ゾーンは経度で6°単位の幅を持ちます。また、日本はゾーン51～56までの座標系を使用しています。
平面直角座標系
平面直角座標系は、狭い範囲を対象とした測量や大縮尺地図に使われる座標系です。地球自体は楕円ですが、範囲を狭くする事で平面上に投影し、計算を簡単にしています。平面直角座標系は日本を19のゾーンに分割し、各ゾーンに原点を設けています。
都道府県を跨ぐような大きな範囲のデータセットであればUTM座標系を選択しますが、もっと狭い範囲のデータであれば平面直角座標系を使用すると良いでしょう。
平面直角座標系の注意点
少しややこしいですが多くの人がイメージする数学のXY座標（横軸がXで縦軸がY）と異なり、横軸がYで縦軸Xとなります。国土交通省では以下の様に定義されています。
座標系のＸ軸は、座標系原点において子午線に一致する軸とし、真北に向う値を正とし、座標系のＹ軸は、座標系原点において座標系のＸ軸に直交する軸とし、真東に向う値を正とする。
つまり北方向は+Y、南方向は-Y、東方向は+X、西方向は-Yとなります。
実行
from pprint import pprint
import pyproj
IN_EPSG = 'EPSG:4326'
OUT_EPSG = 'EPSG:6678'
WKT -CRS (Well-Known text Representation of Coordinate Reference System) の取得
WKTとはベクター形式のデータを表現しているマークアップ言語の事ですが、CRSを表現することもできます。
pprint(
pyproj
.CRS(IN_EPSG)
.to_wkt(pretty=True)
)
('GEOGCRS["WGS 84",\n'
' ENSEMBLE["World Geodetic System 1984 ensemble",\n'
' MEMBER["World Geodetic System 1984 (Transit)"],\n'
' MEMBER["World Geodetic System 1984 (G730)"],\n'
' MEMBER["World Geodetic System 1984 (G873)"],\n'
' MEMBER["World Geodetic System 1984 (G1150)"],\n'
' MEMBER["World Geodetic System 1984 (G1674)"],\n'
' MEMBER["World Geodetic System 1984 (G1762)"],\n'
' MEMBER["World Geodetic System 1984 (G2139)"],\n'
' ELLIPSOID["WGS 84",6378137,298.257223563,\n'
' LENGTHUNIT["metre",1]],\n'
' ENSEMBLEACCURACY[2.0]],\n'
' PRIMEM["Greenwich",0,\n'
' ANGLEUNIT["degree",0.0174532925199433]],\n'
' CS[ellipsoidal,2],\n'
' AXIS["geodetic latitude (Lat)",north,\n'
' ORDER[1],\n'
' ANGLEUNIT["degree",0.0174532925199433]],\n'
' AXIS["geodetic longitude (Lon)",east,\n'
' ORDER[2],\n'
' ANGLEUNIT["degree",0.0174532925199433]],\n'
' USAGE[\n'
' SCOPE["Horizontal component of 3D system."],\n'
' AREA["World."],\n'
' BBOX[-90,-180,90,180]],\n'
' ID["EPSG",4326]]')
CF（Climate and Forecast）の取得
pprint(
pyproj
.CRS(IN_EPSG)
.to_cf()
)
{'crs_wkt': 'GEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 '
'ensemble",MEMBER["World Geodetic System 1984 '
'(Transit)"],MEMBER["World Geodetic System 1984 '
'(G730)"],MEMBER["World Geodetic System 1984 '
'(G873)"],MEMBER["World Geodetic System 1984 '
'(G1150)"],MEMBER["World Geodetic System 1984 '
'(G1674)"],MEMBER["World Geodetic System 1984 '
'(G1762)"],MEMBER["World Geodetic System 1984 '
'(G2139)"],ELLIPSOID["WGS '
'84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic '
'latitude '
'(Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic '
'longitude '
'(Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],USAGE[SCOPE["Horizontal '
'component of 3D '
'system."],AREA["World."],BBOX[-90,-180,90,180]],ID["EPSG",4326]]',
'geographic_crs_name': 'WGS 84',
'grid_mapping_name': 'latitude_longitude',
'horizontal_datum_name': 'World Geodetic System 1984 ensemble',
'inverse_flattening': 298.257223563,
'longitude_of_prime_meridian': 0.0,
'prime_meridian_name': 'Greenwich',
'reference_ellipsoid_name': 'WGS 84',
'semi_major_axis': 6378137.0,
'semi_minor_axis': 6356752.314245179}
CRSオブジェクトからEPSGコードを取得する
crs = pyproj.CRS(IN_EPSG)
print(F"EPSG: {crs.to_epsg()}")
EPSG: 4326
投影変換
実際に投影変換してみましょう。
lon = 140.46548476359524
lat = 40.60794254882885
transformer = (
pyproj
.Transformer
.from_crs(
pyproj.CRS(IN_EPSG),
pyproj.CRS(OUT_EPSG),
always_xy=True
)
)
print('<<< Transformer >>>')
pprint(transformer)
x, y = transformer.transform(lon, lat)
print(f"""
{IN_EPSG} : {(lon, lat)}
{OUT_EPSG}: {(x, y)}
""")
<<< Transformer >>>
<Concatenated Operation Transformer: pipeline>
Description: axis order change (2D) + Inverse of JGD2011 to WGS 84 (1) + Japan Plane Rectangular CS zone X + axis order change (2D)
Area of Use:
- name: Japan - onshore and offshore.
- bounds: (122.38, 17.09, 157.65, 46.05)
EPSG:4326 : (140.46548476359524, 40.60794254882885)
EPSG:6678: (-31128.58688942413, 67564.52258165448) 5
複数の座標を変換する場合も同様に変換する事が出来ます。
lons = [lon, lon + 0.001]
lats = [lat, lat + 0.001]
xs, ys = transformer.transform(lons, lats)
print(f"""
{IN_EPSG} : {(lons, lats)}
{OUT_EPSG}: {(xs, ys)}
""")
EPSG:4326 : ([140.46548476359524, 140.46648476359525], [40.60794254882885, 40.60894254882885])
EPSG:6678: ([-31128.58688942413, -31043.500596195187], [67564.52258165448, 67675.20503145826])
おわりに
今回は pyproj での投影変換について簡単に解説しました。この記事を作成する際に使用したNotebookはGitHub上で公開しています。
Discussion