Index: wflow-py/Scripts/wflow_prepare_step1.py =================================================================== diff -u -r58f139b2afc2424fc9f8796934878bd6e49c31de -r12ea40dc08628f654753679e0972e87b7bb12f7a --- wflow-py/Scripts/wflow_prepare_step1.py (.../wflow_prepare_step1.py) (revision 58f139b2afc2424fc9f8796934878bd6e49c31de) +++ wflow-py/Scripts/wflow_prepare_step1.py (.../wflow_prepare_step1.py) (revision 12ea40dc08628f654753679e0972e87b7bb12f7a) @@ -40,16 +40,17 @@ import os import os.path import getopt -import configparser +import ConfigParser import sys import gc import numpy as np def usage(*args): sys.stdout = sys.stderr - for msg in args: print(msg) - print(__doc__) + for msg in args: + print msg + print __doc__ sys.exit(0) @@ -61,20 +62,20 @@ try: ret = config.get(section, var) except: - print("returning default (" + default + ") for " + section + ":" + var) + print "returning default (" + default + ") for " + section + ":" + var ret = default return ret def OpenConf(fn): - config = configparser.SafeConfigParser() + config = ConfigParser.SafeConfigParser() config.optionxform = str if os.path.exists(fn): config.read(fn) else: - print("Cannot open config file: " + fn) + print "Cannot open config file: " + fn sys.exit(1) return config @@ -99,8 +100,8 @@ snapgaugestoriver = False try: - opts, args = getopt.getopt(sys.argv[1:], 'W:hI:f') - except getopt.error as msg: + opts, args = getopt.getopt(sys.argv[1:], "W:hI:f") + except getopt.error, msg: usage(msg) for o, a in opts: @@ -127,7 +128,7 @@ 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") + print "gauges_x and gauges_y are required entries in the ini file" sys.exit(1) step1dir = configget(config, "directories", "step1dir", "step1") @@ -153,8 +154,8 @@ 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 + ")") + exec "X=np.array(" + gauges_x + ")" + exec "Y=np.array(" + gauges_y + ")" tr.Verbose = 1 @@ -165,11 +166,19 @@ os.makedirs(step2dir) if initialscale > 1: - print("Initial scaling of DEM...") - os.system("resample -r " + str(initialscale) + " " + masterdem + " " + step1dir + "/dem_scaled.map") - print("Reading dem...") - dem = tr.readmap(step1dir + "/dem_scaled.map") - ldddem=dem + print "Initial scaling of DEM..." + os.system( + "resample -r " + + str(initialscale) + + " " + + masterdem + + " " + + step1dir + + "/dem_scaled.map" + ) + print ("Reading dem...") + dem = tr.readmap(step1dir + "/dem_scaled.map") + ldddem = dem else: print ("Reading dem...") dem = tr.readmap(masterdem) @@ -178,22 +187,22 @@ try: catchmask = config.get("files", "catchment_mask") except: - print("No catchment mask...") + 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) + print "clipping DEM with mask....." + 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 try: 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")) + print "no river file specified" + outletpointX = float(configget(config, "settings", "outflowpointX", "0.0")) + outletpointY = float(configget(config, "settings", "outflowpointY", "0.0")) else: - print("river file specified.....") + print "river file specified....." try: outletpointX = float(configget(config, "settings", "outflowpointX", "0.0")) outletpointY = float(configget(config, "settings", "outflowpointY", "0.0")) @@ -269,11 +278,10 @@ 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) + print "Snapping gauges to nearest river cells..." + 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') @@ -283,23 +291,31 @@ try: catchmask = config.get("files", "catchment_mask") except: - print("No catchment mask, finding outlet") + 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") else: - print("reading and converting catchment mask.....") - os.system("resample -r " + str(initialscale) + " " + catchmask + " " + step1dir + "/catchment_overall.map") + print "reading and converting catchment mask....." + 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") tr.setglobaloption("unitcell") - print("Upscalefactor: " + str(upscalefactor)) - + print "Upscalefactor: " + str(upscalefactor) + if upscalefactor > 1: gc.collect() print ("upscale river length1 (checkerboard map)...") @@ -338,53 +354,109 @@ 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: - print("no landuse map...creating uniform map") - #clone=tr.readmap(step2dir + "/wflow_dem.map") + 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") 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__":