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