Pythonでの測地線のバッファリング

Question

土地ポリゴンをShapelyMultiPolygonと指定すると、たとえば、海岸線の周りに12海里のバッファ。

ユークリッド計算を使用しているため、Shapely bufferメソッドを使用すると機能しません。

Pythonで測地線バッファを計算する方法を教えてもらえますか？

Answer 1

これは、ライブラリが平面計算用であることを明示的にドキュメントに明示しているので、滑らかな問題ではありません。それにもかかわらず、あなたの質問に答えるには、マルチポリゴンに使用している座標系を指定する必要があります。 WGS84投影（lat、lon）を使用していると仮定すると、これは別のSO質問（fix-up-shapely-polygon-object-when-discontinuous-after-map-projection）で見つけたレシピです。ライブラリpyprojが必要です。

import pyproj 
from shapely.geometry import MultiPolygon, Polygon 
from shapely.ops import transform as sh_transform 
from functools import partial 

wgs84_globe = pyproj.Proj(proj='latlong', ellps='WGS84') 

def pol_buff_on_globe(pol, radius): 
    _lon, _lat = pol.centroid.coords[0] 
    aeqd = pyproj.Proj(proj='aeqd', ellps='WGS84', datum='WGS84', 
         lat_0=_lat, lon_0=_lon) 
    project_pol = sh_transform(partial(pyproj.transform, wgs84_globe, aeqd), pol) 
    return sh_transform(partial(pyproj.transform, aeqd, wgs84_globe), 
          project_pol.buffer(radius)) 

def multipol_buff_on_globe(multipol, radius): 
    return MultiPolygon([pol_buff_on_globe(g, radius) for g in multipol]) 

pol_buff_on_globe機能は次のとおりです。まず、ポリゴン重心を中心とした方位角の等距離投影を作成します。次に、ポリゴンの座標系をその投影に変更します。その後、そこにバッファを作成し、バッファされたポリゴンの座標系をWGS84座標系に変更します。

いくつかの特別な注意が必要です：

• あなたはaeqd投影に使用される距離にしたい距離を変換する方法を見つける必要があります。
• 極を含めてバッファリングしないように注意してください（上記のSOの質問を参照してください）。
• ポリゴンの重心を使用して投影を中心に置いているという事実は、答えが十分であることを保証するはずですが、精度要件がある場合は、この解決策を使用しないでください。使用している典型的なポリゴンです。