6

I'm working with GeoDjango and have a polygon in the WGS84 projection. I need to find its area in square meters. When I call polygon.area, I get some float out and I have no idea what the units of it are.

For example, Central Park in New York (843 acres according to Wikipedia) has roughly these coordinates (WGS84):

MULTIPOLYGON (((-73.9582363278800017 40.8002491235089977, -73.9489666135250019 40.7967404799830007, -73.9729992062979989 40.7639841980140005, -73.9820972592760029 40.7681446212439980, -73.9582363278800017 40.8002491235089977)))

Accessing the "area" property by calling polygon.area yields 0.00038966160697

When I plug the polygon WKT into PostGIS directly and call ST_Area (as suggested here: How can I measure area from geographic coordinates?), I get 3650892.09953911 (which is in square meters), which comes out to about 920 acres, which is about right (I didn't digitize Central Park perfectly). I'd prefer not to have to do raw SQL calls in the code to avoid potential SQL injections. Is there a way to get to this number straight from GEOS?

Improve this question
3
  • For the record: the correct area of this polygon is about 3,652,600 square meters which is approximately 902.6 acres. As a check, multiplying the result of polyon.area (which is in squared degrees) by the cosine of the mid-latitude, roughly cos(40.78), gives 900.1 acres.
    – whuber
    Commented May 14, 2011 at 18:24
  • @whuber - how are you determining the "correct" area of the polygon? Are you using the local projection for NY? I'm trying to be able to do this across the world - is there generally a good way to do this? I'm surprised that the PostGIS ST_Area function isn't more accurate...
    – Max
    Commented May 16, 2011 at 0:41
  • I'm using the method I recommended in a comment to @scw's answer; namely, I projected with an accurate equal-area projection (the standard Albers for the conterminous US). Allow me to repeat myself: for global datasets, project using a global equal-area projection and compute the projected area. You can use something as simple (and crude) as a cylindrical equal-area projection, which is fine away from the poles, and for polar work use a polar aspect of any local equal-area projection.
    – whuber
    Commented May 16, 2011 at 13:04

1 Answer 1

Reset to default
6

A couple of ways you can solve this: you could switch out your geometry field for a geography field within GeoDjango, then any native area calls should return values in meters which can easily be converted to the units you're interested in.

If you want to stick to storing things in geographic space as plain geometries, then you'll need to do conversion somewhere: you can't calculate area from geographic coordinates. You could instead add a helper method which is in a projected coordinate space to do the transformation:

def get_acres(self): 
    """ 
    Returns the area in acres. 
    """ 
    # Convert our geographic polygons (in WGS84)
    # into a local projection for New York (here EPSG:32118)
    self.polygon.transform(32118)
    meters_sq = self.polygon.area
    
    acres = meters_sq * 0.000247105381  # meters^2 to acres

    return acres

Which projection you should use depends on the extent of the data and how accurate you need the results: here I've illustrated with a specific projection for part of New York, but if your data isn't particularly accurate, you could easily substitute a global projection or just use a simple formula as @whuber mentioned in his comment.

Improve this answer
5
  • +1 In principle any equal-area projection appropriate for your datum will give the same (correct) area, up to the inherent accuracy of its calculation.
    – whuber
    Commented May 14, 2011 at 21:47
  • 1
    and @whuber - thanks for your quick replies! Like scw suggested, I tried the geography field approach (although I would have preferred to avoid creating additional models in the database) like this: class DummyBounds(models.Model): bounds = models.MultiPolygonField(geography=True) objects = models.GeoManager() In my code I write: db = DummyBounds(bounds=polygon) sq_meter_area = db.bounds.area print sq_meter_area And it appears that my "sq_meter_area" variable is actually still in square degrees (I get 0.000386078620723). Am I doing something wrong?
    – Max
    Commented May 15, 2011 at 2:53
  • @Max, I've updated my example -- could you try calling area.sq_m to see if that produces correct results?
    – scw
    Commented May 15, 2011 at 5:39
  • 6
    The area property of a GEOS object is a float (not an Area object), or at least that's what I'm seeing (so you can't call .sq_m).
    – Max
    Commented May 16, 2011 at 0:40
  • I have a geography field but still needed the transform for the area to be correct.
    – Akaisteph7
    Commented Apr 4 at 14:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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