2

I'm writing a small ogr/gdal script to get the area of census blocks:

in_file = data/zips/tabblock2010_01_pophu.zip
ds = ogr.Open(input_file)
lyr = ds.GetLayerByIndex(0)
lyr.ResetReading()
for feat in lyr:
  geom = feat.GetGeometryRef()
  print geom.GetArea()

This gives me numbers like 0.0003370560154999829 for block 010310106002015. Is this an area in meaningful units? If so, which ones? Alternatively, is it something that doesn't make sense like area in lat-long space, where I need to do some kind of projection before I can get the area in real distance units?

EDIT: If I call DumpReadable() I get:

POLYGON ((-86.121456 31.303452,-86.121207 31.304302,-86.120986 31.305055, ...

so it does look like this is in lat/long units, and so probably gives me area in lat-long space. How do I get real area?

  • 1
    Reproject data first into some projected coordinate system that is reasonable around -86.121456 31.303452. – user30184 Aug 28 '16 at 23:58
  • I would suggest WGS84 UTM Zone 16 North spatialreference.org/ref/epsg/wgs-84-utm-zone-16n (I assume your boundaries are in WGS84 Geographic). After geom = feat.GetGeometryRef() use geom.transformTo gdal.org/… to convert the geometry to the OGRSpatialRef there's a tutorial pcjericks.github.io/py-gdalogr-cookbook/projection.html that is fairly easy to follow (note EPSG is 32616 for the suggested spatial reference). – Michael Stimson Aug 29 '16 at 0:22
  • Is there a way to find a good projection for a small polygon in an automated way? (I want to find the area of every census tract in the US.) – Jeff Kaufman Aug 29 '16 at 0:52
  • There is a formula stackoverflow.com/questions/9186496/…. The zones start from -180 and are 6 degrees wide, discern the zone from this add 32600 to the zone to get the EPSG code. Note you can go outside the zone by a little bit (half a degree or so) before the data begins to distort noticeably, depending on your accuracy requirement and data coverage you might be able to get away with a single spatial reference like Lamberts for area calcs to within a few percent see gothos.info/2011/04/common-map-projection-definitions. – Michael Stimson Aug 29 '16 at 1:02
3

Bringing the comments together...

srcSR = osr.SpatialReference()
srcSR.ImportFromEPSG(4326) # WGS84 Geographic
destSR   = osr.SpatialReference()

lyr.ResetReading()
for feat in lyr:
    geom = feat.GetGeometryRef()
    if not geom.IsEmpty():                 # make sure the geometry isn't empty
        geom.AssignSpatialReference(srcSR) # you only need to do this if the shapefile isn't set or is set wrong
        env = geom.GetEnvelope()           # get the Xmin, Ymin, Xmax, Ymax bounds
        CentX = ( env[0] + env[2] ) / 2    # calculate the centre X of the whole geometry
        Zone  = int((CentX + 180)/6) + 1   # see http://gis.stackexchange.com/questions/13291/computing-utm-zone-from-lat-long-point
        EPSG  = 32600 + Zone               # get the EPSG code from the zone and the constant 32600 (all WGS84 UTM North start with 326)
        destSR.ImportFromEPSG(EPSG)        # create the 'to' spatial reference
        geom.TransformTo(destSR)           # project the geometry
        print geom.GetArea()               # get the area in square metres

The shapefile should have the spatial reference defined already, you can use QGIS to check or open the prj file in your favorite text editor (notepad would do) provided you can decipher the codes. If the shapefile does not have a spatial reference defined you will need to enforce the spatial reference with assignSpatialReference as shown in the code but be careful you're using the right one.. assumption can be dangerous. If your spatial reference is defined and is compatible with WGS84 (doesn't need a special transformation) you can comment out that line.

I've put in a check to ensure only valid geometries are calculated; sometimes geometries with no vertices sneak into the data and need to be checked for with the IsEmpty() function as projecting empty geometries may fail.

I'm calculating the middle X coordinate from the envelope, alternately you could try using the centroid of the geometry. You might also want to check the width (env[2] - env[0]) for very large features - you can't get that from a centroid!

  • Awesome, thanks! When I run this, though, I get "AttributeError: getEnvelope". If I rewrite the code to avoid getEnvelope by manually averaging the x-coordinates, then I get AttributeError: transformTo. – Jeff Kaufman Aug 29 '16 at 9:21
  • Um, it's because the methods don't match the docs (the C docs that is).. python is a case sensitive language so transformTo (from C) is not the same as TransformTo. I should have tested it first; I've modified the code changing GetEnvelope, AssignSpatialReference and TransformTo to their correct case (for python). Now I've tested it and it should work fine. – Michael Stimson Aug 29 '16 at 22:44

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.