I am looping over the features of a shapefile containing multiple polygons. If certain criteria are met, I want to save the current feature in a new shapefile.

I am doing this using GDAL 1.10.0 with Python 2.7.3

The code runs, creates the files I want, calculates the area correctly. In the shapefile that is created, the attribute table contains the propoer fields ('name' and FID). Two things are wrong:

  1. No polygons are visible (in ArcMap). This seems to me particularly weird, because the polygons seem to me proper, in the sense that they are closed (and the area is calculated properly).
  2. The coordinate system seems off (in Arc Catalog there appears in the bottom -20,9667 0,0000 (are these degrees??)), despite the fact that I am creating a projection file (they are in German national grid 3)

Any hints as to what the problems could be? I tried to put the relevant excerpt of the code below:

from osgeo import gdal, gdalnumeric, ogr, osr, gdalconst

### prepare shapefile into which to save
driver = ogr.GetDriverByName('ESRI Shapefile')
path_matched = r'.../.../path'
datasource = driver.CreateDataSource(path_matched)
lyr_matched = 'matched_%s' % cur_id  # from above
# Remove output shapefile if it already exists
name_shpfile = lyr_matched + '.shp'
if os.path.exists(os.path.join(path_matched, name_shpfile)):
    driver.DeleteDataSource(os.path.join(path_matched, name_shpfile))
layer = datasource.CreateLayer(lyr_matched, geom_type=ogr.wkbPolygon)
field_name = ogr.FieldDefn("name", ogr.OFTString)

### prepare file in which I want to write the areas of matched polygons
path_area = path_matched
file_area = 'matched_area_%s.txt' % (cur_id)
fobj_area = open(os.path.join(path_area, file_area), 'w')

### original shape file containing all polygons
shp = ogr.Open(shp_all)  # os fobj from above
lyr = shp.GetLayer()

# loop over polygons in original shameful
for FID in featList:
    feat = lyr.GetFeature(FID)

    # this is from the shapefile 
    cur_name = feat.GetField(1)

    # ... work on conditions ...

    if condition = true:
        myPoly = ogr.Geometry(type=ogr.wkbPolygon)
        cur_poly = feat.GetGeometryRef()
        cur_spat_ref = cur_poly.GetSpatialReference()

        cur_area =  cur_poly.GetArea()
        cur_area_string = "%s \t %f \n" % (gw_name_EZG, cur_area)

        featur_dfn = layer.GetLayerDefn()
        feature_abc = ogr.Feature(featur_dfn)
        feature_abc.SetField("name", cur_name)

# generate projection file
prj_filename = lyr_matched + '.prj'
file = open(os.path.join(path_matched, prj_filename), 'w')
targetSR = cur_spat_ref

  • It's quite tedious to debug your code without sample data. Although I see at least two bugs on the fly: cur_id (line 7) is not defined and if condition = true: (line 35) should be if condition == true: instead. Hope this helps. Dec 17, 2013 at 9:02
  • @afalciano thanks, but unfortunately, these are typos (cur_id see comment next to it). Is there a way to upload files somewhere so code can really be debugged?
    – Claus
    Dec 17, 2013 at 9:11
  • There are a lot of file sharing websites today, then it's very easy to post a link in your question. Dec 17, 2013 at 9:16
  • @afalciano here are all the files, the script is called test_filter_shp.py
    – Claus
    Dec 17, 2013 at 10:09
  • Zipping shared files is always welcome! ;) Dec 17, 2013 at 10:47

1 Answer 1


The issue happens because you're saving null geometries in the output shapefile. Here's the solution (commented rows are your wrong ones) in reference to your shared files:

if FID == 1:
        #myPoly = ogr.Geometry(type=ogr.wkbPolygon)
        cur_poly = feat.GetGeometryRef()
        cur_spat_ref = cur_poly.GetSpatialReference()

        cur_area =  cur_poly.GetArea()
        cur_area_string = "%s\n" % (cur_area)
        print cur_area_string

        featur_dfn = layer.GetLayerDefn()
        feature_abc = ogr.Feature(featur_dfn)
        feature_abc.SetGeometry(cur_poly) #feature_abc.SetGeometry(myPoly)
        feature_abc.SetField("name", cur_name)
  • thanks! It does work now. What is a 'null geometry' and when do I need to destroy features/polygons/datasources?
    – Claus
    Dec 17, 2013 at 11:40
  • A null (or empty) geometry is a not (or bad) initialized geometry. You usually need to destroy the datasources after you have done your job. It's a kind of "stop editing" features. Dec 17, 2013 at 11:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.