The script below automates the production of a basic topographic map with contours, streams and slopes.
It builds on all of the steps on the Mapping from NSW Lidar page, so if you haven't worked through them, it will probably just confuse.
This uses PyQGIS, and can be run in the PyQGIS Console (Ctrl+Alt+P). Note that it can't be run as a standalone Python script, at least not without a lot of changes.
It requires some basic understanding of PyQGIS and coding, as there are a number of hard-coded file paths at the top of the script that you will need to replace with your own.
Key steps are:
This version has been tested in QGIS 3.34, with SAGA Next Gen and Whitebox Tools plugins installed, under Windows.
import os
import time
ts = str(int(time.time()))
user = {
'bounds':'C:/Users/brennant/Downloads/data.geojson',
'save_data':True,
'contour_sink_processing':True,
'fill_nodata':True,
'output_dir':'E:/geodata/blue_mountains/',
'output_file':'grose_valley'}
static = {
'input_dem':'E:/geodata/_data/56h.vrt',
'temp_folder':'C:/Users/brennant/AppData/Local/Temp/'
}
style = {
'slope':'E:/geodata/style/slope_cliffline_category.qml',
'contour':'E:/geodata/style/contours_5m_scale.qml',
'stream':'E:/geodata/style/basic_watercourse.qml'
}
params = {}
params['path_to_gpkg'] = user['output_dir'] + user['output_file'] + '.gpkg'
params['working_temp_folder'] = static['temp_folder'] + 'qgis_' + ts + '_' + user['output_file'] + '/'
os.mkdir(params['working_temp_folder'])
resultClip = processing.run("gdal:cliprasterbymasklayer", { 'ALPHA_BAND' : False, 'CROP_TO_CUTLINE' : True, 'DATA_TYPE' : 0, 'EXTRA' : '', 'INPUT' : static['input_dem'], 'KEEP_RESOLUTION' : False, 'MASK' : user['bounds'], 'MULTITHREADING' : False, 'NODATA' : None, 'OPTIONS' : '', 'OUTPUT' : 'TEMPORARY_OUTPUT', 'SET_RESOLUTION' : False, 'SOURCE_CRS' : None, 'TARGET_CRS' : None, 'X_RESOLUTION' : None, 'Y_RESOLUTION' : None })['OUTPUT']
resultFillDEM = processing.run("sagang:fillsinksxxlwangliu", { 'ELEV' : resultClip, 'FILLED' : 'TEMPORARY_OUTPUT', 'MINSLOPE' : 0.1 })['FILLED']
#### Slope
resultSlope1 = processing.run("native:slope",{'INPUT' : resultFillDEM, 'OUTPUT' : 'TEMPORARY_OUTPUT', 'Z_FACTOR' : 1 })['OUTPUT']
resultSlope = processing.run("gdal:rastercalculator",{ 'BAND_A' : 1, 'BAND_B' : None, 'BAND_C' : None, 'BAND_D' : None, 'BAND_E' : None, 'BAND_F' : None, 'EXTRA' : '', 'FORMULA' : '(A>40)*(A<=50)*1 + (A>50)*(A<=60)*2 + (A>60)*(A<=70)*3 + (A>70)*4', 'INPUT_A' : resultSlope1, 'INPUT_B' : None, 'INPUT_C' : None, 'INPUT_D' : None, 'INPUT_E' : None, 'INPUT_F' : None, 'NO_DATA' : None, 'OPTIONS' : '', 'OUTPUT' : 'TEMPORARY_OUTPUT', 'RTYPE' : 5 })['OUTPUT']
if user['save_data']:
processing.run("gdal:translate",{ 'COPY_SUBDATASETS' : False, 'DATA_TYPE' : 0, 'EXTRA' : '-co APPEND_SUBDATASET=YES -co RASTER_TABLE=slope', 'INPUT' : resultSlope, 'NODATA' : None, 'OPTIONS' : '', 'OUTPUT' : params['path_to_gpkg'], 'TARGET_CRS' : None })
rlSlope = iface.addRasterLayer("GPKG:" + params['path_to_gpkg'] + ":slope", "slope")
else:
rlSlope = iface.addRasterLayer(resultSlope)
rlSlope.setName('slope')
rlSlope.loadNamedStyle(style['slope'])
#### Contours
resultContour = processing.run("gdal:contour", { 'BAND' : 1, 'CREATE_3D' : False, 'EXTRA' : '', 'FIELD_NAME' : 'ELEV', 'IGNORE_NODATA' : False, 'INPUT' : resultFillDEM, 'INTERVAL' : 5, 'NODATA' : None, 'OFFSET' : 0, 'OUTPUT' : 'TEMPORARY_OUTPUT' })['OUTPUT']
resultContourSimplified = processing.run("native:simplifygeometries", { 'INPUT' : resultContour, 'METHOD' : 0, 'OUTPUT' : 'memory:buffer', 'TOLERANCE' : 1.0 })['OUTPUT']
vlContour = resultContourSimplified
pr = vlContour.dataProvider()
pr.addAttributes([QgsField("length", QVariant.Double)])
vlContour.updateFields()
expression1 = QgsExpression('$length')
context = QgsExpressionContext()
context.appendScopes(QgsExpressionContextUtils.globalProjectLayerScopes(vlContour))
with edit(vlContour):
for f in vlContour.getFeatures():
context.setFeature(f)
f['length'] = expression1.evaluate(context)
x = vlContour.updateFeature(f)
with edit(vlContour):
for f in vlContour.getFeatures():
if f['length'] < 25:
x = vlContour.deleteFeature(f.id())
if user['save_data']:
optionsContours = QgsVectorFileWriter.SaveVectorOptions()
optionsContours.layerName = 'contours'
optionsContours.actionOnExistingFile = QgsVectorFileWriter.CreateOrOverwriteLayer
writer = QgsVectorFileWriter.writeAsVectorFormat(vlContour, params['path_to_gpkg'], optionsContours)
vlContour = QgsVectorLayer(params['path_to_gpkg']+'|layername=contours', 'contours', 'ogr')
QgsProject.instance().addMapLayer(vlContour)
vlContour.setName('contours')
vlContour.loadNamedStyle(style['contour'])
#### Streams
resultFillDEM2 = processing.run("gdal:translate", { 'COPY_SUBDATASETS' : False, 'DATA_TYPE' : 0, 'EXTRA' : '', 'INPUT' : resultFillDEM, 'NODATA' : None, 'OPTIONS' : '', 'OUTPUT' : 'TEMPORARY_OUTPUT', 'TARGET_CRS' : None })['OUTPUT']
resultD8Flow = processing.run("wbt:D8FlowAccumulation", { 'clip' : False, 'esri_pntr' : False, 'input' : resultFillDEM2, 'log' : True, 'out_type' : 1, 'output' : params['working_temp_folder'] + 'd8flow' + ts +'.tif', 'pntr' : False })['output']
resultStreamNetwork = processing.run("sagang:channelnetwork", { 'CHNLNTWRK' : 'TEMPORARY_OUTPUT', 'CHNLROUTE' : 'TEMPORARY_OUTPUT', 'DIV_CELLS' : 10, 'DIV_GRID' : None, 'ELEVATION' : resultFillDEM2, 'INIT_GRID' : resultD8Flow, 'INIT_METHOD' : 2, 'INIT_VALUE' : 8.7, 'MINLEN' : 40, 'SHAPES' : 'TEMPORARY_OUTPUT', 'SINKROUTE' : None, 'TRACE_WEIGHT' : None })['SHAPES']
if user['save_data']:
processing.run("gdal:convertformat", { 'INPUT' : resultStreamNetwork, 'OPTIONS' : '-update -nln streams', 'OUTPUT' : params['path_to_gpkg'] })
vlStreams = QgsVectorLayer(params['path_to_gpkg']+'|layername=streams', 'streams', 'ogr')
else:
vlStreams = QgsVectorLayer(resultStreamNetwork, 'streams', 'ogr')
QgsProject.instance().addMapLayer(vlStreams)
vlStreams.loadNamedStyle(style['stream'])
#### Save project
if user['save_data']:
QgsProject.instance().write(user['output_dir'] + user['output_file'] +'.qgz')
# TODO: make vector layers read only so that they don't get accidentally selected in QField
This version worked in QGIS 3.16, but changes to SAGA since then mean it does not work in current (3.36) QGIS. It may be of use if you are still running very old versions.
params = {
'bounds':'C:/Users/brennant/Downloads/data.geojson',
'save_data':True,
'output_dir':'E:/geodata/wollemi/',
'output_file':'red_rocks'}
path_to_gpkg = params['output_dir'] + params['output_file'] + '.gpkg'
resultClip = processing.run("gdal:cliprasterbymasklayer", { 'ALPHA_BAND' : False, 'CROP_TO_CUTLINE' : True, 'DATA_TYPE' : 0, 'EXTRA' : '', 'INPUT' : 'E:/geodata/nsw100k.vrt', 'KEEP_RESOLUTION' : False, 'MASK' : params['bounds'], 'MULTITHREADING' : False, 'NODATA' : None, 'OPTIONS' : '', 'OUTPUT' : 'TEMPORARY_OUTPUT', 'SET_RESOLUTION' : False, 'SOURCE_CRS' : None, 'TARGET_CRS' : None, 'X_RESOLUTION' : None, 'Y_RESOLUTION' : None })
#resultFillDEM = processing.runAndLoadResults("saga:fillsinksxxlwangliu", { 'ELEV' : resultClip['OUTPUT'], 'FILLED' : 'TEMPORARY_OUTPUT', 'MINSLOPE' : 0.1 })
resultFillDEM = processing.run("saga:fillsinksxxlwangliu", { 'ELEV' : resultClip['OUTPUT'], 'FILLED' : 'TEMPORARY_OUTPUT', 'MINSLOPE' : 0.1 })
#### Slope
resultSlope = processing.run("native:slope",{'INPUT' : resultClip['OUTPUT'], 'OUTPUT' : 'TEMPORARY_OUTPUT', 'Z_FACTOR' : 1 })
if params['save_data']:
processing.run("gdal:translate",{ 'COPY_SUBDATASETS' : False, 'DATA_TYPE' : 0, 'EXTRA' : '-co APPEND_SUBDATASET=YES -co RASTER_TABLE=slope', 'INPUT' : resultSlope['OUTPUT'], 'NODATA' : None, 'OPTIONS' : '', 'OUTPUT' : path_to_gpkg, 'TARGET_CRS' : None })
rlSlope = iface.addRasterLayer("GPKG:" + path_to_gpkg + ":slope", "slope")
else:
rlSlope = iface.addRasterLayer(resultSlope['OUTPUT'])
rlSlope.setName('slope')
rlSlope.loadNamedStyle('E:/geodata/style/slope_cliffline.qml')
#### Contours
resultContour = processing.run("gdal:contour", { 'BAND' : 1, 'CREATE_3D' : False, 'EXTRA' : '', 'FIELD_NAME' : 'ELEV', 'IGNORE_NODATA' : False, 'INPUT' : resultFillDEM['FILLED'], 'INTERVAL' : 5, 'NODATA' : None, 'OFFSET' : 0, 'OUTPUT' : 'TEMPORARY_OUTPUT' })
resultContourSimplified = processing.run("native:simplifygeometries", { 'INPUT' : resultContour['OUTPUT'], 'METHOD' : 0, 'OUTPUT' : 'memory:buffer', 'TOLERANCE' : 1.0 })
vlContour = resultContourSimplified['OUTPUT']
pr = vlContour.dataProvider()
pr.addAttributes([QgsField("length", QVariant.Double)])
vlContour.updateFields()
expression1 = QgsExpression('$length')
context = QgsExpressionContext()
context.appendScopes(QgsExpressionContextUtils.globalProjectLayerScopes(vlContour))
with edit(vlContour):
for f in vlContour.getFeatures():
context.setFeature(f)
f['length'] = expression1.evaluate(context)
x = vlContour.updateFeature(f)
with edit(vlContour):
for f in vlContour.getFeatures():
if f['length'] < 25:
x = vlContour.deleteFeature(f.id())
if params['save_data']:
optionsContours = QgsVectorFileWriter.SaveVectorOptions()
optionsContours.layerName = 'contours'
optionsContours.actionOnExistingFile = QgsVectorFileWriter.CreateOrOverwriteLayer
writer = QgsVectorFileWriter.writeAsVectorFormat(vlContour, params['output_dir'] + params['output_file'] + '.gpkg', optionsContours)
vlContour = QgsVectorLayer(path_to_gpkg+'|layername=contours', 'contours', 'ogr')
QgsProject.instance().addMapLayer(vlContour)
vlContour.setName('contours')
vlContour.loadNamedStyle('E:/geodata/style/5m_contours.qml')
#### Streams
resultCatchment = processing.run("saga:catchmentarea", { 'ELEVATION' : resultFillDEM['FILLED'], 'FLOW' : 'TEMPORARY_OUTPUT', 'METHOD' : 5 })
resultStreamNetwork = processing.run("saga:channelnetwork", { 'CHNLNTWRK' : 'TEMPORARY_OUTPUT', 'CHNLROUTE' : 'TEMPORARY_OUTPUT', 'DIV_CELLS' : 10, 'DIV_GRID' : None, 'ELEVATION' : resultFillDEM['FILLED'], 'INIT_GRID' : resultCatchment['FLOW'], 'INIT_METHOD' : 2, 'INIT_VALUE' : 5000, 'MINLEN' : 40, 'SHAPES' : 'TEMPORARY_OUTPUT', 'SINKROUTE' : None, 'TRACE_WEIGHT' : None })
if params['save_data']:
processing.run("gdal:convertformat", { 'INPUT' : resultStreamNetwork['SHAPES'], 'OPTIONS' : '-update -nln streams', 'OUTPUT' : path_to_gpkg })
vlStreams = QgsVectorLayer(path_to_gpkg+'|layername=streams', 'streams', 'ogr')
else:
vlStreams = QgsVectorLayer(resultStreamNetwork['SHAPES'], 'streams', 'ogr')
QgsProject.instance().addMapLayer(vlStreams)
vlStreams.loadNamedStyle('E:/geodata/style/basic_watercourse.qml')
#### Save project
if params['save_data']:
QgsProject.instance().write(params['output_dir'] + params['output_file'] +'.qgz')
# TODO: make vector layers read only so that they don't get accidentally selected in QField
# TODO: save depressionless DEM as separate GPKG in case needed? (or should this just be regenerated if required?)
# TODO: set up static parameters