191 lines
8.6 KiB
Python
191 lines
8.6 KiB
Python
#!/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'}) |