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