CEN38/Python_scripts
2
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

init depot renew

Browse Source
This commit is contained in:
Colas
2026-09-22 11:58:57 +02:00
commit e9fa6a7309
205 changed files with 39721 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
hydro_analyse.py
+191
View File
@@ -0,0 +1,191 @@
#!/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/cgeier/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/cgeier/SRV/FICHIERS/SITES/DISTRICTS NATURELS/BASSE VALLEE DE L'ISERE/SONE/TUFI_Sone-à-Soi/SIG/MNT_1m.tif",'INPUT2':0,'OUTPUT':'TEMPORARY_OUTPUT'})