I was writing a little utility to process a large number of polygons in a bunch of shape files, and wanted to have the area of each polygon. I thought I will do it using Python and OGR library. All my files are in WGS84 Latitude Longitude projection.
For experiment I tried the following code:
from osgeo import ogr
driver = ogr.GetDriverByName('ESRI Shapefile')
hinputfile = driver.Open(r'c:\GeoData\AZ.shp', 0)
in_layer = hinputfile.GetLayer(0)
feature = in_layer.GetFeature(0)
geom = feature.GetGeometryRef()
print(geom.GetArea())
This would give the area in square degrees, which I definitely not want. This would not be a good indicator, since the actual area would get smaller and smaller for the same value the bigger the latitude gets. I wanted to get the area in m² or km². This can be easily achieved by using a projection with units as meters. For example the «google» projection (EPSG:900913 - gOOglE get it?) which is now EPSG:3857 which is in meters. Here is the WKT for it.
PROJCS["WGS 84 / Pseudo-Mercator",
GEOGCS["Popular Visualisation CRS",
DATUM["Popular_Visualisation_Datum",
SPHEROID["Popular Visualisation Sphere",6378137,0,
AUTHORITY["EPSG","7059"]],
TOWGS84[0,0,0,0,0,0,0],
AUTHORITY["EPSG","6055"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.01745329251994328,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4055"]],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
PROJECTION["Mercator_1SP"],
PARAMETER["central_meridian",0],
PARAMETER["scale_factor",1],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
AUTHORITY["EPSG","3785"],
AXIS["X",EAST],
AXIS["Y",NORTH]]
So the code to calculate the area of each polygon in a shapefile in sq km would look something like the following:
from osgeo import ogr
from osgeo import osr
driver = ogr.GetDriverByName('ESRI Shapefile')
hinputfile = driver.Open(r'c:\GeoData\AZ.shp', 0)
in_layer = hinputfile.GetLayer(0)
# Prepare a transformation between projections
src_srs = in_layer.GetSpatialRef()
tgt_srs = osr.SpatialReference()
tgt_srs.ImportFromEPSG(3857)
transform = osr.CoordinateTransformation(src_srs, tgt_srs)
for feature in in_layer:
geom = feature.GetGeometryRef()
geom2 = geom.Clone()
geom2.Transform(transform)
area_in_sq_m = geom2.GetArea()
area_in_sq_km = area_in_sq_m / 1000000
print('Area in sq km: ', area_in_sq_km)
hinputfile.Destroy()
Hope this helps.