2

I have a PostGIS table holding cadastral (property) data stored in SRID 3857 (Web Mercator). In Python I want to calculate the length of polygon lines and the area in metres and sq metres respectively.

I'm extracting the geometry using geoalchemy into a shapely object. However, the length and area calculations are grossly incorrect.

EPSG:3857 says unit of measurement is metre.

Official records says the shape should be rectangle approx 15.24m x 45.72m with an area of 665.6sqm. When I load this feature in QGIS and verify using the ruler (old official records with their paper origins can be slightly out), I get similar measurements.

I have a feeling I need to transform from EPSG:3857 to euclidean plane but I have no idea how to do it.

Anyway my code below produces area of 1049.54sqm with ring lines of 55.73m, 18.95m, 55.70m & 18.72m.

Python:

# Get feature from database
property = db.session.query(UrbanProperty).filter(UrbanProperty.rid == item_id).first()
# Use geoalchemy2.shape.to_shape to get shapely object
property_shape = to_shape(property.shape)

# Print area
print("Area: {}".format(property_shape.area))

# Get ring lines
ring = property_shape.exterior.coords
last_point = None
lines = []
i = 0
for this_point in ring:
    if last_point is not None:
        i += 1
        line = LineString([last_point, this_point])
        print("Line {}: {}".format(i, line.length))
    last_point = this_point

Output:

Area: 1049.5458170682496 #should be ~700sqm
Line 1: 55.73504821619876 #should be ~46m
Line 2: 18.950509980859465 #should be ~15m
Line 3: 55.70698981912214 #should be ~46m
Line 4: 18.72151526076074 #should be ~15m

EPSG:3857 WKT Polygon:

POLYGON((16827043.551375 -4003534.72,
16827019.718375 -4003484.337625,
16827002.651875 -4003492.575125,
16827026.6795 -4003542.833875,
16827043.551375 -4003534.72))
4

The Web Mercator coordinate system distorts area and distance a lot. You're going to have a bad time trying to calculate these values in this CRS.

Your best bet is to convert to whatever your local UTM projection is and then calculate area. You may be able to do that using SQLAlchemy, but I don't know how, sorry. You could also get the geometry into Python, convert with Proj4, and then calculate areas on the new geometry...

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.