1

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)
#layer.ResetReading()
field_name = ogr.FieldDefn("name", ogr.OFTString)
layer.CreateField(field_name)

### 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()
        myPoly.AddGeometry(cur_poly)
        cur_spat_ref = cur_poly.GetSpatialReference()

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

        featur_dfn = layer.GetLayerDefn()
        feature_abc = ogr.Feature(featur_dfn)
        feature_abc.SetGeometry(myPoly)
        feature_abc.SetField("name", cur_name)
        feature_abc.SetFID(FID)
        layer.CreateFeature(feature_abc)
        myPoly.Destroy()
        feature.Destroy()


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


datasource.Destroy()
  • 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. – Antonio Falciano Dec 17 '13 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 '13 at 9:11
  • There are a lot of file sharing websites today, then it's very easy to post a link in your question. – Antonio Falciano Dec 17 '13 at 9:16
  • @afalciano here are all the files, the script is called test_filter_shp.py – Claus Dec 17 '13 at 10:09
  • Zipping shared files is always welcome! ;) – Antonio Falciano Dec 17 '13 at 10:47
0

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()
        #myPoly.AddGeometry(cur_poly)
        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)
        feature_abc.SetFID(FID)
        layer.CreateFeature(feature_abc)
        #myPoly.Destroy()
        #feature_abc.Destroy()
  • thanks! It does work now. What is a 'null geometry' and when do I need to destroy features/polygons/datasources? – Claus Dec 17 '13 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. – Antonio Falciano Dec 17 '13 at 11:53

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.