У меня есть GeoJSON
{
"type": "FeatureCollection",
"features": [
{
"type": "Feature",
"properties": {},
"geometry": {
"type": "Polygon",
"coordinates": [
[[13.65374516425911, 52.38533382814119], [13.65239769133293, 52.38675829106993], [13.64970274383571, 52.38675829106993], [13.64835527090953, 52.38533382814119], [13.64970274383571, 52.38390931824483], [13.65239769133293, 52.38390931824483], [13.65374516425911, 52.38533382814119]]
]
}
}
]
}
который http://geojson.io отображается как
Я хотел бы рассчитать его площадь (87106,33 м ^ 2) с помощью Python. Как я могу это сделать?
Что я пробовал
# core modules
from functools import partial
# 3rd pary modules
from shapely.geometry import Polygon
from shapely.ops import transform
import pyproj
l = [[13.65374516425911, 52.38533382814119, 0.0], [13.65239769133293, 52.38675829106993, 0.0], [13.64970274383571, 52.38675829106993, 0.0], [13.64835527090953, 52.38533382814119, 0.0], [13.64970274383571, 52.38390931824483, 0.0], [13.65239769133293, 52.38390931824483, 0.0], [13.65374516425911, 52.38533382814119, 0.0]]
polygon = Polygon(l)
print(polygon.area)
proj = partial(pyproj.transform, pyproj.Proj(init='epsg:4326'),
pyproj.Proj(init='epsg:3857'))
print(transform(proj, polygon).area)
Это дает 1.1516745933889345e-05
и 233827.03300877335
- что первое не имеет никакого смысла, как ожидалось, но как мне исправить второе? (Я понятия не имею, как установить параметр инициализации pyproj.Proj
)
Думаю, epsg:4326
имеет смысл, поскольку это WGS84 (источник), но для epsg:3857 Я не уверен.
Лучшие результаты
Следующее намного ближе:
# core modules
from functools import partial
# 3rd pary modules
import pyproj
from shapely.geometry import Polygon
import shapely.ops as ops
l = [[13.65374516425911, 52.38533382814119, 0],
[13.65239769133293, 52.38675829106993, 0],
[13.64970274383571, 52.38675829106993, 0],
[13.64835527090953, 52.38533382814119, 0],
[13.64970274383571, 52.38390931824483, 0],
[13.65239769133293, 52.38390931824483, 0],
[13.65374516425911, 52.38533382814119, 0]]
polygon = Polygon(l)
print(polygon.area)
geom_area = ops.transform(
partial(
pyproj.transform,
pyproj.Proj(init='EPSG:4326'),
pyproj.Proj(
proj='aea',
lat1=polygon.bounds[1],
lat2=polygon.bounds[3])),
polygon)
print(geom_area.area)
это дает 87254,7 м ^ 2 - это все еще 148 м ^ 2 отличается от того, что говорит geojson.io. Почему это так?