Index: wflow-py/Scripts/wflow_prepare_step1.py =================================================================== diff -u -r073214c7c3a46cfcc2fdc2b27c5a9d101685c031 -r3e39e84af48f1bcb5ec0d243748147be223674f2 --- wflow-py/Scripts/wflow_prepare_step1.py (.../wflow_prepare_step1.py) (revision 073214c7c3a46cfcc2fdc2b27c5a9d101685c031) +++ wflow-py/Scripts/wflow_prepare_step1.py (.../wflow_prepare_step1.py) (revision 3e39e84af48f1bcb5ec0d243748147be223674f2) @@ -46,25 +46,25 @@ import numpy as np - def usage(*args): sys.stdout = sys.stderr - for msg in args: print(msg) + for msg in args: + print(msg) print(__doc__) sys.exit(0) -def configget(config,section,var,default): +def configget(config, section, var, default): """ gets parameter from config file and returns a default value if the parameter is not found """ try: - ret = config.get(section,var) + ret = config.get(section, var) except: print("returning default (" + default + ") for " + section + ":" + var) ret = default - + return ret @@ -77,13 +77,10 @@ else: print("Cannot open config file: " + fn) sys.exit(1) - + return config - - - def main(): """ @@ -96,267 +93,370 @@ strRiver = 8 masterdem = "dem.map" step1dir = "step1" - step2dir="step2" + step2dir = "step2" workdir = "." inifile = "wflow_prepare.ini" recreate = False snapgaugestoriver = False - + try: - opts, args = getopt.getopt(sys.argv[1:], 'W:hI:f') + opts, args = getopt.getopt(sys.argv[1:], "W:hI:f") except getopt.error as msg: usage(msg) - for o, a in opts: - if o == '-W': workdir = a - if o == '-I': inifile = a - if o == '-h': usage() - if o == '-f': recreate = True + if o == "-W": + workdir = a + if o == "-I": + inifile = a + if o == "-h": + usage() + if o == "-f": + recreate = True - tr.setglobaloption("unitcell") - os.chdir(workdir) + tr.setglobaloption("unitcell") + os.chdir(workdir) - config=OpenConf(workdir + "/" + inifile) + config = OpenConf(workdir + "/" + inifile) - masterdem = configget(config,"files","masterdem","dem.map") + masterdem = configget(config, "files", "masterdem", "dem.map") tr.setclone(masterdem) + strRiver = int(configget(config, "settings", "riverorder", "4")) - strRiver = int(configget(config,"settings","riverorder","4")) - try: - gauges_x = config.get("settings","gauges_x") - gauges_y = config.get("settings","gauges_y") + gauges_x = config.get("settings", "gauges_x") + gauges_y = config.get("settings", "gauges_y") except: print("gauges_x and gauges_y are required entries in the ini file") sys.exit(1) - step1dir = configget(config,"directories","step1dir","step1") - step2dir = configget(config,"directories","step2dir","step2") - #upscalefactor = float(config.get("settings","upscalefactor")) - - corevolume = float(configget(config,"settings","corevolume","1E35")) - catchmentprecipitation = float(configget(config,"settings","catchmentprecipitation","1E35")) - corearea = float(configget(config,"settings","corearea","1E35")) - outflowdepth = float(configget(config,"settings","lddoutflowdepth","1E35")) - - initialscale = int(configget(config,"settings","initialscale","1")) - csize= float(configget(config,"settings","cellsize","1")) + step1dir = configget(config, "directories", "step1dir", "step1") + step2dir = configget(config, "directories", "step2dir", "step2") + # upscalefactor = float(config.get("settings","upscalefactor")) - snapgaugestoriver=bool(int(configget(config,"settings","snapgaugestoriver","1"))) - lddglobaloption=configget(config,"settings","lddglobaloption","lddout") + corevolume = float(configget(config, "settings", "corevolume", "1E35")) + catchmentprecipitation = float( + configget(config, "settings", "catchmentprecipitation", "1E35") + ) + corearea = float(configget(config, "settings", "corearea", "1E35")) + outflowdepth = float(configget(config, "settings", "lddoutflowdepth", "1E35")) + + initialscale = int(configget(config, "settings", "initialscale", "1")) + csize = float(configget(config, "settings", "cellsize", "1")) + + snapgaugestoriver = bool( + int(configget(config, "settings", "snapgaugestoriver", "1")) + ) + lddglobaloption = configget(config, "settings", "lddglobaloption", "lddout") tr.setglobaloption(lddglobaloption) - lu_water= configget(config,"files","lu_water","") - lu_paved= configget(config,"files","lu_paved","") - + lu_water = configget(config, "files", "lu_water", "") + lu_paved = configget(config, "files", "lu_paved", "") + # X/Y coordinates of the gauges the system exec("X=np.array(" + gauges_x + ")") exec("Y=np.array(" + gauges_y + ")") - tr.Verbose=1 + tr.Verbose = 1 # make the directories to save results in - if not os.path.isdir(step1dir +"/"): + if not os.path.isdir(step1dir + "/"): os.makedirs(step1dir) if not os.path.isdir(step2dir): - os.makedirs(step2dir) + os.makedirs(step2dir) - if initialscale > 1: print("Initial scaling of DEM...") - os.system("resample -r " + str(initialscale) + " " + masterdem + " " + step1dir + "/dem_scaled.map") + os.system( + "resample -r " + + str(initialscale) + + " " + + masterdem + + " " + + step1dir + + "/dem_scaled.map" + ) print("Reading dem...") - dem = tr.readmap(step1dir + "/dem_scaled.map") - ldddem=dem + dem = tr.readmap(step1dir + "/dem_scaled.map") + ldddem = dem else: print("Reading dem...") dem = tr.readmap(masterdem) - ldddem=dem - - + ldddem = dem + try: - catchmask = config.get("files","catchment_mask") + catchmask = config.get("files", "catchment_mask") except: print("No catchment mask...") else: print("clipping DEM with mask.....") - mask=tr.readmap(catchmask) - ldddem = tr.ifthen(tr.boolean(mask),ldddem) - dem = tr.ifthen(tr.boolean(mask),dem) + mask = tr.readmap(catchmask) + ldddem = tr.ifthen(tr.boolean(mask), ldddem) + dem = tr.ifthen(tr.boolean(mask), dem) - - # See if there is a shape file of the river to burn in + # See if there is a shape file of the river to burn in try: - rivshp = config.get("files","river") + rivshp = config.get("files", "river") except: print("no river file specified") - outletpointX = float(configget(config,"settings","outflowpointX","0.0")) - outletpointY = float(configget(config,"settings","outflowpointY","0.0")) + outletpointX = float(configget(config, "settings", "outflowpointX", "0.0")) + outletpointY = float(configget(config, "settings", "outflowpointY", "0.0")) else: print("river file specified.....") try: - outletpointX = float(configget(config,"settings","outflowpointX","0.0")) - outletpointY = float(configget(config,"settings","outflowpointY","0.0")) + outletpointX = float(configget(config, "settings", "outflowpointX", "0.0")) + outletpointY = float(configget(config, "settings", "outflowpointY", "0.0")) except: - print("Need to specify the river outletpoint (a point at the end of the river within the current map)") + print( + "Need to specify the river outletpoint (a point at the end of the river within the current map)" + ) exit(1) - - outletpointmap = tr.points_to_map(dem,outletpointX,outletpointY,0.5) - tr.report(outletpointmap,step1dir + "/outletpoint.map") - #rivshpattr = config.get("files","riverattr") - tr.report(dem * 0.0,step1dir + "/nilmap.map") - thestr = "gdal_translate -of GTiff " + step1dir + "/nilmap.map " + step1dir + "/riverburn.tif" + + outletpointmap = tr.points_to_map(dem, outletpointX, outletpointY, 0.5) + tr.report(outletpointmap, step1dir + "/outletpoint.map") + # rivshpattr = config.get("files","riverattr") + tr.report(dem * 0.0, step1dir + "/nilmap.map") + thestr = ( + "gdal_translate -of GTiff " + + step1dir + + "/nilmap.map " + + step1dir + + "/riverburn.tif" + ) os.system(thestr) rivshpattr = os.path.splitext(os.path.basename(rivshp))[0] - os.system("gdal_rasterize -burn 1 -l " + rivshpattr + " " + rivshp + " " + step1dir + "/riverburn.tif") - thestr = "gdal_translate -of PCRaster " + step1dir + "/riverburn.tif " + step1dir + "/riverburn.map" + os.system( + "gdal_rasterize -burn 1 -l " + + rivshpattr + + " " + + rivshp + + " " + + step1dir + + "/riverburn.tif" + ) + thestr = ( + "gdal_translate -of PCRaster " + + step1dir + + "/riverburn.tif " + + step1dir + + "/riverburn.map" + ) os.system(thestr) riverburn = tr.readmap(step1dir + "/riverburn.map") # Determine regional slope assuming that is the way the river should run # Determine regional slope assuming that is the way the river should run - #tr.setglobaloption("unitcell") - #demregional=tr.windowaverage(dem,100) - ldddem = tr.ifthenelse(riverburn >= 1.0, dem -1000 , dem) - - tr.setglobaloption("unittrue") - upscalefactor=int(csize/tr.celllength()) + # tr.setglobaloption("unitcell") + # demregional=tr.windowaverage(dem,100) + ldddem = tr.ifthenelse(riverburn >= 1.0, dem - 1000, dem) + tr.setglobaloption("unittrue") + upscalefactor = int(csize / tr.celllength()) + print("Creating ldd...") - ldd=tr.lddcreate_save(step1dir +"/ldd.map",ldddem, recreate, outflowdepth=outflowdepth,corevolume=corevolume,catchmentprecipitation=catchmentprecipitation,corearea=corearea) + ldd = tr.lddcreate_save( + step1dir + "/ldd.map", + ldddem, + recreate, + outflowdepth=outflowdepth, + corevolume=corevolume, + catchmentprecipitation=catchmentprecipitation, + corearea=corearea, + ) print("Determining streamorder...") - stro=tr.streamorder(ldd) - tr.report(stro,step1dir + "/streamorder.map") + stro = tr.streamorder(ldd) + tr.report(stro, step1dir + "/streamorder.map") strdir = tr.ifthen(stro >= strRiver, stro) - tr.report(strdir,step1dir + "/streamorderrive.map") - tr.report(tr.boolean(tr.ifthen(stro >= strRiver, stro)),step1dir + "/rivers.map") + tr.report(strdir, step1dir + "/streamorderrive.map") + tr.report(tr.boolean(tr.ifthen(stro >= strRiver, stro)), step1dir + "/rivers.map") tr.setglobaloption("unittrue") # outlet (and other gauges if given) - #TODO: check is x/y set if not skip this + # TODO: check is x/y set if not skip this print("Outlet...") - outlmap = tr.points_to_map(dem,X,Y,0.5) + outlmap = tr.points_to_map(dem, X, Y, 0.5) if snapgaugestoriver: print("Snapping gauges to nearest river cells...") - tr.report(outlmap,step1dir + "/orggauges.map") - outlmap= tr.snaptomap(outlmap,strdir) + tr.report(outlmap, step1dir + "/orggauges.map") + outlmap = tr.snaptomap(outlmap, strdir) + # noutletmap = tr.points_to_map(dem,XX,YY,0.5) + # tr.report(noutletmap,'noutlet.map') - #noutletmap = tr.points_to_map(dem,XX,YY,0.5) - #tr.report(noutletmap,'noutlet.map') + tr.report(outlmap, step1dir + "/gauges.map") - - tr.report(outlmap,step1dir + "/gauges.map") - - # check if there is a pre-define catchment map + # check if there is a pre-define catchment map try: - catchmask = config.get("files","catchment_mask") + catchmask = config.get("files", "catchment_mask") except: print("No catchment mask, finding outlet") # Find catchment (overall) outlet = tr.find_outlet(ldd) - sub = tr.subcatch(ldd,outlet) - tr.report(sub,step1dir + "/catchment_overall.map") + sub = tr.subcatch(ldd, outlet) + tr.report(sub, step1dir + "/catchment_overall.map") else: print("reading and converting catchment mask.....") - os.system("resample -r " + str(initialscale) + " " + catchmask + " " + step1dir + "/catchment_overall.map") + os.system( + "resample -r " + + str(initialscale) + + " " + + catchmask + + " " + + step1dir + + "/catchment_overall.map" + ) sub = tr.readmap(step1dir + "/catchment_overall.map") print("Scatch...") - sd = tr.subcatch(ldd,tr.ifthen(outlmap>0,outlmap)) - tr.report(sd,step1dir + "/scatch.map") + sd = tr.subcatch(ldd, tr.ifthen(outlmap > 0, outlmap)) + tr.report(sd, step1dir + "/scatch.map") tr.setglobaloption("unitcell") print("Upscalefactor: " + str(upscalefactor)) - + if upscalefactor > 1: gc.collect() print("upscale river length1 (checkerboard map)...") - ck = tr.checkerboard(dem,upscalefactor) - tr.report(ck,step1dir + "/ck.map") - tr.report(dem,step1dir + "/demck.map") + ck = tr.checkerboard(dem, upscalefactor) + tr.report(ck, step1dir + "/ck.map") + tr.report(dem, step1dir + "/demck.map") print("upscale river length2...") - fact = tr.area_riverlength_factor(ldd, ck,upscalefactor) - tr.report(fact,step1dir + "/riverlength_fact.map") - - #print("make dem statistics...") - dem_ = tr.areaaverage(dem,ck) - tr.report(dem_,step1dir + "/demavg.map") - + fact = tr.area_riverlength_factor(ldd, ck, upscalefactor) + tr.report(fact, step1dir + "/riverlength_fact.map") + + # print("make dem statistics...") + dem_ = tr.areaaverage(dem, ck) + tr.report(dem_, step1dir + "/demavg.map") + print("Create DEM statistics...") - dem_ = tr.areaminimum(dem,ck) - tr.report(dem_,step1dir + "/demmin.map") - dem_ = tr.areamaximum(dem,ck) - tr.report(dem_,step1dir + "/demmax.map") + dem_ = tr.areaminimum(dem, ck) + tr.report(dem_, step1dir + "/demmin.map") + dem_ = tr.areamaximum(dem, ck) + tr.report(dem_, step1dir + "/demmax.map") # calculate percentiles - order = tr.areaorder(dem,ck) - n = tr.areatotal(tr.spatial(tr.scalar(1.0)),ck) + order = tr.areaorder(dem, ck) + n = tr.areatotal(tr.spatial(tr.scalar(1.0)), ck) #: calculate 25 percentile - perc = tr.area_percentile(dem,ck,n,order,25.0) - tr.report(perc,step1dir + "/dem25.map") - perc = tr.area_percentile(dem,ck,n,order,10.0) - tr.report(perc,step1dir + "/dem10.map") - perc = tr.area_percentile(dem,ck,n,order,50.0) - tr.report(perc,step1dir + "/dem50.map") - perc = tr.area_percentile(dem,ck,n,order,33.0) - tr.report(perc,step1dir + "/dem33.map") - perc = tr.area_percentile(dem,ck,n,order,66.0) - tr.report(perc,step1dir + "/dem66.map") - perc = tr.area_percentile(dem,ck,n,order,75.0) - tr.report(perc,step1dir + "/dem75.map") - perc = tr.area_percentile(dem,ck,n,order,90.0) - tr.report(perc,step1dir + "/dem90.map") + perc = tr.area_percentile(dem, ck, n, order, 25.0) + tr.report(perc, step1dir + "/dem25.map") + perc = tr.area_percentile(dem, ck, n, order, 10.0) + tr.report(perc, step1dir + "/dem10.map") + perc = tr.area_percentile(dem, ck, n, order, 50.0) + tr.report(perc, step1dir + "/dem50.map") + perc = tr.area_percentile(dem, ck, n, order, 33.0) + tr.report(perc, step1dir + "/dem33.map") + perc = tr.area_percentile(dem, ck, n, order, 66.0) + tr.report(perc, step1dir + "/dem66.map") + perc = tr.area_percentile(dem, ck, n, order, 75.0) + tr.report(perc, step1dir + "/dem75.map") + perc = tr.area_percentile(dem, ck, n, order, 90.0) + tr.report(perc, step1dir + "/dem90.map") else: - print("No fancy scaling done. Going strait to step2....") - tr.report(dem,step1dir + "/demavg.map") - Xul = float(config.get("settings","Xul")) - Yul = float(config.get("settings","Yul")) - Xlr = float(config.get("settings","Xlr")) - Ylr = float(config.get("settings","Ylr")) - gdalstr = "gdal_translate -projwin " + str(Xul) + " " + str(Yul) + " " +str(Xlr) + " " +str(Ylr) + " -of PCRaster " - #gdalstr = "gdal_translate -a_ullr " + str(Xul) + " " + str(Yul) + " " +str(Xlr) + " " +str(Ylr) + " -of PCRaster " - print(gdalstr) - tr.report(tr.cover(1.0),step1dir + "/wflow_riverlength_fact.map") - # Now us gdat tp convert the maps - os.system(gdalstr + step1dir + "/wflow_riverlength_fact.map" + " " + step2dir + "/wflow_riverlength_fact.map") - os.system(gdalstr + step1dir + "/demavg.map" + " " + step2dir + "/wflow_dem.map") - os.system(gdalstr + step1dir + "/demavg.map" + " " + step2dir + "/wflow_demmin.map") - os.system(gdalstr + step1dir + "/demavg.map" + " " + step2dir + "/wflow_demmax.map") - os.system(gdalstr + step1dir + "/gauges.map" + " " + step2dir + "/wflow_gauges.map") - os.system(gdalstr + step1dir + "/rivers.map" + " " + step2dir + "/wflow_river.map") - os.system(gdalstr + step1dir + "/streamorder.map" + " " + step2dir + "/wflow_streamorder.map") - os.system(gdalstr + step1dir + "/gauges.map" + " " + step2dir + "/wflow_outlet.map") - os.system(gdalstr + step1dir + "/scatch.map" + " " + step2dir + "/wflow_catchment.map") - os.system(gdalstr + step1dir + "/ldd.map" + " " + step2dir + "/wflow_ldd.map") - os.system(gdalstr + step1dir + "/scatch.map" + " " + step2dir + "/wflow_subcatch.map") - - if lu_water: - os.system(gdalstr + lu_water + " " + step2dir + "/WaterFrac.map") + print("No fancy scaling done. Going strait to step2....") + tr.report(dem, step1dir + "/demavg.map") + Xul = float(config.get("settings", "Xul")) + Yul = float(config.get("settings", "Yul")) + Xlr = float(config.get("settings", "Xlr")) + Ylr = float(config.get("settings", "Ylr")) + gdalstr = ( + "gdal_translate -projwin " + + str(Xul) + + " " + + str(Yul) + + " " + + str(Xlr) + + " " + + str(Ylr) + + " -of PCRaster " + ) + # gdalstr = "gdal_translate -a_ullr " + str(Xul) + " " + str(Yul) + " " +str(Xlr) + " " +str(Ylr) + " -of PCRaster " + print(gdalstr) + tr.report(tr.cover(1.0), step1dir + "/wflow_riverlength_fact.map") + # Now us gdat tp convert the maps + os.system( + gdalstr + + step1dir + + "/wflow_riverlength_fact.map" + + " " + + step2dir + + "/wflow_riverlength_fact.map" + ) + os.system( + gdalstr + step1dir + "/demavg.map" + " " + step2dir + "/wflow_dem.map" + ) + os.system( + gdalstr + step1dir + "/demavg.map" + " " + step2dir + "/wflow_demmin.map" + ) + os.system( + gdalstr + step1dir + "/demavg.map" + " " + step2dir + "/wflow_demmax.map" + ) + os.system( + gdalstr + step1dir + "/gauges.map" + " " + step2dir + "/wflow_gauges.map" + ) + os.system( + gdalstr + step1dir + "/rivers.map" + " " + step2dir + "/wflow_river.map" + ) + os.system( + gdalstr + + step1dir + + "/streamorder.map" + + " " + + step2dir + + "/wflow_streamorder.map" + ) + os.system( + gdalstr + step1dir + "/gauges.map" + " " + step2dir + "/wflow_outlet.map" + ) + os.system( + gdalstr + step1dir + "/scatch.map" + " " + step2dir + "/wflow_catchment.map" + ) + os.system(gdalstr + step1dir + "/ldd.map" + " " + step2dir + "/wflow_ldd.map") + os.system( + gdalstr + step1dir + "/scatch.map" + " " + step2dir + "/wflow_subcatch.map" + ) - if lu_paved: - os.system(gdalstr + lu_paved + " " + step2dir + "/PathFrac.map") + if lu_water: + os.system(gdalstr + lu_water + " " + step2dir + "/WaterFrac.map") - try: - lumap = config.get("files","landuse") - except: + if lu_paved: + os.system(gdalstr + lu_paved + " " + step2dir + "/PathFrac.map") + + try: + lumap = config.get("files", "landuse") + except: print("no landuse map...creating uniform map") - #clone=tr.readmap(step2dir + "/wflow_dem.map") + # clone=tr.readmap(step2dir + "/wflow_dem.map") tr.setclone(step2dir + "/wflow_dem.map") - tr.report(tr.nominal(1),step2dir + "/wflow_landuse.map") - else: - os.system("resample --clone " + step2dir + "/wflow_dem.map " + lumap + " " + step2dir + "/wflow_landuse.map") + tr.report(tr.nominal(1), step2dir + "/wflow_landuse.map") + else: + os.system( + "resample --clone " + + step2dir + + "/wflow_dem.map " + + lumap + + " " + + step2dir + + "/wflow_landuse.map" + ) - try: - soilmap = config.get("files","soil") - except: - print("no soil map..., creating uniform map") - tr.setclone(step2dir + "/wflow_dem.map") - tr.report(tr.nominal(1),step2dir + "/wflow_soil.map") - else: - os.system("resample --clone " + step2dir + "/wflow_dem.map " + soilmap + " " + step2dir + "/wflow_soil.map") + try: + soilmap = config.get("files", "soil") + except: + print("no soil map..., creating uniform map") + tr.setclone(step2dir + "/wflow_dem.map") + tr.report(tr.nominal(1), step2dir + "/wflow_soil.map") + else: + os.system( + "resample --clone " + + step2dir + + "/wflow_dem.map " + + soilmap + + " " + + step2dir + + "/wflow_soil.map" + ) if __name__ == "__main__":