#!/usr/bin/env python3 # -*- coding: UTF-8 -*- # from osgeo import gdal, gdalconst # from osgeo_utils import gdal_polygonize # import rasterio # import numpy as np # from pcraster import setclone,lddcreate,catchment,readmap,__version__,streamorder,spatial,threading,report # from pcraster.multicore._operations import set_nr_worker_threads,nr_worker_threads reclassif_flux = { 5:{'min':1,'max':1}, 4:{'min':2,'max':2}, 3:{'min':3,'max':3}, 6:{'min':4,'max':4}, 255:{'min':5,'max':5}, 2:{'min':6,'max':6}, 7:{'min':7,'max':7}, 0:{'min':8,'max':8}, 1:{'min':9,'max':9}, } def createldd( src_ds:str,ldd_file:str,outflowdepth:int=9999999,corevolume:int=9999999,corearea:int=9999999, catchmentprecipitation:int=9999999,lddin=False,unitcell=False): """ Transformation du GeoTIFF au format PCRaster par l'utilisation de la fonction gdal.Translate . Calcul de la direction des flux. Création d'un fichier .map représentant la direction des flux. Parameters ---------- src_ds : Elevation (MNT/DEM) au format GeoTIFF. ldd_file : Fichier de sortie de direction des flux au format ``.map`` . outflowdepth : spatial, non spatial scalar corevolume : spatial, non spatial scalar corearea : spatial, non spatial scalar catchmentprecipitation : spatial, non spatial scalar lddin : spatial, non spatial scalar unitcell : spatial, non spatial scalar See Also -------- https://pcraster.geo.uu.nl/pcraster/4.4.0/documentation/pcraster_manual/sphinx/op_lddcreate.html#index-0 """ from pcraster import lddcreate,report,setglobaloption,setclone from osgeo import gdalconst,gdal from osgeo_utils import gdal_calc from datetime import datetime as dt from os import sched_getaffinity,listdir,path from pcraster.multicore._operations import set_nr_worker_threads,nr_worker_threads def _lddcreate(dst_file,outflowdepth,corevolume,corearea,catchmentprecipitation): yield lddcreate(dst_file,outflowdepth,corevolume,corearea,catchmentprecipitation) cpu_dipo = len(sched_getaffinity(0))-1 # if nr_worker_threads() != cpu_dipo: # set_nr_worker_threads(cpu_dipo) timeinit = dt.now() print(timeinit) OUTPATH = ldd_file.rsplit('/',1)[0] pcr_temp = 'flow_direction_'+ldd_file.rsplit('/',1)[1] dst_file = path.join(OUTPATH,pcr_temp) if pcr_temp in listdir(OUTPATH): print(" PCRaster format already exist.") else: # Convertir GeoTIFF au format PCRaster print(" INIT gdal.Translate to format PCRaster.") src_ds = gdal.Open(raster50) gdal.Translate(dst_file, src_ds, format='PCRaster', outputType=gdalconst.GDT_Float32, metadataOptions="VS_SCALAR",callback=gdal.TermProgress) print(" END gdal.Translate to format PCRaster :",dt.now()-timeinit) # gdal_calc.py -A input_raster.tif --calc="(A==max(A.flatten()))*1" --outfile=max_pixel.tif --overwrite if ldd_file.rsplit('.',1)[1] != 'map': ldd_file += '.map' # Calcule de la direction du flux setglobaloption("lddin") if lddin else setglobaloption("lddout") setglobaloption("unitcell") if unitcell else setglobaloption("unittrue") setclone(dst_file) print(" INIT create local drain direction map with flow directions.") ldd = _lddcreate(dst_file, outflowdepth, corevolume, corearea, catchmentprecipitation) # Save result try : report(next(ldd),ldd_file) print(" END create local drain direction map with flow directions :",dt.now()-timeinit) print(" END :",dt.now()-timeinit) return ldd except Exception as e: print(" END :",dt.now()-timeinit) print(e) return ldd if __name__ == "__main__": # if nr_worker_threads != 7: # set_nr_worker_threads(7) PATH = "/media/colas/SRV/FICHIERS/SITES/DISTRICTS NATURELS/BASSE VALLEE DE L'ISERE/SONE/TUFI_Sone-à-Soi/SIG/" raster = PATH+"MNT_1m.tif" dst_filename = '/home/colas/Documents/tmp/OUTPT.map' xyz_filename = '/home/colas/Documents/tmp/OUTPT.xyz' ldd_filename = '/home/colas/Documents/tmp/ldd.map' strahler_file = '/home/colas/Documents/tmp/strahler.map' localPATH = "/home/colas/Documents/tmp/LOUIS/sone_tuff/" raster50 = PATH+"MNT_50cm.tif" dst_file50 = localPATH+'ldd_50cm.map' createldd( src_ds = raster50, ldd_file = dst_file50 ) # # src = rasterio.open(raster) # src_ds = gdal.Open(raster50) # srcband = src_ds.GetRasterBand(1) # # srcdata = srcband.ReadAsArray() # stats = srcband.GetStatistics(True, True) # min_alti = srcband.ComputeRasterMinMax()[0] # print("[ STATS ] = Minimum=%i, Maximum=%i, Mean=%i, StdDev=%i"%(stats[0], stats[1], stats[2], stats[3])) # # Calcule du pixel ayant la plus basse altitude # # gdal_calc.Calc(A=dst_file,calc="A==min(A.flatten())",outfile=OUTPATH+'/min_pixel.tif',overwrite=True) # gdal_calc.Calc(A=dst_file,calc="(A==numpy.min(A.flatten()))*1",outfile=OUTPATH+'/min_pixel.tif',overwrite=True) # # Convertir GeoTIFF au format PCRaster # gdal.Translate(dst_file50, src_ds, format='PCRaster', outputType=gdalconst.GDT_Float32, # metadataOptions="VS_SCALAR") # # gdal.Translate(xyz_filename, src_ds, format='XYZ') # # gdal.Translate(xyz_filename, src_ds, format='XYZ') # # gdal.Info(src_ds) # # gdal.Info(raster) # setclone(dst_filename) # Dem = readmap(dst_filename) # # Calculer la direction du flux # ldd = lddcreate(Dem, 9999999,9999999,9999999,9999999) # ldd = lddcreate(dst_filename, 1e31, 1e31, 1e31, 1e31) # 4h # ldd = readmap(ldd_filename) # # Détermination des flux / Calculer les commandes Strahler # # Plus l'ordre est élevé, plus le flux est important. # strahler = streamorder(ldd) # report(strahler,strahler_file) # save result # # Calculer le réseau de canaux # from pcraster import spatial,boolean,ordinal # input_nonspatial = 8 # setclone(strahler_file) # SpatialResultB = spatial(boolean(input_nonspatial)) # SpatialResultO = spatial(ordinal(input_nonspatial)) # strahler = readmap(strahler_file) # ResultComparison = strahler >= SpatialResultO # # Définir le point de sortie # points = '' # catchm = catchment(ldd, points) # # processing.run("native:reclassifybytable", {'INPUT_RASTER':'/home/colas/Documents/tmp/ldd.map','RASTER_BAND':1,'TABLE':['1','1','5','2','2','4','3','3','3','4','4','6','5','5','255','6','6','2','7','7','7','8','8','0','9','9','1'],'NO_DATA':-9999,'RANGE_BOUNDARIES':2,'NODATA_FOR_MISSING':False,'DATA_TYPE':5,'OUTPUT':'TEMPORARY_OUTPUT'}) # # processing.run("pcraster:spatial", {'INPUT':8,'INPUT1':0,'INPUT2':'/home/colas/Documents/tmp/ldd.map','OUTPUT':'TEMPORARY_OUTPUT'}) # # qgis_process run pcraster:spatial --distance_units=meters --area_units=m2 --ellipsoid=EPSG:7019 --INPUT=8 --INPUT1=0 --INPUT2=/home/colas/Documents/tmp/ldd.map --OUTPUT=TEMPORARY_OUTPUT # # qgis_process run pcraster:comparisonoperators --distance_units=meters --area_units=m2 --ellipsoid=EPSG:7019 --INPUT='%3Fcrs%3DEPSG%3A2154%26extent%3D868999.5%2C6447000.5%2C885999.5%2C6467000.5%26width%3D17000%26height%3D20000%26formula%3D%2522strahler%401%2522%2520%253E%253D5%2520%26strahler%3Auri%3D%2Fhome%2Fcolas%2FDocuments%2Ftmp%2Fstrahler.map%26strahler%3Aprovider%3Dgdal' --INPUT1=0 --INPUT2='%3Fcrs%3DEPSG%3A2154%26extent%3D868999.5%2C6447000.5%2C885999.5%2C6467000.5%26width%3D17000%26height%3D20000%26formula%3D%2522strahler%401%2522%2520%253E%253D5%2520%26strahler%3Auri%3D%2Fhome%2Fcolas%2FDocuments%2Ftmp%2Fstrahler.map%26strahler%3Aprovider%3Dgdal' --OUTPUT=TEMPORARY_OUTPUT # # processing.run("pcraster:comparisonoperators", {'INPUT':'%3Fcrs%3DEPSG%3A2154%26extent%3D868999.5%2C6447000.5%2C885999.5%2C6467000.5%26width%3D17000%26height%3D20000%26formula%3D%2522strahler%401%2522%2520%253E%253D5%2520%26strahler%3Auri%3D%2Fhome%2Fcolas%2FDocuments%2Ftmp%2Fstrahler.map%26strahler%3Aprovider%3Dgdal','INPUT1':0,'INPUT2':'%3Fcrs%3DEPSG%3A2154%26extent%3D868999.5%2C6447000.5%2C885999.5%2C6467000.5%26width%3D17000%26height%3D20000%26formula%3D%2522strahler%401%2522%2520%253E%253D5%2520%26strahler%3Auri%3D%2Fhome%2Fcolas%2FDocuments%2Ftmp%2Fstrahler.map%26strahler%3Aprovider%3Dgdal','OUTPUT':'TEMPORARY_OUTPUT'}) # # processing.run("pcraster:col2map", {'INPUT':'/home/colas/Documents/tmp/bilan_sites.csv','INPUT1':"/media/colas/SRV/FICHIERS/SITES/DISTRICTS NATURELS/BASSE VALLEE DE L'ISERE/SONE/TUFI_Sone-à-Soi/SIG/MNT_1m.tif",'INPUT2':0,'OUTPUT':'TEMPORARY_OUTPUT'})