====== Basic Automation ======

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 [[nsw_lidar|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:
  * take a bounds file and clip a raster from a large virtual raster
  * create a depressionless DEM
  * create a slope raster and style it
  * create contours
    * remove short (<25m) contours
    * style contours
  * create and style streams (vector)
  * save project (if flag is set)

===== New version =====

This version has been tested in QGIS 3.34, with [[install_plugins|SAGA Next Gen and Whitebox Tools plugins]] installed, under Windows.

<code py basic_map.py>
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
</code>

===== Previous version =====

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.

<code py basic_map_old.py>
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
</code>