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