10

Python 3.4 と shapely 1.3.2 を使用して、長い/緯度の座標ペアのリストから Polygon オブジェクトを作成し、それらを解析するために既知のテキスト文字列に変換しています。このような Polygon は次のようになります。

POLYGON ((-116.904 43.371, -116.823 43.389, -116.895 43.407, -116.908 43.375, -116.904 43.371))

shapely は射影を処理せず、デカルト空間のすべてのジオメトリ オブジェクトを実装するため、次のようにそのポリゴンで area メソッドを呼び出します。

poly.area

その多角形の面積を平方度の単位で教えてくれます。平方メートルなどの平面単位で面積を取得するには、別の投影法 (どれ?) を使用してポリゴンの座標を変換する必要があると思います。

pyproj ライブラリがこれを行う方法を提供する必要があることを何度か読みました。pyprojを使用して、形の整ったPolygonオブジェクト全体を別の投影に変換してから面積を計算する方法はありますか?

私は自分のポリゴンを使って他のことを行います (あなたが今考えていることではありません)。特定の場合にのみ、面積を計算する必要があります。

これまでのところ、次の例しか見つかりませんでした: http://all-geo.org/volcan01010/2012/11/change-coordinates-with-pyproj/

これは、各 Polygon オブジェクトを外側のリングに分割し、存在する場合は内側のリングに分割し、座標を取得し、座標の各ペアを別の投影に変換して Polygon オブジェクトを再構築し、その面積を計算することを意味します (とにかく、単位は何ですか?)。これは解決策のように見えますが、あまり実用的ではありません。

より良いアイデアはありますか？

4

3 に答える 3

11

最終的に、matplotlib ライブラリの Basemap ツールキットを使用して作成しました。それがどのように機能するかを説明します。おそらく、これはいつか誰かに役立つでしょう。

1. システムに matplotlib ライブラリをダウンロードしてインストールします。 http://matplotlib.org/downloads.html

Windows バイナリについては、次のページを使用することをお勧めします: http://www.lfd.uci.edu/~gohlke/pythonlibs/#matplotlib 次のヒントに注意してください。

numpy、dateutil、pytz、pyparsing、six、およびオプションで pillow、pycairo、tornado、wxpython、pyside、pyqt、ghostscript、miktex、ffmpeg、mencoder、avconv、または imagemagick が必要です。

したがって、システムにまだインストールされていない場合は、matplotlib が適切に機能するように、numpy、dateutil、pytz、pyparsing、および six をダウンロードしてインストールする必要があります (Windows ユーザーの場合: Linux ユーザーの場合は、これらすべてがページにあります)。 、pythonパッケージマネージャー「pip」がうまくいくはずです）。

2. matplotlib の「basemap」ツールキットをダウンロードしてインストールします。http://matplotlib.org/basemap/users/installing.htmlから、 または Windows バイナリの場合はこちらから: http://www.lfd.uci.edu/~gohlke/pythonlibs/#basemap

3. Python コードでプロジェクションを実行します。

#Import necessary libraries
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np

#Coordinates that are to be transformed
long = -112.367
lat = 41.013

#Create a basemap for your projection. Which one you use is up to you.
#Some examples can be found at http://matplotlib.org/basemap/users/mapsetup.html
m = Basemap(projection='sinu',lon_0=0,resolution='c')

#Project the coordinates:
projected_coordinates = m(long,lat)

出力：

投影座標 (10587117.191355567、14567974.064658936)

そのような単純な。ここで、投影された座標を使用して shapely でポリゴンを構築し、その後 shapely の面積法で面積を計算すると、(使用した投影法に従って) 平方メートル単位で面積が得られます。平方キロメートルを求めるには、10^6 で割ります。

編集:単一の座標のみを変換するのではなく、ポリゴンのようなジオメトリ オブジェクト全体を変換するのに苦労しました。これは、多くのコードを書くことを意味しました

  • 多角形の外輪の座標を取得する
  • ポリゴンに穴があるかどうかを判断し、ある場合は各穴を個別に処理します
  • 外輪と任意の穴の座標の各ペアを変換します
  • 全体を元に戻し、投影された座標でポリゴン オブジェクトを作成します
  • それはポリゴン専用です...マルチポリゴンとジオメトリコレクションのために、これにさらにループを追加します

それから、物事をより簡単にする形の良いドキュメントのこの部分に出くわしました: http://toblerity.org/shapely/manual.html#shapely.ops.transform

たとえば、上記のように投影マップが設定されている場合:

m = Basemap(width=1,height=1, resolution='l',projection='laea',lat_ts=50,lat_0=50,lon_0=-107.)

次に、次の方法でこの射影を使用して、形の整ったジオメトリ オブジェクトを変換できます。

from shapely.ops import transform
projected_geometry = transform(m,geometry_object)
于 2014-05-19T09:12:49.830 に答える
9

測地線を計算します。これは非常に正確で、楕円体のみが必要です (投影ではありません)。これは、pyproj 2.3.0 以降で実行できます。

from pyproj import Geod
from shapely import wkt

# specify a named ellipsoid
geod = Geod(ellps="WGS84")

poly = wkt.loads('''\
POLYGON ((-116.904 43.371, -116.823 43.389, -116.895 43.407,
-116.908 43.375, -116.904 43.371))''')

area = abs(geod.geometry_area_perimeter(poly)[0])

print('# Geodesic area: {:.3f} m^2'.format(area))

# # Geodesic area: 13205034.647 m^2

abs()正の領域のみを返すために使用されます。ポリゴンの巻き方向によっては、負の領域を返す場合があります。

于 2020-10-02T01:15:44.040 に答える