Extracting features from shapefiles

I tend to forget that ESRI shapefiles are more like databases than just bags of geometry. ogrinfo can extract features.

ogrinfo -geom=no -sql "select HUC_12,HU_10_NAME,HU_12_Name from WBD_Subwatershed" WBD_Subwatershed.shp

unfortunately the output is some terrible text format.

OGRFeature(WBD_Subwatershed):3
  HUC_12 (String) = 180703030303
  HU_10_NAME (String) = Lower San Luis Rey River
  HU_12_Name (String) = Pilgrim Creek

However it turns out using GDAL/OGR’s bindings to Python is relatively straightforward, even if the API is basically a straight-up rip of the C API. Decent tutorial is here.

import ogr

shapefile = ogr.Open("NHDPlusCA/NHDPlus18/WBDSnapshot/WBD/WBD_Subwatershed.shp")
layer = shapefile.GetLayerByName("WBD_Subwatershed")
for feature in layer:
    featureDefinition = layer.GetLayerDefn()
    for i in range(featureDefinition.GetFieldCount()):
        fieldDefinition = featureDefinition.GetFieldDefn(i)
        print fieldDefinition.GetName(),
        print "%s" % feature.GetFieldAsString(i)