diff --git a/.gitignore b/.gitignore index fd9da1c..bdbcab7 100644 --- a/.gitignore +++ b/.gitignore @@ -39,3 +39,10 @@ .idea **.idea .idea/* + +# Temporary files generated by ArcGIS +*.vat.cpg +*.vat.dbf +*.tif.aux.xml +*.tfw +*.tif.ovr \ No newline at end of file diff --git a/py_main/Config.py b/py_main/Config.py index 826eed3..2636e5e 100644 --- a/py_main/Config.py +++ b/py_main/Config.py @@ -102,15 +102,13 @@ VlySrc = None else: raise ValueError("The vlysrc %s is not existed or have no access permission!" % VlySrc) -else: - VlySrcCal = VlySrc + if not isFileExists(RdgSrc): if RdgSrc.lower() == 'none' or RdgSrc.lower() == '': RdgSrc = None else: raise ValueError("The RdgSrc %s is not existed or have no access permission!" % RdgSrc) -else: - RdgSrcCal = RdgSrc + if not isFileExists(rpiFile): if rpiFile.lower() == 'none' or rpiFile.lower() == '': rpiFile = None diff --git a/py_main/PreProcessing.py b/py_main/PreProcessing.py index 95f2b02..558edd0 100644 --- a/py_main/PreProcessing.py +++ b/py_main/PreProcessing.py @@ -1,171 +1,181 @@ -#! /usr/bin/env python -# coding=utf-8 -# @Description: Calculate terrain attributes from DEM and other optional inputs for deriving slope position. -# Slope, Curvature, RPI, HAND, etc. -# @Author: Liang-Jun Zhu -# -import time -from shutil import copy2 -import TauDEM -from Nomenclature import * -from Util import * -from RidgeExtraction import findRidge - - -def PreProcessing(model): - startT = time.time() - logStatus = open(log_preproc, 'w') - if model == 0: - logStatus.write("Preprocessing based on D8 flow model.\n") - elif model == 1: - logStatus.write("Preprocessing based on D-infinity flow model.\n") - logStatus.flush() - logStatus.write("[Preprocessing] [1/7] Converting DEM file format to GeoTiff...\n") - logStatus.flush() - Raster2GeoTIFF(rawdem, dem) - logStatus.write("[Preprocessing] [2/7] Removing pits...\n") - logStatus.flush() - TauDEM.pitremove(dem, inputProc, demfil, mpiexeDir = mpiexeDir, exeDir = exeDir, - hostfile = hostfile) # pitremove in TauDEM - logStatus.write("[Preprocessing] [3/7] Flow direction and slope in radian...\n") - logStatus.flush() - TauDEM.D8FlowDir(demfil, inputProc, D8FlowDir, D8Slp, mpiexeDir = mpiexeDir, exeDir = exeDir, hostfile = hostfile) - - if model == 1: - TauDEM.DinfFlowDir(demfil, inputProc, DinfFlowDir, DinfSlp, mpiexeDir = mpiexeDir, exeDir = exeDir, - hostfile = hostfile) - logStatus.write("[Preprocessing] [4/7] Generating flow accumulation...\n") - logStatus.flush() - # flow accumulation without weight grid or outlet - TauDEM.AreaD8(D8FlowDir, '', '', 'false', inputProc, D8ContriArea, mpiexeDir = mpiexeDir, exeDir = exeDir, - hostfile = hostfile) - global D8StreamThreshold - if D8StreamThreshold > 0: - TauDEM.Threshold(D8ContriArea, '', D8StreamThreshold, inputProc, D8Stream, mpiexeDir = mpiexeDir, - exeDir = exeDir, hostfile = hostfile) - else: - # initial stream - maxAccum, minAccum, meanAccum, STDAccum = GetRasterStat(D8ContriArea) - TauDEM.Threshold(D8ContriArea, '', meanAccum, inputProc, D8Stream, mpiexeDir = mpiexeDir, - exeDir = exeDir, hostfile = hostfile) - global outlet - if outlet is None: - TauDEM.ConnectDown(D8ContriArea, outletpre, inputProc, mpiexeDir = mpiexeDir, exeDir = exeDir, - hostfile = hostfile) - outlet = outletpre - TauDEM.MoveOutletsToStreams(D8FlowDir, D8Stream, outlet, maxMoveDist, inputProc, outletM, mpiexeDir = mpiexeDir, - exeDir = exeDir, hostfile = hostfile) - if model == 0: - logStatus.write("[Preprocessing] [5/7] Generating stream source raster based on Drop Analysis...\n") - elif model == 1: - logStatus.write( - "[Preprocessing] [5/7] Generating stream source raster based on Threshold derived from D8 flow model drop analysis or assigned...\n") - logStatus.flush() - if model == 1: - TauDEM.AreaDinf(DinfFlowDir, '', '', 'false', inputProc, DinfContriArea, mpiexeDir = mpiexeDir, exeDir = exeDir, - hostfile = hostfile) - # ReCalculate flow accumulation using PkrDglStream as weight grid - if D8StreamThreshold <= 0: - TauDEM.StreamSkeleton(demfil, PkrDglStream, inputProc, mpiexeDir = mpiexeDir, exeDir = exeDir, - hostfile = hostfile) - TauDEM.AreaD8(D8FlowDir, '', PkrDglStream, 'false', inputProc, D8ContriArea, mpiexeDir = mpiexeDir, - exeDir = exeDir, hostfile = hostfile) - maxAccum, minAccum, meanAccum, STDAccum = GetRasterStat(D8ContriArea) - - if meanAccum < STDAccum: - minthresh = meanAccum - else: - minthresh = meanAccum - STDAccum - maxthresh = meanAccum + STDAccum - - TauDEM.DropAnalysis(demfil, D8FlowDir, D8ContriArea, D8ContriArea, outletM, minthresh, maxthresh, numthresh, - logspace, inputProc, drpFile, mpiexeDir = mpiexeDir, exeDir = exeDir, - hostfile = hostfile) - drpf = open(drpFile, "r") - tempContents = drpf.read() - (beg, d8drpThreshold) = tempContents.rsplit(' ', 1) - drpf.close() - D8StreamThreshold = d8drpThreshold - - TauDEM.Threshold(D8ContriArea, '', D8StreamThreshold, inputProc, D8Stream, - mpiexeDir = mpiexeDir, exeDir = exeDir, hostfile = hostfile) - if model == 1: - global DinfStreamThreshold - if DinfStreamThreshold == 0: - DinfStreamThreshold = D8StreamThreshold - TauDEM.Threshold(DinfContriArea, '', DinfStreamThreshold, inputProc, DinfStream, mpiexeDir = mpiexeDir, - exeDir = exeDir, hostfile = hostfile) - logStatus.write("[Preprocessing] [6/7] Calculating RPI(Relative Position Index)...\n") - logStatus.flush() - if model == 0: # D8 model - # HAND - TauDEM.D8DistDownToStream(D8FlowDir, demfil, D8Stream, D8DistDown_V, 'vertical', D8StreamTag, inputProc, - mpiexeDir = mpiexeDir, exeDir = exeDir, hostfile = hostfile) - if rpiMethod == 1: # calculate RPI based on hydrological proximity measures (Default). - TauDEM.D8DistDownToStream(D8FlowDir, demfil, D8Stream, D8DistDown, D8DownMethod, D8StreamTag, inputProc, - mpiexeDir = mpiexeDir, exeDir = exeDir, hostfile = hostfile) - TauDEM.D8DistUpToRidge(D8FlowDir, demfil, D8DistUp, D8UpMethod, D8UpStats, inputProc, rdg = RdgSrc, - mpiexeDir = mpiexeDir, exeDir = exeDir, hostfile = hostfile) - TauDEM.SimpleCalculator(D8DistDown, D8DistUp, RPID8, 4, inputProc, mpiexeDir = mpiexeDir, exeDir = exeDir, - hostfile = hostfile) - elif model == 1: # Dinf model - # HAND - TauDEM.DinfDistDown(DinfFlowDir, demfil, D8Stream, DinfDownStat, 'vertical', 'false', DinfDistDownWG, - inputProc, DinfDistDown_V, mpiexeDir = mpiexeDir, exeDir = exeDir, hostfile = hostfile) - if rpiMethod == 1: - TauDEM.DinfDistDown(DinfFlowDir, demfil, D8Stream, DinfDownStat, DinfDownMethod, 'false', DinfDistDownWG, - inputProc, DinfDistDown, mpiexeDir = mpiexeDir, exeDir = exeDir, hostfile = hostfile) - TauDEM.DinfDistUpToRidge(DinfFlowDir, demfil, DinfSlp, propthresh, DinfUpStat, DinfUpMethod, 'false', - inputProc, DinfDistUp, rdg = RdgSrc, mpiexeDir = mpiexeDir, exeDir = exeDir, - hostfile = hostfile) - TauDEM.SimpleCalculator(DinfDistDown, DinfDistUp, RPIDinf, 4, inputProc, mpiexeDir = mpiexeDir, - exeDir = exeDir, hostfile = hostfile) - if rpiMethod == 0: # calculate RPI based on Skidmore's method - TauDEM.StreamNet(demfil, D8FlowDir, D8ContriArea, D8Stream, outletM, D8StreamOrd, NetTree, NetCoord, - D8StreamNet, SubBasin, inputProc, mpiexeDir = mpiexeDir, exeDir = exeDir, hostfile = hostfile) - if VlySrcCal is None or not isFileExists(VlySrcCal): - copy2(D8Stream, VlySrcCal) - if RdgSrcCal is None or not isFileExists(RdgSrcCal): - # C++ version - angfile = D8FlowDir - elevfile = D8DistDown_V - if model == 1: # D-inf model - angfile = DinfFlowDir - elevfile = DinfDistDown_V - TauDEM.ExtractRidge(angfile, elevfile, RdgSrcCal, inputProc, mpiexeDir = mpiexeDir, - exeDir = exeDir, hostfile = hostfile) - # findRidge(1, RdgSrcCal) # Python-version - TauDEM.RPISkidmore(VlySrcCal, RdgSrcCal, RPISkidmore, inputProc, 1, 1, dist2Vly, dist2Rdg, - mpiexeDir = mpiexeDir, exeDir = exeDir, hostfile = hostfile) - logStatus.write("[Preprocessing] [7/7] Calculating Plan Curvature and Profile Curvature...\n") - logStatus.flush() - TauDEM.Curvature(inputProc, demfil, prof = ProfC, horiz = HorizC, mpiexeDir = mpiexeDir, exeDir = exeDir, - hostfile = hostfile) - - if model == 0: - slopeTrans(D8Slp, Slope) - if rpiMethod == 1: - copy2(RPID8, RPI) - else: - copy2(RPISkidmore, RPI) - copy2(D8DistDown_V, HAND) - elif model == 1: - slopeTrans(DinfSlp, Slope) - if rpiMethod == 1: - copy2(RPIDinf, RPI) - else: - copy2(RPISkidmore, RPI) - copy2(DinfDistDown_V, HAND) - logStatus.write("[Preprocessing] Preprocessing succeed!\n") - logStatus.flush() - endT = time.time() - cost = (endT - startT) / 60. - logStatus.write("Time consuming: %.2f min.\n" % cost) - logStatus.flush() - logStatus.close() - - -if __name__ == '__main__': - ini, proc, root = GetInputArgs() - LoadConfiguration(ini, proc, root) - PreProcessing(FlowModel) +#! /usr/bin/env python +# coding=utf-8 +# @Description: Calculate terrain attributes from DEM and other optional inputs for deriving slope position. +# Slope, Curvature, RPI, HAND, etc. +# @Author: Liang-Jun Zhu +# +import time +from shutil import copy2 +import TauDEM +from Nomenclature import * +from Util import * +from RidgeExtraction import findRidge + + +def PreProcessing(model): + startT = time.time() + logStatus = open(log_preproc, 'w') + # if valley and ridge are provided + global VlySrcCal + global RdgSrcCal + if VlySrc is not None: + VlySrcCal = VlySrc + if RdgSrc is not None: + RdgSrcCal = RdgSrc + if model == 0: + logStatus.write("Preprocessing based on D8 flow model.\n") + elif model == 1: + logStatus.write("Preprocessing based on D-infinity flow model.\n") + logStatus.flush() + logStatus.write("[Preprocessing] [1/7] Converting DEM file format to GeoTiff...\n") + logStatus.flush() + Raster2GeoTIFF(rawdem, dem) + logStatus.write("[Preprocessing] [2/7] Removing pits...\n") + logStatus.flush() + TauDEM.pitremove(dem, inputProc, demfil, mpiexeDir = mpiexeDir, exeDir = exeDir, + hostfile = hostfile) # pitremove in TauDEM + logStatus.write("[Preprocessing] [3/7] Flow direction and slope in radian...\n") + logStatus.flush() + TauDEM.D8FlowDir(demfil, inputProc, D8FlowDir, D8Slp, mpiexeDir = mpiexeDir, exeDir = exeDir, hostfile = hostfile) + + if model == 1: + TauDEM.DinfFlowDir(demfil, inputProc, DinfFlowDir, DinfSlp, mpiexeDir = mpiexeDir, exeDir = exeDir, + hostfile = hostfile) + logStatus.write("[Preprocessing] [4/7] Generating flow accumulation...\n") + logStatus.flush() + # flow accumulation without weight grid or outlet + TauDEM.AreaD8(D8FlowDir, '', '', 'false', inputProc, D8ContriArea, mpiexeDir = mpiexeDir, exeDir = exeDir, + hostfile = hostfile) + global D8StreamThreshold + if D8StreamThreshold > 0: + TauDEM.Threshold(D8ContriArea, '', D8StreamThreshold, inputProc, D8Stream, mpiexeDir = mpiexeDir, + exeDir = exeDir, hostfile = hostfile) + else: + # initial stream + maxAccum, minAccum, meanAccum, STDAccum = GetRasterStat(D8ContriArea) + TauDEM.Threshold(D8ContriArea, '', meanAccum, inputProc, D8Stream, mpiexeDir = mpiexeDir, + exeDir = exeDir, hostfile = hostfile) + global outlet + if outlet is None: + TauDEM.ConnectDown(D8ContriArea, outletpre, inputProc, mpiexeDir = mpiexeDir, exeDir = exeDir, + hostfile = hostfile) + outlet = outletpre + TauDEM.MoveOutletsToStreams(D8FlowDir, D8Stream, outlet, maxMoveDist, inputProc, outletM, mpiexeDir = mpiexeDir, + exeDir = exeDir, hostfile = hostfile) + if model == 0: + logStatus.write("[Preprocessing] [5/7] Generating stream source raster based on Drop Analysis...\n") + elif model == 1: + logStatus.write( + "[Preprocessing] [5/7] Generating stream source raster based on Threshold derived from D8 flow model drop analysis or assigned...\n") + logStatus.flush() + if model == 1: + TauDEM.AreaDinf(DinfFlowDir, '', '', 'false', inputProc, DinfContriArea, mpiexeDir = mpiexeDir, exeDir = exeDir, + hostfile = hostfile) + # ReCalculate flow accumulation using PkrDglStream as weight grid + if D8StreamThreshold <= 0: + TauDEM.StreamSkeleton(demfil, PkrDglStream, inputProc, mpiexeDir = mpiexeDir, exeDir = exeDir, + hostfile = hostfile) + TauDEM.AreaD8(D8FlowDir, '', PkrDglStream, 'false', inputProc, D8ContriArea, mpiexeDir = mpiexeDir, + exeDir = exeDir, hostfile = hostfile) + maxAccum, minAccum, meanAccum, STDAccum = GetRasterStat(D8ContriArea) + + if meanAccum < STDAccum: + minthresh = meanAccum + else: + minthresh = meanAccum - STDAccum + maxthresh = meanAccum + STDAccum + + TauDEM.DropAnalysis(demfil, D8FlowDir, D8ContriArea, D8ContriArea, outletM, minthresh, maxthresh, numthresh, + logspace, inputProc, drpFile, mpiexeDir = mpiexeDir, exeDir = exeDir, + hostfile = hostfile) + drpf = open(drpFile, "r") + tempContents = drpf.read() + (beg, d8drpThreshold) = tempContents.rsplit(' ', 1) + drpf.close() + D8StreamThreshold = d8drpThreshold + + TauDEM.Threshold(D8ContriArea, '', D8StreamThreshold, inputProc, D8Stream, + mpiexeDir = mpiexeDir, exeDir = exeDir, hostfile = hostfile) + if model == 1: + global DinfStreamThreshold + if DinfStreamThreshold == 0: + DinfStreamThreshold = D8StreamThreshold + TauDEM.Threshold(DinfContriArea, '', DinfStreamThreshold, inputProc, DinfStream, mpiexeDir = mpiexeDir, + exeDir = exeDir, hostfile = hostfile) + logStatus.write("[Preprocessing] [6/7] Calculating RPI(Relative Position Index)...\n") + logStatus.flush() + StreamSource = D8Stream + if VlySrcCal is not None and isFileExists(VlySrcCal): + StreamSource = VlySrcCal + if model == 0: # D8 model + # HAND + TauDEM.D8DistDownToStream(D8FlowDir, demfil, StreamSource, D8DistDown_V, 'vertical', D8StreamTag, inputProc, + mpiexeDir = mpiexeDir, exeDir = exeDir, hostfile = hostfile) + if rpiMethod == 1: # calculate RPI based on hydrological proximity measures (Default). + TauDEM.D8DistDownToStream(D8FlowDir, demfil, StreamSource, D8DistDown, D8DownMethod, D8StreamTag, inputProc, + mpiexeDir = mpiexeDir, exeDir = exeDir, hostfile = hostfile) + TauDEM.D8DistUpToRidge(D8FlowDir, demfil, D8DistUp, D8UpMethod, D8UpStats, inputProc, rdg = RdgSrc, + mpiexeDir = mpiexeDir, exeDir = exeDir, hostfile = hostfile) + TauDEM.SimpleCalculator(D8DistDown, D8DistUp, RPID8, 4, inputProc, mpiexeDir = mpiexeDir, exeDir = exeDir, + hostfile = hostfile) + elif model == 1: # Dinf model + # HAND + TauDEM.DinfDistDown(DinfFlowDir, demfil, StreamSource, DinfDownStat, 'vertical', 'false', DinfDistDownWG, + inputProc, DinfDistDown_V, mpiexeDir = mpiexeDir, exeDir = exeDir, hostfile = hostfile) + if rpiMethod == 1: + TauDEM.DinfDistDown(DinfFlowDir, demfil, StreamSource, DinfDownStat, DinfDownMethod, 'false', DinfDistDownWG, + inputProc, DinfDistDown, mpiexeDir = mpiexeDir, exeDir = exeDir, hostfile = hostfile) + TauDEM.DinfDistUpToRidge(DinfFlowDir, demfil, DinfSlp, propthresh, DinfUpStat, DinfUpMethod, 'false', + inputProc, DinfDistUp, rdg = RdgSrc, mpiexeDir = mpiexeDir, exeDir = exeDir, + hostfile = hostfile) + TauDEM.SimpleCalculator(DinfDistDown, DinfDistUp, RPIDinf, 4, inputProc, mpiexeDir = mpiexeDir, + exeDir = exeDir, hostfile = hostfile) + if rpiMethod == 0: # calculate RPI based on Skidmore's method + TauDEM.StreamNet(demfil, D8FlowDir, D8ContriArea, D8Stream, outletM, D8StreamOrd, NetTree, NetCoord, + D8StreamNet, SubBasin, inputProc, mpiexeDir = mpiexeDir, exeDir = exeDir, hostfile = hostfile) + if VlySrcCal is None or not isFileExists(VlySrcCal): + copy2(D8Stream, VlySrcCal) + if RdgSrcCal is None or not isFileExists(RdgSrcCal): + # C++ version + angfile = D8FlowDir + elevfile = D8DistDown_V + if model == 1: # D-inf model + angfile = DinfFlowDir + elevfile = DinfDistDown_V + TauDEM.ExtractRidge(angfile, elevfile, RdgSrcCal, inputProc, mpiexeDir = mpiexeDir, + exeDir = exeDir, hostfile = hostfile) + # findRidge(1, RdgSrcCal) # Python-version + TauDEM.RPISkidmore(VlySrcCal, RdgSrcCal, RPISkidmore, inputProc, 1, 1, dist2Vly, dist2Rdg, + mpiexeDir = mpiexeDir, exeDir = exeDir, hostfile = hostfile) + logStatus.write("[Preprocessing] [7/7] Calculating Plan Curvature and Profile Curvature...\n") + logStatus.flush() + TauDEM.Curvature(inputProc, demfil, prof = ProfC, horiz = HorizC, mpiexeDir = mpiexeDir, exeDir = exeDir, + hostfile = hostfile) + + if model == 0: + slopeTrans(D8Slp, Slope) + if rpiMethod == 1: + copy2(RPID8, RPI) + else: + copy2(RPISkidmore, RPI) + copy2(D8DistDown_V, HAND) + elif model == 1: + slopeTrans(DinfSlp, Slope) + if rpiMethod == 1: + copy2(RPIDinf, RPI) + else: + copy2(RPISkidmore, RPI) + copy2(DinfDistDown_V, HAND) + logStatus.write("[Preprocessing] Preprocessing succeed!\n") + logStatus.flush() + endT = time.time() + cost = (endT - startT) / 60. + logStatus.write("Time consuming: %.2f min.\n" % cost) + logStatus.flush() + logStatus.close() + + +if __name__ == '__main__': + ini, proc, root = GetInputArgs() + LoadConfiguration(ini, proc, root) + PreProcessing(FlowModel) diff --git a/py_main/pbs_dgpm.sh b/py_main/pbs_dgpm.sh new file mode 100644 index 0000000..dd262fa --- /dev/null +++ b/py_main/pbs_dgpm.sh @@ -0,0 +1,16 @@ +#!/usr/bin/env bash +## This is an example for dgpm-cluster with CentOS 6.2. +#PBS -N pitremoveJob +#PBS -l nodes=1:ppn=4 +#PBS -j oe +#PBS -m abe -M zlj@lreis.ac.cn + +# Path to your executable. +cd $HOME/data + +# Add any addition to your environment variables like PATH. +# For example, if your local MPI installation is in $HOME/mpich/bin +export PATH=$HOME/mpich/bin:$PATH + +# Launch program using the hosts +mpiexec -n 4 $HOME/TauDEM5.3/exec/pitremove -z "pvdem.tif" -fel "pvdemfel.tif" \ No newline at end of file