User Tools

Site Tools


qgis_basic_automation

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 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 that you will need to replace with your own (I will probably move these into static variables at some point).

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)
basic_map.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
qgis_basic_automation.txt · Last modified: 2021/02/19 20:31 by bushwalking