2

I am having difficulties after loading a shapefile with a sql statement with ogr. The layer and its features seem to be correctly loaded. Still I am not able to write them into a new shapefile or pass the layer to the gdal_rasterize function. The shapefile remains an empty file and the raster has the typical no input/something wrong happened values from 1.79..e+308 to -1.79..e+308 (not visualized in QGIS). Code example for the creation of the shapefile (paths are omitted):

file = ogr.Open(file_path)
sql = 'Select * From {} where {} <= 5'.format(layer_name, "CS")
tmp_layer = file.ExecuteSQL(sql)
print(tmp_layer)
>>><osgeo.ogr.Layer; proxy of <Swig Object of type 'OGRLayerShadow *' at 0x00000214D99617E0> >
for feat in tmp_layer:
    print(feat.GetField('CS'))
    >>>(prints the field of the subset of features as expected)
drv = ogr.GetDriverByName("ESRI Shapefile")
out_ds = drv.CreateDataSource(output_folder + '/asd.shp')
proj = tmp_layer.GetSpatialRef()
layer = out_ds.CreateLayer('asd', proj, ogr.wkbMultiPolygon)
# copy the schema of the original shapefile to the destination shapefile
lyr_def = tmp_layer.GetLayerDefn()
for i in range(lyr_def.GetFieldCount()):
    layer.CreateField(lyr_def.GetFieldDefn(i))
#write features from tmp_layer to new layer
for feat in tmp_layer:
    layer.CreateFeature(feat)
out_ds = layer = None

What is my mistake?

  • Does it work without the filtering? Does it work with different data (try naturalearthdata.com since that is free for everyone to use)? – bugmenot123 Jun 18 at 15:38
  • 1
    Oh damn... I tried so much and I dont know what wrent wrong before - because I am quite sure I tested it without the iteration over the features. It seems that while iterating over the features these get consumed and are not available anymore... And thats the whole problem! (So it didnt work without the sql statement) – gilla Jun 18 at 16:16
2

When you iterate over the layer it essentially calls lyr.GetNextFeature(). After the last feature the pointer remains at the end which is why when you iterate over the layer the second time it doesn't find any features.

If you plan on iterating over the layer more than once you need to make sure you call lyr.ResetReading() to reset the pointer to the first feature.

#write features from tmp_layer to new layer
tmp_layer.ResetReading()
for feat in tmp_layer:
    layer.CreateFeature(feat)
1

The iteration over the features seem to consume these features like in an interator. Thus omitting the iteration for printing will solve the problem.

  • You can accept your answer so future readers know what the solution is. – Marcelo Villa Jun 18 at 19:21

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.