Fisheye: Tag 97c6554831d691a3f3d1298dfb70c96808dae45e refers to a dead (removed) revision in file `wflow-py/Scripts/wtools_py/CheckInput.py'. Fisheye: No comparison available. Pass `N' to diff? Fisheye: Tag 97c6554831d691a3f3d1298dfb70c96808dae45e refers to a dead (removed) revision in file `wflow-py/Scripts/wtools_py/CreateGrid.py'. Fisheye: No comparison available. Pass `N' to diff? Fisheye: Tag 97c6554831d691a3f3d1298dfb70c96808dae45e refers to a dead (removed) revision in file `wflow-py/Scripts/wtools_py/StaticMaps.py'. Fisheye: No comparison available. Pass `N' to diff? Index: wflow-py/Scripts/wtools_py/create_grid.py =================================================================== diff -u -r405114bf53fffa2637a5383483a1fb99fbc73f39 -r97c6554831d691a3f3d1298dfb70c96808dae45e --- wflow-py/Scripts/wtools_py/create_grid.py (.../create_grid.py) (revision 405114bf53fffa2637a5383483a1fb99fbc73f39) +++ wflow-py/Scripts/wtools_py/create_grid.py (.../create_grid.py) (revision 97c6554831d691a3f3d1298dfb70c96808dae45e) @@ -36,28 +36,27 @@ import wflow.wflowtools_lib as wt -def main(): - ### Read input arguments ##### - logfilename = 'wtools_create_grid.log' +def parse_args(): + ### Read input arguments ##### parser = OptionParser() usage = "usage: %prog [options]" parser = OptionParser(usage=usage) parser.add_option('-q', '--quiet', dest='verbose', default=True, action='store_false', help='do not print status messages to stdout') - parser.add_option('-f', '--file', dest='inputfile', nargs=1, + parser.add_option('-f', '--file', dest='inputfile', nargs=1, help='file of which extent will be read. Most logically the catchment layer\nformat: ESRI Shapefile or any gdal supported raster format (preferred GeoTiff)') - parser.add_option('-e', '--extent', - nargs=4, dest='extent', type='float', - help='extent in WGS 1984 lat-lon (xmin, ymin, xmax, ymax)') +# parser.add_option('-e', '--extent', +# nargs=4, dest='extent', type='float', +# help='extent in WGS 1984 lat-lon (xmin, ymin, xmax, ymax)') parser.add_option('-l', '--logfile', dest='logfilename', default='wtools_create_grid.log', help='log file name') parser.add_option('-p', '--projection', dest='projection', default='EPSG:4326', help='Only used if no file is provided, either of type EPSG:<####> or +proj...') parser.add_option('-c', '--cellsize', type='float', - nargs=1, dest='cellsize', default=0.005, + nargs=1, dest='cellsize', help='extent') parser.add_option('-s', '--snap', dest='snap', default=False, action='store_true', @@ -67,29 +66,50 @@ help='Destination folder (default=./wflow)') (options, args) = parser.parse_args() - - ##### Preprocessing ##### + + print options.__dict__.items() + + + ##### Preprocessing ##### # check if either a file or an extent is provided. If not, sys.exit - if not np.logical_or(os.path.isfile(options.inputfile), options.extent is not None): - parser.error('No file or extent given') + + + if options.inputfile is None and options.extent is None: + parser.error('No input file (-f filename) or extent (-e (xmin, ymin, xmax, ymax)) given') parser.print_help() sys.exit(1) + + if not options.inputfile is None: + if not os.path.exists(options.inputfile): + parser.error('input file provided but not found, please check path') + parser.print_help() + sys.exit(1) + + if options.cellsize is None: + parser.error('no cell size (-c cellsize) provided') + parser.print_help() + sys.exit(1) + + + return options +def main(logfilename,destination,inputfile,projection,cellsize,snap=False,verbose=True,locationid='wflow_mask'): + # open a logger, dependent on verbose print to screen or not - logger, ch = wt.setlogger(logfilename, 'WTOOLS', options.verbose) + logger, ch = wt.setlogger(logfilename, 'WTOOLS', verbose) # delete old files - if os.path.isdir(options.destination): - shutil.rmtree(options.destination) - os.makedirs(options.destination) + if os.path.isdir(destination): + shutil.rmtree(destination) + os.makedirs(destination) ### Get information #### - if options.inputfile is not None: + if inputfile is not None: # retrieve extent from input file. Check if projection is provided - file_ext = os.path.splitext(os.path.basename(options.inputfile))[1] + file_ext = os.path.splitext(os.path.basename(inputfile))[1] if file_ext == '.shp': - file_att = os.path.splitext(os.path.basename(options.inputfile))[0] - ds = ogr.Open(options.inputfile) + file_att = os.path.splitext(os.path.basename(inputfile))[0] + ds = ogr.Open(inputfile) # read the extent of the shapefile lyr = ds.GetLayerByName(file_att) extent = lyr.GetExtent() @@ -99,16 +119,16 @@ else: # Read extent from a GDAL compatible file try: - extent_in = wt.get_extent(options.inputfile) + extent_in = wt.get_extent(inputfile) except: msg = 'Input file {:s} not a shape or gdal file'.format( - options.inputfile) + inputfile) wt.close_with_error(logger, ch, msg) sys.exit(1) # # get spatial reference from grid file try: - srs = wt.get_projection(options.inputfile) + srs = wt.get_projection(inputfile) except: logger.warning( 'No projection found, assuming WGS 1984 lat long') @@ -128,20 +148,20 @@ # srs = osr.SpatialReference() # srs.ImportFromWkt(WktString) else: - lonmin, latmin, lonmax, latmax = options.extent + lonmin, latmin, lonmax, latmax = extent srs_4326 = osr.SpatialReference() srs_4326.ImportFromEPSG(4326) srs = osr.SpatialReference() - if options.projection is not None: + if projection is not None: # import projection as an srs object - if options.projection.lower()[0:4] == 'epsg': + if projection.lower()[0:4] == 'epsg': # make a proj4 string - srs.ImportFromEPSG(int(options.projection[5:])) - elif options.projection.lower()[0:5] == '+proj': - srs.ImportFromProj4(options.projection) + srs.ImportFromEPSG(int(projection[5:])) + elif projection.lower()[0:5] == '+proj': + srs.ImportFromProj4(projection) else: msg = 'Projection "{:s}" is not a valid projection'.format( - options.projection) + projection) wt.close_with_error(logger, ch, msg) else: logger.warning('No projection found, assuming WGS 1984 lat long') @@ -166,56 +186,56 @@ else: geodatum = 'WGS 1984' - if options.snap: + if snap: logger.info('Snapping raster') - snap = len(str(options.cellsize - np.floor(options.cellsize))) - 2 - extent_out = wt.round_extent(extent_in, options.cellsize, snap) + snap = len(str(cellsize - np.floor(cellsize))) - 2 + extent_out = wt.round_extent(extent_in, cellsize, snap) else: extent_out = extent_in - cols = int((extent_out[2] - extent_out[0]) / options.cellsize) # +2) - rows = int((extent_out[3] - extent_out[1]) / options.cellsize) # +2) + cols = int((extent_out[2] - extent_out[0]) / cellsize) # +2) + rows = int((extent_out[3] - extent_out[1]) / cellsize) # +2) cells = rows * cols - xorg = extent_out[0] # -options.cellsize - yorg = extent_out[3] # +options.cellsize + xorg = extent_out[0] # -cellsize + yorg = extent_out[3] # +cellsize # create clone raster print('rows: {0} cols: {1}'.format(rows, cols)) dummy_raster = np.zeros((rows, cols)) - 9999. clone_file_map = os.path.abspath( - os.path.join(options.destination, 'mask.map')) + os.path.join(destination, 'mask.map')) clone_file_tif = os.path.abspath( - os.path.join(options.destination, 'mask.tif')) + os.path.join(destination, 'mask.tif')) logger.info('Writing PCRaster clone to {:s}'.format(clone_file_map)) gis.gdal_writemap(clone_file_map, 'PCRaster', xorg, yorg, dummy_raster, - -9999., resolution=options.cellsize, + -9999., resolution=cellsize, srs=srs) logger.info('Writing Geotiff clone to {:s}'.format(clone_file_tif)) gis.gdal_writemap(clone_file_tif, 'GTiff', xorg, yorg, dummy_raster, - -9999., resolution=options.cellsize, + -9999., resolution=cellsize, zlib=True, srs=srs) # create grid.xml - root = etree.Element('regular', locationId='wflow_mask') + root = etree.Element('regular', locationId=locationid) etree.SubElement(root, 'rows').text = str(rows) etree.SubElement(root, 'columns').text = str(cols) etree.SubElement(root, 'geoDatum').text = geodatum etree.SubElement(root, 'firstCellCenter') - etree.SubElement(root[3], 'x').text = str(xorg + 0.5 * options.cellsize) - etree.SubElement(root[3], 'y').text = str(yorg - 0.5 * options.cellsize) - etree.SubElement(root, 'xCellSize').text = str(options.cellsize) - etree.SubElement(root, 'yCellSize').text = str(options.cellsize) - xml_file = os.path.abspath(os.path.join(options.destination, 'grid.xml')) + etree.SubElement(root[3], 'x').text = str(xorg + 0.5 * cellsize) + etree.SubElement(root[3], 'y').text = str(yorg - 0.5 * cellsize) + etree.SubElement(root, 'xCellSize').text = str(cellsize) + etree.SubElement(root, 'yCellSize').text = str(cellsize) + xml_file = os.path.abspath(os.path.join(destination, 'grid.xml')) logger.info('Writing FEWS grid definition to {:s}'.format(xml_file)) gridxml = open(xml_file, 'w+') gridxml.write(etree.tostring(root, pretty_print=True)) gridxml.close() # create shape file Driver = ogr.GetDriverByName("ESRI Shapefile") - shp_file = os.path.abspath(os.path.join(options.destination, 'mask.shp')) + shp_file = os.path.abspath(os.path.join(destination, 'mask.shp')) logger.info('Writing shape of clone to {:s}'.format(shp_file)) shp_att = os.path.splitext(os.path.basename(shp_file))[0] shp = Driver.CreateDataSource(shp_file) @@ -225,11 +245,11 @@ lyr.CreateField(fieldDef) ring = ogr.Geometry(ogr.wkbLinearRing) ring.AddPoint(xorg, yorg) - ring.AddPoint(xorg + cols * options.cellsize, + ring.AddPoint(xorg + cols * cellsize, yorg) - ring.AddPoint(xorg + cols * options.cellsize, - yorg - rows * options.cellsize) - ring.AddPoint(xorg, yorg - rows * options.cellsize) + ring.AddPoint(xorg + cols * cellsize, + yorg - rows * cellsize) + ring.AddPoint(xorg, yorg - rows * cellsize) ring.AddPoint(xorg, yorg) ring.CloseRings polygon = ogr.Geometry(ogr.wkbPolygon) @@ -248,8 +268,9 @@ 'With this amount of cells your model will run slow.\nConsider a larger cell-size. Fast models run with < 1,000,000 cells') logger, ch = wt.closeLogger(logger, ch) del logger, ch - sys.exit(1) + #sys.exit(1) if __name__ == "__main__": - main() + argdict = parse_args() + main(**vars(argdict)) \ No newline at end of file Index: wflow-py/Scripts/wtools_py/static_maps.py =================================================================== diff -u -r405114bf53fffa2637a5383483a1fb99fbc73f39 -r97c6554831d691a3f3d1298dfb70c96808dae45e --- wflow-py/Scripts/wtools_py/static_maps.py (.../static_maps.py) (revision 405114bf53fffa2637a5383483a1fb99fbc73f39) +++ wflow-py/Scripts/wtools_py/static_maps.py (.../static_maps.py) (revision 97c6554831d691a3f3d1298dfb70c96808dae45e) @@ -35,17 +35,19 @@ import pcraster as pcr import wflow.wflowtools_lib as wt + +import fiona +from rasterio import features +import rasterio + driver = ogr.GetDriverByName("ESRI Shapefile") def clip_catchment_by_cell(cell_geom, catchment_geom): return cell_geom.intersection(catchment_geom) - -def main(): - - ### Read input arguments ##### - logfilename = 'wtools_static_maps.log' +def parse_args(): + ### Read input arguments ##### parser = OptionParser() usage = "usage: %prog [options]" parser = OptionParser(usage=usage) @@ -57,6 +59,9 @@ parser.add_option('-s', '--source', dest='source', default='wflow', help='Source folder containing clone (default=./wflow)') + parser.add_option('-f', '--forceoutlet', + dest='forceoutlet', default=False, action='store_true', + help='Force the outlet to river cell closest to DEM-edge') parser.add_option('-d', '--destination', dest='destination', default='staticmaps', help='Destination folder (default=./staticmaps)') @@ -90,20 +95,32 @@ parser.add_option('-A', '--alltouch', dest='alltouch', default=False, action='store_true', help='option to burn catchments "all touching".\nUseful when catchment-size is small compared to cellsize') + parser.add_option('-l', '--logfile', + dest='logfilename', default='wtools_static_maps.log', + help='log file name') + (options, args) = parser.parse_args() + return options + +def main(source,destination,inifile,dem_in,rivshp,catchshp,gaugeshp=None,landuse=None,soil=None,lai=None,other_maps=None,logfilename='wtools_static_maps.log',verbose=True,clean=True,alltouch=False,forceoutlet=False): # parse other maps into an array - options.other_maps = options.other_maps.replace( - ' ', '').replace('[', '').replace(']', '').split(',') + if not other_maps == None: + if type(other_maps) == str: + print other_maps + other_maps = other_maps.replace( + ' ', '').replace('[', '').replace(']', '').split(',') + + + + source = os.path.abspath(source) + clone_map = os.path.join(source, 'mask.map') + clone_shp = os.path.join(source, 'mask.shp') + clone_prj = os.path.join(source, 'mask.prj') - options.source = os.path.abspath(options.source) - clone_map = os.path.join(options.source, 'mask.map') - clone_shp = os.path.join(options.source, 'mask.shp') - clone_prj = os.path.join(options.source, 'mask.prj') - - if None in (options.inifile, - options.rivshp, - options.catchshp, - options.dem_in): + if None in (inifile, + rivshp, + catchshp, + dem_in): msg = """The following files are compulsory: - ini file - DEM (raster) @@ -113,31 +130,31 @@ print(msg) parser.print_help() sys.exit(1) - if not os.path.exists(options.inifile): + if not os.path.exists(inifile): print 'path to ini file cannot be found' sys.exit(1) - if not os.path.exists(options.rivshp): + if not os.path.exists(rivshp): print 'path to river shape cannot be found' sys.exit(1) - if not os.path.exists(options.catchshp): + if not os.path.exists(catchshp): print 'path to catchment shape cannot be found' sys.exit(1) - if not os.path.exists(options.dem_in): + if not os.path.exists(dem_in): print 'path to DEM cannot be found' sys.exit(1) # open a logger, dependent on verbose print to screen or not - logger, ch = wt.setlogger(logfilename, 'WTOOLS', options.verbose) + logger, ch = wt.setlogger(logfilename, 'WTOOLS', verbose) # create directories # TODO: check if workdir is still necessary, try to # keep in memory as much as possible # delete old files (when the source and destination folder are different) - if np.logical_and(os.path.isdir(options.destination), - options.destination is not options.source): - shutil.rmtree(options.destination) - if options.destination is not options.source: - os.makedirs(options.destination) + if np.logical_and(os.path.isdir(destination), + destination is not source): + shutil.rmtree(destination) + if destination is not source: + os.makedirs(destination) # Read mask if not(os.path.exists(clone_map)): @@ -160,7 +177,7 @@ # READ CONFIG FILE # open config-file - config = wt.OpenConf(options.inifile) + config = wt.OpenConf(inifile) # read settings snapgaugestoriver = wt.configget(config, 'settings', @@ -192,10 +209,7 @@ minorder = wt.configget(config, 'parameters', 'riverorder_min', 3, datatype='int') - percentiles = np.array( - config.get('parameters', 'statisticmaps', '0, 100').replace( - ' ', '').split(','), dtype='float') - + percentiles = np.array(config.get('parameters', 'statisticmaps', '0, 100').replace(' ', '').split(','), dtype='float') # read the parameters for generating a temporary very high resolution grid if unit_clone == 'degree': cellsize_hr = wt.configget(config, 'parameters', @@ -210,7 +224,7 @@ rows_hr = int((float(ymax) - float(ymin)) / cellsize_hr + 2) hr_trans = (float(xmin), cellsize_hr, float(0), float(ymax), 0, -cellsize_hr) - clone_hr = os.path.join(options.destination, 'clone_highres.tif') + clone_hr = os.path.join(destination, 'clone_highres.tif') # make a highres clone as well! wt.CreateTif(clone_hr, rows_hr, cols_hr, hr_trans, srs, 0) @@ -246,7 +260,7 @@ # read mask location (optional) masklayer = wt.configget( - config, 'mask', 'masklayer', options.catchshp) + config, 'mask', 'masklayer', catchshp) # ???? empty = pcr.ifthen(ones == 0, pcr.scalar(0)) @@ -255,7 +269,7 @@ # in old code) # first add a missing value to dem_in - ds = gdal.Open(options.dem_in, gdal.GA_Update) + ds = gdal.Open(dem_in, gdal.GA_Update) RasterBand = ds.GetRasterBand(1) fill_val = RasterBand.GetNoDataValue() @@ -266,135 +280,106 @@ # reproject to clone map: see http://stackoverflow.com/questions/10454316/how-to-project-and-resample-a-grid-to-match-another-grid-with-gdal-python # resample DEM logger.info('Resampling dem from {:s} to {:s}'.format(os.path.abspath( - options.dem_in), os.path.join(options.destination, dem_map))) - gis.gdal_warp(options.dem_in, clone_map, os.path.join( - options.destination, dem_map), format='PCRaster', gdal_interp=gdalconst.GRA_Average) + dem_in), os.path.join(destination, dem_map))) + gis.gdal_warp(dem_in, clone_map, os.path.join( + destination, dem_map), format='PCRaster', gdal_interp=gdalconst.GRA_Average) # retrieve amount of rows and columns from clone # TODO: make windowstats applicable to source/target with different projections. This does not work yet. # retrieve srs from DEM try: - srs_dem = wt.get_projection(options.dem_in) + srs_dem = wt.get_projection(dem_in) except: logger.warning( 'No projection found in DEM, assuming WGS 1984 lat long') srs_dem = osr.SpatialReference() srs_dem.ImportFromEPSG(4326) clone2dem_transform = osr.CoordinateTransformation(srs, srs_dem) # if srs.ExportToProj4() == srs_dem.ExportToProj4(): - for percentile in percentiles: - if percentile >= 100: - logger.info('computing window maximum') - percentile_dem = os.path.join( - options.destination, 'wflow_dem_max.map') - elif percentile <= 0: - logger.info('computing window minimum') - percentile_dem = os.path.join( - options.destination, 'wflow_dem_min.map') - else: - logger.info( - 'computing window {:d} percentile'.format(int(percentile))) - percentile_dem = os.path.join( - options.destination, 'wflow_dem_{:03d}.map'.format(int(percentile))) + + wt.windowstats(dem_in, len(yax), len(xax),trans, srs,destination,percentiles, transform=clone2dem_transform, logger=logger) - percentile_dem = os.path.join( - options.destination, 'wflow_dem_{:03d}.map'.format(int(percentile))) - stats = wt.windowstats(options.dem_in, len(yax), len(xax), - trans, srs, percentile_dem, percentile, transform=clone2dem_transform, logger=logger) -# else: -# logger.warning('Projections of DEM and clone are different. DEM statistics for different projections is not yet implemented') - - """ - - # burn in rivers - # first convert and clip the river shapefile - # retrieve river shape projection, if not available assume EPSG:4326 - file_att = os.path.splitext(os.path.basename(options.rivshp))[0] - ds = ogr.Open(options.rivshp) - lyr = ds.GetLayerByName(file_att) - extent = lyr.GetExtent() - extent_in = [extent[0], extent[2], extent[1], extent[3]] - try: - # get spatial reference from shapefile - srs_rivshp = lyr.GetSpatialRef() - logger.info('Projection in river shapefile is {:s}'.format(srs_rivshp.ExportToProj4())) - except: - logger.warning('No projection found in {:s}, assuming WGS 1984 lat-lon'.format(options.rivshp)) - srs_rivshp = osr.SpatialReference() - srs_rivshp.ImportFromEPSG(4326) - rivprojshp = os.path.join(options.destination, 'rivshp_proj.shp') - logger.info('Projecting and clipping {:s} to {:s}'.format(options.rivshp, rivprojshp)) - # TODO: Line below takes a very long time to process, the bigger the shapefile, the more time. How do we deal with this? - call(('ogr2ogr','-s_srs', srs_rivshp.ExportToProj4(),'-t_srs', srs.ExportToProj4(), '-clipsrc', '{:f}'.format(xmin), '{:f}'.format(ymin), '{:f}'.format(xmax), '{:f}'.format(ymax), rivprojshp, options.rivshp)) - """ - - # TODO: BURNING!! - - # project catchment layer to projection of clone - file_att = os.path.splitext(os.path.basename(options.catchshp))[0] - print options.catchshp - ds = ogr.Open(options.catchshp) - lyr = ds.GetLayerByName(file_att) - extent = lyr.GetExtent() - extent_in = [extent[0], extent[2], extent[1], extent[3]] - try: - # get spatial reference from shapefile - srs_catchshp = lyr.GetSpatialRef() - logger.info('Projection in catchment shapefile is {:s}'.format( - srs_catchshp.ExportToProj4())) - except: - logger.warning( - 'No projection found in {:s}, assuming WGS 1984 lat-lon'.format(options.catchshp)) - srs_catchshp = osr.SpatialReference() - srs_catchshp.ImportFromEPSG(4326) - catchprojshp = os.path.join(options.destination, 'catchshp_proj.shp') - logger.info('Projecting {:s} to {:s}'.format( - options.catchshp, catchprojshp)) - call(('ogr2ogr', '-s_srs', srs_catchshp.ExportToProj4(), '-t_srs', srs.ExportToProj4(), '-clipsrc', - '{:f}'.format(xmin), '{:f}'.format(ymin), '{:f}'.format(xmax), '{:f}'.format(ymax), catchprojshp, options.catchshp)) - - # + ## read catchment shape-file to create catchment map + src = rasterio.open(clone_map) + shapefile = fiona.open(catchshp,"r") + catchment_shapes = [feature["geometry"] for feature in shapefile] + image = features.rasterize(catchment_shapes,out_shape=src.shape,all_touched=True,transform=src.transform) + catchment_domain = pcr.numpy2pcr(pcr.Ordinal,image.copy(),0) + pcr.report(catchment_domain,catchment_map) + + ## read river shape-file and create burn layer + shapefile = fiona.open(rivshp,"r") + river_shapes = [feature["geometry"] for feature in shapefile] + image = features.rasterize(river_shapes,out_shape=src.shape,all_touched=False,transform=src.transform) + rivers = pcr.numpy2pcr(pcr.Nominal,image.copy(),0) + riverdem = pcr.scalar(rivers) * pcr.readmap(os.path.join(destination, dem_map)) + pcr.setglobaloption("lddin") + riverldd = pcr.lddcreate(riverdem, 1e35, 1e35, 1e35, 1e35) + if forceoutlet: + riverorder = pcr.streamorder(riverldd) + riverorder = pcr.ifthen(pcr.pcrnot(riverorder == pcr.mapmaximum(riverorder)),pcr.nominal(1)) + riverorder = pcr.cover(riverorder,pcr.ifthen(pcr.scalar(riverldd) == 5,pcr.nominal(1))) + riverclump = pcr.clump(riverorder) + riverpitsus = pcr.ifthen(pcr.downstream(riverldd,pcr.scalar(riverldd)) == 5, pcr.scalar(1)) + outsidedem = pcr.ifthen(pcr.pcrnot(pcr.defined(pcr.readmap(os.path.join(destination, dem_map)))),pcr.boolean(1)) + dist2outside = pcr.spread(pcr.cover(outsidedem,0),0,1) + mindistselec = pcr.ifthen(pcr.areaminimum(dist2outside,riverclump) == dist2outside,pcr.scalar(1)) + pitsectionselec = pcr.ifthen(pcr.areamaximum(pcr.cover(riverpitsus,0),riverclump) ==1,pcr.scalar(1)) + riveroutlet = pcr.cover(mindistselec * pitsectionselec * 1000,0) + riverldd = pcr.lddcreate(riverdem - riveroutlet, 1e35, 1e35, 1e35, 1e35) + else: riveroutlet = pcr.cover(pcr.ifthen(pcr.scalar(riverldd) == 5, pcr.scalar(1000)),0) + burn_layer = pcr.cover((pcr.scalar(pcr.ifthen(pcr.streamorder(riverldd) > 1, pcr.streamorder(riverldd) ))-1 )*1000 + riveroutlet,0) + + ## create ldd per catchment logger.info('Calculating ldd') - ldddem = pcr.readmap(os.path.join(options.destination, dem_map)) - ldd_select = pcr.lddcreate(ldddem, 1e35, 1e35, 1e35, 1e35) - pcr.report(ldd_select, os.path.join(options.destination, 'wflow_ldd.map')) + wflow_ldd = pcr.ldd(clone_map) + for idx, shape in enumerate(catchment_shapes): + logger.info('Computing ldd for catchment ' + str(idx) + '/' + str(len(catchment_shapes))) + image = features.rasterize([shape],out_shape=src.shape,all_touched=True,transform=src.transform) + catchment = pcr.numpy2pcr(pcr.Scalar,image.copy(),0) + ldddem = (pcr.readmap(os.path.join(destination, dem_map)) * pcr.scalar(catchment_domain) * catchment) - burn_layer + ldd_select = pcr.lddcreate(ldddem, 1e35, 1e35, 1e35, 1e35) + wflow_ldd=pcr.cover(wflow_ldd,ldd_select) + + pcr.report(wflow_ldd, os.path.join(destination, 'wflow_ldd.map')) + # compute stream order, identify river cells - streamorder = pcr.ordinal(pcr.streamorder(ldd_select)) + streamorder = pcr.ordinal(pcr.streamorder(wflow_ldd)) river = pcr.ifthen(streamorder >= pcr.ordinal(minorder), pcr.boolean(1)) # find the minimum value in the DEM and cover missing values with a river with this value. Effect is none!! so now left out! - # mindem = int(np.min(pcr.pcr2numpy(pcr.ordinal(os.path.join(options.destination, dem_map)),9999999))) - # dem_resample_map = pcr.cover(os.path.join(options.destination, dem_map), pcr.scalar(river)*0+mindem) - # pcr.report(dem_resample_map, os.path.join(options.destination, dem_map)) - pcr.report(streamorder, os.path.join(options.destination, streamorder_map)) - pcr.report(river, os.path.join(options.destination, river_map)) + # mindem = int(np.min(pcr.pcr2numpy(pcr.ordinal(os.path.join(destination, dem_map)),9999999))) + # dem_resample_map = pcr.cover(os.path.join(destination, dem_map), pcr.scalar(river)*0+mindem) + # pcr.report(dem_resample_map, os.path.join(destination, dem_map)) + pcr.report(streamorder, os.path.join(destination, streamorder_map)) + pcr.report(river, os.path.join(destination, river_map)) # deal with your catchments - if options.gaugeshp == None: + if gaugeshp == None: logger.info('No gauges defined, using outlets instead') gauges = pcr.ordinal( pcr.uniqueid( pcr.boolean( - pcr.ifthen(pcr.scalar(ldd_select) == 5, + pcr.ifthen(pcr.scalar(wflow_ldd) == 5, pcr.boolean(1) ) ) ) ) - pcr.report(gauges, os.path.join(options.destination, gauges_map)) + pcr.report(gauges, os.path.join(destination, gauges_map)) # TODO: Add the gauge shape code from StaticMaps.py (line 454-489) # TODO: add river length map (see SticMaps.py, line 492-499) # report river length # make a high resolution empty map - dem_hr_file = os.path.join(options.destination, 'dem_highres.tif') - burn_hr_file = os.path.join(options.destination, 'burn_highres.tif') - demburn_hr_file = os.path.join(options.destination, 'demburn_highres.map') - riv_hr_file = os.path.join(options.destination, 'riv_highres.map') - gis.gdal_warp(options.dem_in, clone_hr, dem_hr_file) + dem_hr_file = os.path.join(destination, 'dem_highres.tif') + burn_hr_file = os.path.join(destination, 'burn_highres.tif') + demburn_hr_file = os.path.join(destination, 'demburn_highres.map') + riv_hr_file = os.path.join(destination, 'riv_highres.map') + gis.gdal_warp(dem_in, clone_hr, dem_hr_file) # wt.CreateTif(riv_hr, rows_hr, cols_hr, hr_trans, srs, 0) - file_att = os.path.splitext(os.path.basename(options.rivshp))[0] + file_att = os.path.splitext(os.path.basename(rivshp))[0] # open the shape layer - ds = ogr.Open(options.rivshp) + ds = ogr.Open(rivshp) lyr = ds.GetLayerByName(file_att) gis.ogr_burn(lyr, clone_hr, -100, file_out=burn_hr_file, format='GTiff', gdal_type=gdal.GDT_Float32, fill_value=0) @@ -409,75 +394,76 @@ xax_hr, yax_hr, demburn_hr, -9999.) pcr.setclone(demburn_hr_file) demburn_hr = pcr.readmap(demburn_hr_file) + + logger.info('Calculating ldd to determine river length') ldd_hr = pcr.lddcreate(demburn_hr, 1e35, 1e35, 1e35, 1e35) - pcr.report(ldd_hr, os.path.join(options.destination, 'ldd_hr.map')) + pcr.report(ldd_hr, os.path.join(destination, 'ldd_hr.map')) pcr.setglobaloption('unitcell') riv_hr = pcr.scalar(pcr.streamorder(ldd_hr) >= minorder) * pcr.downstreamdist(ldd_hr) pcr.report(riv_hr, riv_hr_file) pcr.setglobaloption('unittrue') pcr.setclone(clone_map) logger.info('Computing river length') - #riverlength = wt.windowstats(riv_hr,clone_rows,clone_columns,clone_trans,srs_clone,resultdir,'frac',clone2dem_transform) - riverlength = wt.windowstats(riv_hr_file, len(yax), len(xax), - trans, srs, os.path.join(options.destination, riverlength_fact_map), stat='fact', logger=logger) + wt.windowstats(riv_hr_file, len(yax), len(xax),trans, srs, destination, stat='fact', transform=clone2dem_transform, logger=logger) + # TODO: nothing happends with the river lengths yet. Need to decide how to # use these # report outlet map - pcr.report(pcr.ifthen(pcr.ordinal(ldd_select) == 5, pcr.ordinal(1)), - os.path.join(options.destination, outlet_map)) + pcr.report(pcr.ifthen(pcr.ordinal(wflow_ldd) == 5, pcr.ordinal(1)), + os.path.join(destination, outlet_map)) # report subcatchment map - subcatchment = pcr.subcatchment(ldd_select, gauges) + subcatchment = pcr.subcatchment(wflow_ldd, gauges) pcr.report(pcr.ordinal(subcatchment), os.path.join( - options.destination, subcatch_map)) + destination, subcatch_map)) # Report land use map - if options.landuse == None: + if landuse == None: logger.info('No land use map used. Preparing {:s} with only ones.'. - format(os.path.join(options.destination, landuse_map))) + format(os.path.join(destination, landuse_map))) pcr.report(pcr.nominal(ones), os.path.join( - options.destination, landuse_map)) + destination, landuse_map)) else: logger.info('Resampling land use from {:s} to {:s}'. - format(os.path.abspath(options.landuse), - os.path.join(options.destination, os.path.abspath(landuse_map)))) - gis.gdal_warp(options.landuse, + format(os.path.abspath(landuse), + os.path.join(destination, os.path.abspath(landuse_map)))) + gis.gdal_warp(landuse, clone_map, - os.path.join(options.destination, landuse_map), + os.path.join(destination, landuse_map), format='PCRaster', gdal_interp=gdalconst.GRA_Mode, gdal_type=gdalconst.GDT_Int32) # report soil map - if options.soil == None: + if soil == None: logger.info('No soil map used. Preparing {:s} with only ones.'. - format(os.path.join(options.destination, soil_map))) + format(os.path.join(destination, soil_map))) pcr.report(pcr.nominal(ones), os.path.join( - options.destination, soil_map)) + destination, soil_map)) else: logger.info('Resampling soil from {:s} to {:s}'. - format(os.path.abspath(options.soil), - os.path.join(options.destination, os.path.abspath(soil_map)))) - gis.gdal_warp(options.soil, + format(os.path.abspath(soil), + os.path.join(destination, os.path.abspath(soil_map)))) + gis.gdal_warp(soil, clone_map, - os.path.join(options.destination, soil_map), + os.path.join(destination, soil_map), format='PCRaster', gdal_interp=gdalconst.GRA_Mode, gdal_type=gdalconst.GDT_Int32) - if options.lai == None: + if lai == None: logger.info('No vegetation LAI maps used. Preparing default maps {:s} with only ones.'. - format(os.path.join(options.destination, soil_map))) + format(os.path.join(destination, soil_map))) pcr.report(pcr.nominal(ones), os.path.join( - options.destination, soil_map)) + destination, soil_map)) else: - dest_lai = os.path.join(options.destination, 'clim') + dest_lai = os.path.join(destination, 'clim') os.makedirs(dest_lai) for month in range(12): lai_in = os.path.join( - options.lai, 'LAI00000.{:03d}'.format(month + 1)) + lai, 'LAI00000.{:03d}'.format(month + 1)) lai_out = os.path.join( dest_lai, 'LAI00000.{:03d}'.format(month + 1)) logger.info('Resampling vegetation LAI from {:s} to {:s}'. @@ -489,32 +475,33 @@ format='PCRaster', gdal_interp=gdalconst.GRA_Bilinear, gdal_type=gdalconst.GDT_Float32) - - # report soil map - if options.other_maps == None: +# +# # report soil map + if other_maps == None: logger.info('No other maps used. Skipping other maps.') else: logger.info('Resampling list of other maps...') - for map_file in options.other_maps: + for map_file in other_maps: map_name = os.path.split(map_file)[1] logger.info('Resampling a map from {:s} to {:s}'. format(os.path.abspath(map_file), - os.path.join(options.destination, map_name))) + os.path.join(destination, os.path.splitext(os.path.basename(map_file))[0]+'.map'))) gis.gdal_warp(map_file, clone_map, - os.path.join(options.destination, map_name), + os.path.join(destination,os.path.splitext(os.path.basename(map_file))[0]+'.map'), format='PCRaster', gdal_interp=gdalconst.GRA_Mode, gdal_type=gdalconst.GDT_Float32) - if options.clean: - wt.DeleteList(glob.glob(os.path.join(options.destination, '*.xml')), + if clean: + wt.DeleteList(glob.glob(os.path.join(destination, '*.xml')), logger=logger) - wt.DeleteList(glob.glob(os.path.join(options.destination, 'clim', '*.xml')), + wt.DeleteList(glob.glob(os.path.join(destination, 'clim', '*.xml')), logger=logger) - wt.DeleteList(glob.glob(os.path.join(options.destination, '*highres*')), + wt.DeleteList(glob.glob(os.path.join(destination, '*highres*')), logger=logger) if __name__ == "__main__": - main() + argdict = parse_args() + main(**vars(argdict)) Index: wflow-py/Scripts/wtools_py/wflow_fews.py =================================================================== diff -u --- wflow-py/Scripts/wtools_py/wflow_fews.py (revision 0) +++ wflow-py/Scripts/wtools_py/wflow_fews.py (revision 97c6554831d691a3f3d1298dfb70c96808dae45e) @@ -0,0 +1,100 @@ + +import create_grid as cg +import static_maps as sm +import glob +import os +import pcraster as pcr +import zipfile +from lxml import etree + +# prameters From extent_in: + +source = '.\wflow_catchment\extents\extent_in.xml' +#source = 'd:\FEWS\FEWS_Accelerator\Mekong_SA\Modules\wflow\wflow_catchment\extents\extent_in.xml' + +tree = etree.parse(source) +doc = tree.getroot() +model_name = doc.attrib['locationId'] +rows = float(doc[0].text) +cols = float(doc[1].text) +cellsize = float(doc[4].text) +xmin = float(doc[3][0].text) - 0.5 * cellsize +ymin = float(doc[3][1].text) + 0.5 * cellsize +xmax = xmin + cols * cellsize +ymax = ymin + rows * cellsize + + +#rename dir and change current directory +os.rename(os.path.join(os.getcwd(),'wflow_catchment'),os.path.join(os.getcwd(),model_name)) +os.chdir(os.path.join(os.getcwd(),model_name)) + +#get base-data from GEE + +''' +Martijn and Genna, do your thing!!!! + +ToDo: + - Get the .\data\catchments\catchments.shp (I think GEOJSON is also fine) from GEE + - Get the .\data\rivers\rivers.shp (I think GEOJSON is also fine) from GEE + - Get the .\data\dem\dem.tif from GEE + - Get the .\data\parameters\ directory filled from GEE. *.tif is GeoTiff. LAI00000.* is PCRaster map + .\data\parameters\FirstZoneCapacity.tif + .\data\parameters\FirstZoneKsatVer.tif + .\data\parameters\ FirstZoneMinCapacity.tif + .\data\parameters\InfiltCapSoil.tif + .\data\parameters\ M.tif + .\data\parameters\PathFrac.tif + .\data\parameters\thetaS.tif + .\data\parameters\WaterFrac.tif + .\data\parameters\wflow_landuse.tif + .\data\parameters\wflow_soil.tif + .\data\parameters\clim\LAI00000.001 + .\data\parameters\clim\LAI00000.002 + .\data\parameters\clim\LAI00000.003 + .\data\parameters\clim\LAI00000.004 + .\data\parameters\clim\LAI00000.005 + .\data\parameters\clim\LAI00000.006 + .\data\parameters\clim\LAI00000.007 + .\data\parameters\clim\LAI00000.008 + .\data\parameters\clim\LAI00000.009 + .\data\parameters\clim\LAI00000.010 + .\data\parameters\clim\LAI00000.011 + .\data\parameters\clim\LAI00000.012 +''' + +#create grid +logfilename = 'wtools_create_grid.log' +destination = '.\mask' +inputfile = '.\data\catchments\catchments.shp' +projection = 'EPSG:4326' + +cg.main(logfilename,destination,inputfile,projection,cellsize,snap=True,locationid=model_name) + + +# create static maps +source = r'.\mask' +destination = r'.\staticmaps' +inifile = 'staticmaps.ini' +dem_in = r'.\data\dem\dem.tif' +rivshp = r'.\data\rivers\rivers.shp' +catchshp = r'.\data\catchments\catchments.shp' +lai = r'.\data\parameters\clim' +other_maps=glob.glob(os.path.join(os.getcwd(),'data\parameters','*.tif')) + +sm.main(source,destination,inifile,dem_in,rivshp,catchshp,lai=lai,other_maps=other_maps) + +# save default state-files in FEWS-config +state_dir = r'.\outstate' +state_files = ['CanopyStorage.map','SatWaterDepth.map','Snow.map','SnowWater.map','SurfaceRunoff.map','TSoil.map','UStoreLayerDepth_0.map','WaterLevel.map'] +zip_loc = r'..\\..\\..\\Config\\ColdStateFiles\\' + model_name[:2].upper() + model_name[2:] + '_GA_Historical default.zip' + +mask = pcr.readmap(os.path.join(os.getcwd(),source,'mask.map')) + +zf = zipfile.ZipFile(zip_loc, mode='w') + +for state_file in state_files: + state_path = os.path.join(os.getcwd(),state_dir,state_file) + pcr.report(pcr.cover(mask,pcr.scalar(0)),state_path) + zf.write(state_path, state_file, compress_type=zipfile.ZIP_DEFLATED) + +zf.close() \ No newline at end of file Index: wflow-py/_version.py =================================================================== diff -u -ra70af6fb5d306fffe3f818867365fe27ea19e27e -r97c6554831d691a3f3d1298dfb70c96808dae45e --- wflow-py/_version.py (.../_version.py) (revision a70af6fb5d306fffe3f818867365fe27ea19e27e) +++ wflow-py/_version.py (.../_version.py) (revision 97c6554831d691a3f3d1298dfb70c96808dae45e) @@ -2,4 +2,4 @@ VERSION='1.0.master.1' MVERSION='1.0.master' NVERSION='1.0.1' -BUILD='2017-07-20 11:08:33.863000' +BUILD='2017-07-20 13:07:46.763000' Index: wflow-py/make_wflow_exe_cx.py =================================================================== diff -u -ra70af6fb5d306fffe3f818867365fe27ea19e27e -r97c6554831d691a3f3d1298dfb70c96808dae45e --- wflow-py/make_wflow_exe_cx.py (.../make_wflow_exe_cx.py) (revision a70af6fb5d306fffe3f818867365fe27ea19e27e) +++ wflow-py/make_wflow_exe_cx.py (.../make_wflow_exe_cx.py) (revision 97c6554831d691a3f3d1298dfb70c96808dae45e) @@ -155,7 +155,7 @@ packages.append('bmi') #packages.append('pkg_resources') else: - includes = ['wflow.wflow_bmi', 'wflow.wflow_w3ra', 'wflow.wflow_bmi_combined','lxml.etree', 'lxml._elementpath', 'gzip', 'numpy.core._methods', 'numpy.lib.format'] + includes = ["matplotlib.backends.backend_qt4agg", 'wflow.wflow_bmi', 'wflow.wflow_w3ra', 'wflow.wflow_bmi_combined','lxml.etree', 'lxml._elementpath', 'gzip', 'numpy.core._methods', 'numpy.lib.format'] # "include_msvcr": True, options = {"includes": includes, "packages": packages,'include_files': data_files, "build_exe": thename, Index: wflow-py/wflow/__init__.py =================================================================== diff -u -ra70af6fb5d306fffe3f818867365fe27ea19e27e -r97c6554831d691a3f3d1298dfb70c96808dae45e --- wflow-py/wflow/__init__.py (.../__init__.py) (revision a70af6fb5d306fffe3f818867365fe27ea19e27e) +++ wflow-py/wflow/__init__.py (.../__init__.py) (revision 97c6554831d691a3f3d1298dfb70c96808dae45e) @@ -2,7 +2,7 @@ __version__='1.0.master' __release__='1.0.master.1' __versionnr__='1.0.1' -__build__='2017-07-20 11:08:33.863000' +__build__='2017-07-20 13:07:46.763000' import os, sys import osgeo.gdal as gdal Index: wflow-py/wflow/wf_DynamicFramework.py =================================================================== diff -u -ra70af6fb5d306fffe3f818867365fe27ea19e27e -r97c6554831d691a3f3d1298dfb70c96808dae45e --- wflow-py/wflow/wf_DynamicFramework.py (.../wf_DynamicFramework.py) (revision a70af6fb5d306fffe3f818867365fe27ea19e27e) +++ wflow-py/wflow/wf_DynamicFramework.py (.../wf_DynamicFramework.py) (revision 97c6554831d691a3f3d1298dfb70c96808dae45e) @@ -564,7 +564,6 @@ self.logger.error("Variable change string (apply_once) could not be executed: " + execstr) - @profile def wf_updateparameters(self): """ Update the model Parameters (can be used in static and dynamic part of the model) @@ -2421,7 +2420,7 @@ - @profile + def wf_readmap(self, name, default, verbose=True,fail=False,ncfilesource="not set",silent=False): """ Adjusted version of readmapNew. the style variable is used to indicated Index: wflow-py/wflow/wflow_sbm.py =================================================================== diff -u -ra70af6fb5d306fffe3f818867365fe27ea19e27e -r97c6554831d691a3f3d1298dfb70c96808dae45e --- wflow-py/wflow/wflow_sbm.py (.../wflow_sbm.py) (revision a70af6fb5d306fffe3f818867365fe27ea19e27e) +++ wflow-py/wflow/wflow_sbm.py (.../wflow_sbm.py) (revision 97c6554831d691a3f3d1298dfb70c96808dae45e) @@ -1041,8 +1041,6 @@ self.GWScale = (self.DemMax - self.DrainageBase) / self.SoilThickness / self.RunoffGeneratingGWPerc - - @profile def dynamic(self): """ Stuf that is done for each timestep of the model