2017-06-16 29 views
-1

私が知っているように、織物はデカルト座標系のみを使用します。地球上に緯度と経度の2つのポイントがあります。この2つの点の周りに1km半径のバッファを作成し、このバッファが交差するポリゴンを見つける必要があります。地球投影上の2つの滑らかなポリゴンを交差させる

しかしconstrustion

バッファ=ポイント(54.4353,65.87343).buffer(0.001)は、単純な円を作成するが、地球上の投影では楕円となり、私は1キロ半径を有する2つの本当の円を必要とします。

私はバッファを新しい投影に変換してから交差させる必要があると思いますが、正しい方法はまだありません。

答えて

0

あなたが言うことをする必要があります。そのためには、投影を扱うライブラリを使用する必要があります(ここではpyprojが選択です)。機能point_buff_on_globeはあなたにその点を中心に正距方位図法で与えられた点をバッファリングした結果である緯度経度でのポリゴン(あなたができる最善を与えるGeodesic buffering in python

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

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

def point_buff_on_globe(lat, lon, radius): 
    #First, you build the Azimuthal Equidistant Projection centered in the 
    # point given by WGS84 lat, lon coordinates 
    aeqd = pyproj.Proj(proj='aeqd', ellps='WGS84', datum='WGS84', 
         lat_0=lat, lon_0=lon) 
    #You then transform the coordinates of that point in that projection 
    project_coords = pyproj.transform(wgs84_globe, aeqd, lon, lat) 
    # Build a shapely point with that coordinates and buffer it in the aeqd projection 
    aeqd_buffer = Point(project_coords).buffer(radius) 
    # Transform back to WGS84 each coordinate of the aeqd buffer. 
    # Notice the clever use of sh_transform with partial functor, this is 
    # something that I learned here in SO. A plain iteration in the coordinates 
    # will do the job too. 
    projected_pol = sh_transform(partial(pyproj.transform, aeqd, wgs84_globe), 
          aeqd_buffer) 
    return projected_pol 

で同様の質問があります要件2つの観測:。。

  1. 私はradius引数の単位を覚えていない私は、メートルであるので、あなたは10キロのバッファが必要な場合、あなたはそれが10E3渡す必要になると思う。しかし下さい。それを確認してください!
  2. または離れて離れている点を含むことができる。プロジェクションは、プロジェクションを中心にしているポイントにポイントが近くなるとうまく動作します。