Calculating area of a lat,lon polygon in some useful units

Question

How do I tell the default geo spatial reference system? I have this code:

wkt = "POLYGON ((lon1 lat1, lon2 lat2, lon3 lat3))"
poly = ogr.CreateGeometryFromWkt(wkt)
print "Area = %d" % poly.GetArea()

I run it and get back a number. But what are the units? The documentation for OGRGeometry's get_Area() method says:

the area of the feature in square units of the spatial reference system in use.

So.. how do I tell what spatial reference system is currently in use?

Another way to formulate my question is:

I have a polygon defined as a sequence of lon,lat pairs. How do I get its area in sq kms?

So then I say: probably when people talk about areas of polygons on Earth they don't mean integrating the actual surface integral but area of some projection. So let me try setting the projection I want. And I got this code:

wkt = "POLYGON ((lon1 lat1, lon2 lat2, lon3 lat3))"
poly = ogr.CreateGeometryFromWkt(wkt)
sr = osr.SpatialReference()
sr.SetWellKnownGeogCS("NAD83")
sr.SetProjCS("SRS_PT_TRANSVERSE_MERCATOR")
poly.AssignSpatialReference(sr)
poly.GetArea()

but I get the same number as before! Printing out sr object gives me plausible looking output (it says "metre"!)

PROJCS["SRS_PT_TRANSVERSE_MERCATOR",
    GEOGCS["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6269"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4269"]],
    PROJECTION["unnamed"],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH]]

So I tried setting source and target systems and creating a transformation:

poly = ogr.CreateGeometryFromWkt(wkt)
sr = osr.SpatialReference()
sr_proj = osr.SpatialReference()
sr.SetWellKnownGeogCS("EPSG:4326"))
sr_proj.ImportFromEPSG(6991)
print(sr)
print(sr_proj)
poly.AssignSpatialReference(sr)
transform = osr.CoordinateTransformation(sr, sr_proj)
poly.Transform(transform)
print("Area = %f" % poly.GetArea())

gets me square meters! Yay. (but the problem now is... how do I pick which EPSG to use?? I don't know where the area is in advance!)

Answer 1

Nowadays, the most commonly used lat/long coordinate system is WGS 84. If you received lat/long coordinates of unknown/unmentioned coordinate system, you can assume that it was in WGS 84 (for example, GPS receiver work in WGS 84 by default). In fact, there are hundred of geographic (aka lat/long coordinate system, so using WGS84 could be inaccurate but you can't do anything about that). Fortunately, you focus on area and not on position. The shift in position could be something like 100 m, but all points are shifted and the impact on the area is much smaller. Note that this will distort the shape of your polygon, so make sure to densify the vertices of your polygon before projecting (if you have large polygons, or if you are close to the pole).

Conclusion : EPSG:4326 as input

For the choice of a projected coordinate system, local coordinate systems are designed to minimize the distortions, but if you don't work in a single region you need to check what to use based on the coordinates. For very specific coordinates systems this becomes quite complex. Alternatively, you could select the optimal UTM zone, North or South based on the latitude, then (ceil(longitude / 6)+31 to find the zone (then the EPSG code = 32600 + "zone_number" ). However, if you are interested in computing the area , you could take a global projection cylindrical equal area EPSG:9835 as a target.

Answer 2

I guess I'd recommend reprojecting your coordinates to some kind of planimetric coordinate system (such as a UTM system) before building your polygon. Note that different projections have different errors in computing areas. Also, when setting the projection, I'd either recommend using the projection's WKT if you have it or setting it from an EPSG code using the SpatialRefernce importFromEPSG function instead of trying to set it the way you're doing it. You can get a list of EPSG codes from spatialreference.org or just use Google and see what you find. For example, a commonly used code for UTM 31N with a WGS84 projection is here. The EPSG code for WGS84 (geographic lat/lon) is 4326. So you'd do something like this:

1. Reproject your coordinates using a "from" (source) and a "to" (target) OGRSpatialReference, and then use OGRCreateCoordinateTransformation to create a transformation and transform the needed points using the CoordinateTransformation.Transform function. I've provided links to the C++ docs but the python should be similar.
2. Do what you did in the original posting, but use the projected "planimetric" coordinates instead of the lat/lon coordinates. Also, I don't think it's essential to set the projection of your polygon to compute the area, but it's helpful if you're planning on exporting it at some point.

HTH