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:
- 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).
- 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()