How can I create contours from GeoTIFF and Python (GDAL, Rasterio, etc?) into shapefile with prj file

Question

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

Answer 1

You are not defining a spatial reference for the shapefile in your code. In order to do so, you have to specify it as the second argument when using the CreateLayer() method. For example, if your spatial reference would happen to be WGS84, you could write the following:

from osgeo import osr

sr = osr.SpatialReference()
sr.ImportFromEPSG(4326)
contour_shp = ogr_ds.CreateLayer('contour', sr)

I guess you want the spatial reference to be the same as the raster, so you could do the following as well:

from osgeo import osr

sr = osr.SpatialReference(indataset1.GetProjection())
contour_shp = ogr_ds.CreateLayer('contour', sr)

---

As for the intervals, I am not familiar with the ContourGenerate() function. However, you can easily create a list with intervals using Python's range() function by specifying a step as the third argument. For example:

intervals = list(range(1000, 11000, 1000))  # start, end, step

print(intervals)
# [1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000]

Then, you can simply specify the intervals you just created as the argument.

gdal.ContourGenerate(in1, 0, 0, intervals, 0, 0, contour_shp, 0, 1)

Answer 2

This is what I wound up going with.

The answer shows how to produce contours from a shapefile filled with polygons in a grid formation (think fishnet). The same approach should work for points, perhaps with a minimal amount of change.

Input can be shapefile or geopackage. Might work with other formats.

import os
import geopandas as gpd
from geocube.api.core import make_geocube
from osgeo import gdal
from osgeo.gdalconst import *
from osgeo import osr
from numpy import *
from osgeo import ogr

#mdROS is the field I'm wanting to work with from the original shp,
#so I have followed it all the way through the code. You could variablise this. 

# Make a raster using the shapefile - could go straight from dataframe
    print('Making Raster')
    input_geopackage = os.sep.join([output_dir, os.path.basename(out_dir)]) + 'input_file.shp'
    out_grid = make_geocube(
        vector_data=input_geopackage,
        resolution=(-30, 30 ),
    )
    out_grid["mdROS"].rio.to_raster(os.sep.join([output_dir, os.path.basename(out_dir)]) + 'mdROS.tif')

    #make contours from the raster
    print('Making Contours')
    indataset1 = gdal.Open(os.sep.join([output_dir, os.path.basename(out_dir)]) + 'mdROS.tif', GA_ReadOnly)
    sr = osr.SpatialReference(indataset1.GetProjection())
    in1 = indataset1.GetRasterBand(1)
    dst_filename = os.sep.join([output_dir, os.path.basename(out_dir)]) + 'mdROS_contour.shp'
    ogr_ds = ogr.GetDriverByName("ESRI Shapefile").CreateDataSource(dst_filename)
    #contour_shp = ogr_ds.CreateLayer('contour') The line below outputs a spatially referenced shp
    contour_shp = ogr_ds.CreateLayer('contour', sr)

    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 - ContourGenerate from https://gdal.org/python/osgeo.gdal-module.html#ContourGenerate
    gdal.ContourGenerate(
        in1,    #Band srcBand
        500,      #double contourInterval - This defines contour intervals
        0,      #double contourBase
        [],      #int fixedLevelCount
        0,      #int useNoData
        0,      #double noDataValue
        contour_shp, #Layer dstLayer
        0,      #int idField
        1       #int elevField
        )
    #ogr_ds.Destroy() - Better not to use this method. See here: https://gdal.org/api/python_gotchas.html
    #instead, empty and delete the dataframe
    ogr_ds = None
    del ogr_ds