0

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!)

  • There is no default. You need to know the spatial reference system of your data. Note that trying to get area from a polygon with coordinates in lon, lat will will give you a pretty useless value, the units are decimal degrees. – user2856 Sep 18 at 3:05
  • what are the coordinates usually specified as? I thought using lon, lat was pretty standard for this? the number, btw, looks sort of meaningful -- it appears to be 1/10,000 of sq km. But i am not sure why. Perhaps I should rephrase my question. I have a polygon in lon,lat pairs. How do I get its area in sq kms? – MK. Sep 18 at 3:11
  • I guess I don't understand some fundamental concepts here. Lon, Lat pairs are coordinates on a sphere. WKT format seems to work with these. Is getArea() defined on a sphere or do I need to convert to plane projection coordinates first? – MK. Sep 18 at 3:35
  • okay I understand the 10,000 thing now -- I am giving inputs in degrees, and get output in degrees, and degree is roughly equal to 100km (depending on where you are of course). So the question changes to: what is the proper way of calculating the area based on lat, lon polygon... – MK. Sep 18 at 3:53
  • 2
    Tell first that the lat/lon polygon is in this system epsg.io/4326. Then convert geometry into some projected coordinate system that is reasonable on the area of your data. Finally get the area of the projected polygon. – user30184 Sep 18 at 4:56
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.

| improve this answer | |
  • Is EPSG:9835 preserves areas which it sounds like it does, that's what I need. But when I try sr_proj.ImportFromEPSG(9835) I get "ERROR 1: PROJ: proj_create_from_database: crs not found" – MK. Sep 18 at 14:28
0

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

| improve this answer | |
  • I am not sure I understand how your suggestion is different from the last version of the code I have. I need a source and target EPSGs, right? I think you are saying project individual coordinates instead of the whole poly? But transform from what to what? – MK. Sep 18 at 6:02
  • In the last version of the code, you're still creating the polygon from the geographic (lat/lon) coordinates. You need to create the polygon from some projected coordinates to compute a valid area. You can use the projection mentioned by @radouxju or some whatever projection if you wish, but you have to reproject the coordinates before calling ogr.CreateGeometryFromWkt. – Burrow Sep 18 at 15:58
  • I thought my code says "takes these coordinates, interpret them as some reasonable geo coord system, transform them to a projected system and calculate area. I think you are saying "apply transform before creating geo" but I'm not sure what the difference? – MK. Sep 18 at 18:26

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.