Вычислить площадь многоугольника в плоских единицах (например, квадратных метрах) в Shapely

Я использую Python 3.4 и shapely 1.3.2 для создания объекта Polygon из списка пар координат long / lat, который я преобразовываю в известную текстовую строку, чтобы проанализировать их. Такой многоугольник может выглядеть так:

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, а затем вычисление его площади (в любом случае, что это за единица?). Это похоже на решение, но не очень практично.

Есть идеи получше?


person Dirk    schedule 16.05.2014    source источник


Ответы (3)


Хорошо, наконец-то я сделал это с помощью набора инструментов Basemap библиотеки matplotlib. Я объясню, как это работает, может быть, это кому-нибудь когда-нибудь пригодится.

1. Загрузите и установите библиотеку matplotlib в вашу систему. http://matplotlib.org/downloads.html

Для двоичных файлов Windows я рекомендую использовать эту страницу: http://www.lfd.uci.edu/~gohlke/pythonlibs/#matplotlib Остерегайтесь подсказки, в которой говорится:

Требуется numpy, dateutil, pytz, pyparsing, six и, возможно, подушка, pycairo, tornado, wxpython, pyside, pyqt, ghostscript, miktex, ffmpeg, mencoder, avconv или imagemagick.

Следовательно, если он еще не установлен в вашей системе, вам необходимо загрузить и установить numpy, dateutil, pytz, pyparsing и шесть, чтобы matplotlib работал правильно (для пользователей Windows: все они могут быть найдены на странице, для пользователей Linux , диспетчер пакетов python "pip" должен помочь).

2. Загрузите и установите набор инструментов "basemap" для matplotlib. Либо из 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)

Выход:

projected_coordinates (10587117.191355567, 14567974.064658936)

Просто как тот. Теперь, когда вы используете спроецированные координаты для построения многоугольника с shapely, а затем вычисляете площадь с помощью метода shapely area, вы получаете площадь в квадратных метрах (в соответствии с используемой вами проекцией). Чтобы получить квадратные километры, разделите на 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)
person Dirk    schedule 19.05.2014
comment
См. Также пример преобразования проекции здесь. Обратите внимание, что точность расчета площади зависит от проекции. Например, используйте зону UTM, если все данные хорошо вписываются в одну. - person Mike T; 09.06.2014

Вычислите геодезическую область, которая очень точна и требует только эллипсоида (не проекции). Это можно сделать с помощью 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() используется для возврата только положительных областей. Отрицательная область может быть возвращена в зависимости от направления намотки многоугольника.

person Mike T    schedule 02.10.2020

Перевести в радианы и предположить, что Земля представляет собой идеальную сферу радиусом 6370 км:

p = np.array ([[- 116.904,43.371], [- 116.823, 43.389], [- 116.895,43.407], [- 116.908,43.375], [- 116.904,43.371]])

poly = Многоугольник (np.radians (p))

poly.area = 4.468737548271707e-07

poly.area * 6370 ** 2 = 18,132751662246623

person DST    schedule 12.06.2014