I'm using a few reference StackExchange/Overflow posts here - and manual, but not coming up with what I'm after.
I have a GeoTIFF that I'm importing, and I need to develop out contours at regular intervals.
I've repurposed code from this post, but the two issues I'm having are:
- I can't seem to specify regular intervals, rather I have to load my intervals manually into an array
- The shapefile is created without a prj file, meaning it's not spatially referenced.
My code so far:
from osgeo import gdal from osgeo.gdalconst import * from numpy import * from osgeo import ogr #make contours from the raster print('Making Contours') indataset1 = gdal.Open(os.sep.join([output_dir, os.path.basename(sabre_out_dir)]) + 'mdROS.tif', GA_ReadOnly) in1 = indataset1.GetRasterBand(1) dst_filename = os.sep.join([output_dir, os.path.basename(sabre_out_dir)]) + 'mdROS_contour.shp' 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("mdROS", ogr.OFTReal) contour_shp.CreateField(field_defn) #Generate Contourlines - Instead of manually specifying all intervals I'd like to say "every X interval" gdal.ContourGenerate(in1, 0, 0, [1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000,9000,10000], 0, 0, contour_shp, 0, 1) ogr_ds.Destroy()