I'd like to generate contourlines from a SRTM image within Python. It seems to calculate but if I want to add my contourlines nothing shows up and the attribute table is empty as well. Please take a look at my code:
from osgeo import gdal, gdal_array
from osgeo.gdalconst import *
from numpy import *
from osgeo import ogr
#Read in SRTM data
indataset1 = gdal.Open( src_filename_1, GA_ReadOnly)
in1 = indataset1.GetRasterBand(1)
#Generate layer to save Contourlines in
ogr_ds = ogr.GetDriverByName("ESRI Shapefile").CreateDataSource(dst_filename)
contour_shp = ogr_ds.CreateLayer('contour')
field_defn = ogr.FieldDefn("ID", ogr.OFTInteger)
contour_shp.CreateField(field_defn)
field_defn = ogr.FieldDefn("elev", ogr.OFTReal)
contour_shp.CreateField(field_defn)
#Generate Contourlines
gdal.ContourGenerate(in1, 100, 0, [], 0, 0, contour_shp, 0, 1)
ogr_ds.Destroy()
Field ID and field elevation seem to be empty, but the contour_shape file is fairly huge ~100MB.
Any idea what might went wrong?
Update: I got it! I forgot to close the datasource with: ogr_ds.Destroy()
Don't use the Destroy()
method, as described in the GDAL/OGR Python Gotchas.
To save and close a dataset, dereference the variable, and optionally delete it.
I typically use this at the end to save/close either a GDAL or OGR dataset:
ogr_ds = None
del ogr_ds