Index: examples/wflow_rhine_sbm/wflow_sbm_NC.ini =================================================================== diff -u -ra102d31a70edb76290d7a7ec392483e23c948286 -rced522d9cb1398cc92a7bbbcc69d23b88708d088 --- examples/wflow_rhine_sbm/wflow_sbm_NC.ini (.../wflow_sbm_NC.ini) (revision a102d31a70edb76290d7a7ec392483e23c948286) +++ examples/wflow_rhine_sbm/wflow_sbm_NC.ini (.../wflow_sbm_NC.ini) (revision ced522d9cb1398cc92a7bbbcc69d23b88708d088) @@ -1,9 +1,3 @@ -[API] -IF=0,1 -InwaterMM=2,4 -#FirstZoneDepth=2,1 -#Latitude=1,5 -#Longitude=1,5 # Define the forcings needed for the model here @@ -36,7 +30,7 @@ # - monthlyclim: read a map corresponding to the current month (12 maps in total) # - dailyclim: read a map corresponding to the current day of the year # - hourlyclim: read a map corresponding to the current hour of the day (24 in total) (not implemented yet) -RootingDepth=monthlyclim/ROOT,monthlyclim,75,1 +#RootingDepth=monthlyclim/ROOT,monthlyclim,75,1 [run] # either a runinfo file or a start and end-time are required @@ -78,8 +72,11 @@ # netcdfoutput requires also outputformat = 1 (default) and additionally the name of the file netcdfoutput = outmaps.nc +netcdfstatesoutput = states.nc netcdfinput= inmaps.nc +netcdfstatesinput = instates.nc netcdfwritebuffer=100 +netcdf_least_significant_digit=2 @@ -91,7 +88,7 @@ [outputmaps] -#self.Inflow=iflow +self.Inflow=iflow #self.WaterDem=wat #self.SubCellFrac=scf Index: wflow-py/Scripts/pcr2netcdf.py =================================================================== diff -u -r0171e35884ddb86a9eb8614d8162999700954e38 -rced522d9cb1398cc92a7bbbcc69d23b88708d088 --- wflow-py/Scripts/pcr2netcdf.py (.../pcr2netcdf.py) (revision 0171e35884ddb86a9eb8614d8162999700954e38) +++ wflow-py/Scripts/pcr2netcdf.py (.../pcr2netcdf.py) (revision ced522d9cb1398cc92a7bbbcc69d23b88708d088) @@ -26,7 +26,7 @@ -S startdate in "%d-%m-%Y %H:%M:%S" e.g. 31-12-1990 00:00:00 -E endDate in "%d-%m-%Y %H:%M:%S" -s startstep (in the mapstack, default = 1) - -Y Do not make seperate files per year + -Y Make seperate files per year -P set the EPSG string. default: "EPSG:4326" -N Mapstack-name (prefix) You can sepcify multiple input mapstack to merge them into one netcdf @@ -320,17 +320,17 @@ print "ERROR: Failed to initialize logger with logfile: " + logfilename sys.exit(2) -def date_range(start, end, tdelta="days"): +# def date_range(start, end, tdelta="days"): +# +# +# +# if tdelta == "days": +# r = (end+dt.timedelta(days=1)-start).days +# return [start+dt.timedelta(days=i) for i in range(r)] +# else: +# r = (end+dt.timedelta(days=1)-start).days * 24 +# return [start+dt.timedelta(hours=i) for i in range(r)] - - - if tdelta == "days": - r = (end+dt.timedelta(days=1)-start).days - return [start+dt.timedelta(days=i) for i in range(r)] - else: - r = (end+dt.timedelta(days=1)-start).days * 24 - return [start+dt.timedelta(hours=i) for i in range(r)] - def date_range_peryear(start, end, tdelta="days"): ret = [] @@ -346,6 +346,11 @@ return ret +def date_range(start, end, timestepsecs): + r = int((end + dt.timedelta(seconds=timestepsecs) - start).total_seconds()/timestepsecs) + return [start + dt.timedelta(seconds=(timestepsecs * i)) for i in range(r)] + + def main(argv=None): """ Perform command line execution of the model. @@ -377,7 +382,7 @@ zlib=True least_significant_digit=None startstep = 1 - perYear = True + perYear = False if argv is None: argv = sys.argv[1:] if len(argv) == 0: @@ -393,13 +398,14 @@ for o, a in opts: if o == '-S': startstr = a - if o == '-s': startstep = int(a) + if o == '-s': + startstep = int(a) if o == '-E': endstr = a if o == '-O': ncoutfile = a if o == '-c': inifile = a if o == '-I': mapstackfolder = a if o == '-b': mbuf = int(a) - if o == '-Y': perYear = False + if o == '-Y': perYear = True if o == '-z': zlib=True if o == '-P': EPSG = a if o == '-F': Format=a @@ -432,12 +438,12 @@ if perYear: timeList = date_range_peryear(start, end, tdelta="days") else: - timeList = date_range(start, end, tdelta="days") + timeList = date_range(start, end, timestepsecs) else: if perYear: timeList = date_range_peryear(start, end, tdelta="hours") else: - timeList = date_range(start, end, tdelta="hours") + timeList = date_range(start, end,timestepsecs) if os.path.exists(inifile): inimetadata = getnetcdfmetafromini(inifile) Index: wflow-py/wflow/wf_DynamicFramework.py =================================================================== diff -u -r8ca1b4bc61b3bc0d9c58451a4685378b663c5034 -rced522d9cb1398cc92a7bbbcc69d23b88708d088 --- wflow-py/wflow/wf_DynamicFramework.py (.../wf_DynamicFramework.py) (revision 8ca1b4bc61b3bc0d9c58451a4685378b663c5034) +++ wflow-py/wflow/wf_DynamicFramework.py (.../wf_DynamicFramework.py) (revision ced522d9cb1398cc92a7bbbcc69d23b88708d088) @@ -615,9 +615,11 @@ self.outputFormat = int(configget(self._userModel().config, 'framework', 'outputformat', '1')) self.APIDebug = int(configget(self._userModel().config, 'framework', 'debug', str(self.APIDebug))) self.ncfile = configget(self._userModel().config, 'framework', 'netcdfinput', "None") + self.ncfilestates = configget(self._userModel().config, 'framework', "netcdfstatesinput", "None") self.ncoutfile = configget(self._userModel().config, 'framework', 'netcdfoutput', "None") self.ncoutfilestatic = configget(self._userModel().config, 'framework', 'netcdfstaticoutput', "None") - self.ncinfilestatic = configget(self._userModel().config, 'framework', 'netcdfstaticinput', "None") + self.ncoutfilestate = configget(self._userModel().config, 'framework', 'netcdfstatesoutput', "None") + self.ncfilestatic = configget(self._userModel().config, 'framework', 'netcdfstaticinput', "None") self.EPSG = configget(self._userModel().config, 'framework', 'EPSG', "EPSG:4326") self.ncfileformat = configget(self._userModel().config, 'framework', 'netcdf_format', "NETCDF4") self.ncfilecompression = configget(self._userModel().config, 'framework', 'netcdf_zlib', "True") @@ -680,7 +682,7 @@ self._userModel()._setFirstTimeStep(self._d_firstTimestep) self._d_lastTimestep = int(nrseconds / self.timestepsecs) + self._d_firstTimestep - + # Setup all the netCDF files that may be used for input/output if self.ncfile != "None": mstacks = configsection(self._userModel().config, "inputmapstacks") varlst = [] @@ -689,6 +691,13 @@ self.logger.debug("Found following input variables to get from netcdf file: " + str(varlst)) self.NcInput = netcdfinput(os.path.join(caseName, self.ncfile), self.logger, varlst) + # Setup all the netCDF files that may be used for input/output + if self.ncfilestates != "None": + smaps = self._userModel().stateVariables() + maps = [s + ".map" for s in smaps] + self.logger.debug("Found following input states to get from netcdf file: " + str(maps)) + self.NcInputStates = netcdfinputstates(os.path.join(caseName, self.ncfilestates), self.logger, maps) + if self.ncoutfile != 'None': # Ncoutput buffer = int(configget(self._userModel().config, 'framework', 'netcdfwritebuffer', "50")) meta = {} @@ -714,14 +723,25 @@ maxbuf=1, metadata=meta, EPSG=self.EPSG,Format=self.ncfileformat, zlib=self.ncfilecompression,least_significant_digit=self.ncfiledigits) + if self.ncoutfilestate != 'None': # Ncoutput + meta = {} + meta['caseName'] = caseName + meta['runId'] = runId + meta['wflow_version'] =__version__ + meta['wflow_release'] =__release__ + self.NcOutputState = netcdfoutputstatic(os.path.join(caseName, runId, self.ncoutfilestate), + self.logger, self.datetime_laststep,1,timestepsecs=self.timestepsecs, + maxbuf=1, metadata=meta, EPSG=self.EPSG,Format=self.ncfileformat, + zlib=self.ncfilecompression,least_significant_digit=self.ncfiledigits) - if self.ncinfilestatic != 'None': # Ncoutput + + if self.ncfilestatic != 'None': # Ncoutput meta = {} meta['caseName'] = caseName meta['runId'] = runId meta['wflow_version'] =__version__ meta['wflow_release'] =__release__ - self.NcInputStatic = netcdfinput(os.path.join(caseName, self.ncinfilestatic), self.logger, varlst) + self.NcInputStatic = netcdfinput(os.path.join(caseName, self.ncfilestatic), self.logger, varlst) # Fill the summary (stat) list from the ini file self.statslst = [] @@ -842,13 +862,13 @@ for z in savevar: fname = os.path.join(directory, var + "_" + str(a)).replace("\\", "/") + ".map" # report(z,fname) - self.reportStatic(z, fname, style=1, gzipit=False, longname=fname) + self.reportState(z, fname, style=1, gzipit=False, longname=fname) a = a + 1 except: # execstr = "report(self._userModel()." + var +",\"" + fname + "\")" # exec execstr thevar = eval("self._userModel()." + var) - self.reportStatic(thevar, fname, style=1, gzipit=False, longname=fname) + self.reportState(thevar, fname, style=1, gzipit=False, longname=fname) except: self.logger.warn("Problem saving state variable: " + var) self.logger.warn(execstr) @@ -966,8 +986,11 @@ if nr == 0: try: mpath = os.path.join(directory, var + ".map").replace("\\", "/") - execstr = "self._userModel()." + var + "= readmap(\"" + mpath + "\")" - exec execstr + tvar = self.wf_readmap(mpath,0.0,ncfilesource=self.ncfilestates) + wf_readmtvar = self.wf_readmap(mpath,0.0,ncfilesource=self.ncfilestates) + setattr(self._userModel(), var,tvar) + #execstr = "self._userModel()." + var + "= readmap(\"" + mpath + "\")" + #exec execstr except: self.logger.error( "problem while reading state variable from disk: " + mpath + " Suggest to use the -I uption to restart") @@ -1761,6 +1784,44 @@ elif self.outputFormat == 4: numpy.savetxt(path, pcr2numpy(variable, -999), fmt="%0.6g") + + def reportState(self, variable, name, style=1, gzipit=False, longname=None): + """ + + :param variable: + :param name: + :param style: + :param gzipit: + :param longname: + :return: + """ + if longname == None: + longname = name + path = name + + if self.outputFormat == 1: + if sys.version_info[0] == 2 and sys.version_info[1] >= 6: + try: + import pcraster as PCRaster + except: + import PCRaster as PCRaster + else: + import PCRaster + if not hasattr(self, 'NcOutputState'): + PCRaster.report(variable, path) + if gzipit: + Gzip(path, storePath=True) + else: + self.NcOutputState.savetimestep(1, variable, var=name, name=longname) + + elif self.outputFormat == 2: + numpy.savez(path, pcr2numpy(variable, -999)) + elif self.outputFormat == 3: + numpy.savez(path, pcr2numpy(variable, -999)) + elif self.outputFormat == 4: + numpy.savetxt(path, pcr2numpy(variable, -999), fmt="%0.6g") + + def _reportNew(self, variable, name, style=1, gzipit=False, longname=None): """ outputformat: (set in the [framework] section of the init file). @@ -1876,7 +1937,7 @@ - def wf_readmap(self, name, default, verbose=True,fail=False): + def wf_readmap(self, name, default, verbose=True,fail=False,ncfilesource="not set"): """ Adjusted version of readmapNew. the style variable is used to indicated how the data is read:: @@ -1915,8 +1976,10 @@ newName = name + nameSuffix if self._inResume(): - timestep = self._userModel().firstTimeStep() - newName = generateNameT(name, timestep - 1) + if os.path.splitext(name)[1] == ".map": + newName = name + else: + newName = name + nameSuffix if hasattr(self._userModel(), "_inDynamic"): if self._userModel()._inDynamic() or self._inUpdateWeight(): @@ -1929,8 +1992,11 @@ if self._userModel()._inDynamic(): if self.ncfile != "None": - retval = self.NcInput.gettimestep(self._userModel().currentTimeStep(), self.logger, var=varname) - return retval + retval, succ = self.NcInput.gettimestep(self._userModel().currentTimeStep(), self.logger, var=varname) + if succ: + return retval + else: + return scalar(default) if os.path.isfile(path): mapje = readmap(path) @@ -1943,6 +2009,13 @@ sys.exit(1) return scalar(default) if self._userModel()._inInitial(): + if ncfilesource == self.ncfilestatic: + retval, succ = self.NcInputStatic.gettimestep(1, self.logger, var=varname) + if succ: + return retval + else: + return scalar(default) + if os.path.isfile(path): mapje = readmap(path) return mapje @@ -1954,7 +2027,26 @@ sys.exit(1) return scalar(default) + if self._inResume(): + if ncfilesource == self.ncfilestates and ncfilesource not in 'None': + retval, succ = self.NcInputStates.gettimestep(1, self.logger, var=varname) + if succ: + return retval + else: + return scalar(default) + if os.path.isfile(path): + mapje = readmap(path) + return mapje + else: + if verbose: + self.logger.debug("State input data (" + os.path.abspath(path) + ") not present, returning " + str(default)) + if fail: + self.logger.error("Required map: " + os.path.abspath(path) + " not found, exiting..") + sys.exit(1) + return scalar(default) + + elif style == 2: # Assuming they are set in memory by the API # first get basename (last bit of path) name = os.path.basename(name) Index: wflow-py/wflow/wf_netcdfio.py =================================================================== diff -u -r9d6d3f02533e3452466d6115f29c905c9eb646f2 -rced522d9cb1398cc92a7bbbcc69d23b88708d088 --- wflow-py/wflow/wf_netcdfio.py (.../wf_netcdfio.py) (revision 9d6d3f02533e3452466d6115f29c905c9eb646f2) +++ wflow-py/wflow/wf_netcdfio.py (.../wf_netcdfio.py) (revision ced522d9cb1398cc92a7bbbcc69d23b88708d088) @@ -458,10 +458,71 @@ logging.debug("reading new netcdf data block starting at: " + str(ncindex)) for vars in self.alldat: self.alldat[vars] = self.dataset.variables[vars][ncindex:ncindex + self.maxsteps] + self.fstep = ncindex self.lstep = ncindex + self.maxsteps np_step = self.alldat[var][ncindex - self.fstep, :, :] - return numpy2pcr(Scalar, np_step, 1E31) + miss = float(self.dataset.variables[var]._FillValue) + return numpy2pcr(Scalar, np_step, miss), True else: + #logging.debug("Var (" + var + ") not found returning 0") + return cover(scalar(0.0)), False + + +class netcdfinputstates(): + def __init__(self, netcdffile, logging, vars=[]): + """ + First try to setup a class read netcdf files + (converted with pcr2netcdf.py) + + netcdffile: file to read the forcing data from + logging: python logging object + vars: list of variables to get from file + """ + + if os.path.exists(netcdffile): + self.dataset = netCDF4.Dataset(netcdffile, mode='r') + else: + logging.error(os.path.abspath(netcdffile) + " not found!") + exit(ValueError) + + logging.info("Reading state input from netCDF file: " + netcdffile + ": " + str(self.dataset).replace('\n', ' ')) + self.alldat = {} + a = pcr2numpy(cover(0.0), 0.0).flatten() + # Determine steps to load in mem based on estimated memory usage + floatspermb = 1048576 / 4 + maxmb = 4000 + self.maxsteps = maxmb * len(a) / floatspermb + 1 + self.fstep = 0 + self.lstep = self.fstep + self.maxsteps + + + for var in vars: + try: + self.alldat[var] = self.dataset.variables[var][self.fstep:self.maxsteps] + except: + self.alldat.pop(var, None) + logging.warn("Variable " + var + " not found in netcdf file: " + netcdffile) + + def gettimestep(self, timestep, logging, var='P'): + """ + Gets a map for a single timestep. reads data in blocks assuming sequential access + + timestep: framework timestep (1-based) + logging: python logging object + var: variable to get from the file + """ + ncindex = timestep - 1 + if self.alldat.has_key(var): + if ncindex == self.lstep: # Read new block of data in mem + logging.debug("reading new netcdf data block starting at: " + str(ncindex)) + for vars in self.alldat: + self.alldat[vars] = self.dataset.variables[vars][ncindex:ncindex + self.maxsteps] + self.fstep = ncindex + self.lstep = ncindex + self.maxsteps + np_step = self.alldat[var][ncindex - self.fstep, :, :] + miss = float(self.dataset.variables[var]._FillValue) + return numpy2pcr(Scalar, np_step, miss), True + else: logging.debug("Var (" + var + ") not found returning 0") - return cover(scalar(0.0)) + return cover(scalar(0.0)), False Index: wflow-py/wflow/wflow_sbm.py =================================================================== diff -u -rbdd2b962a90dd8a3b47f0a47506c29092438c0c8 -rced522d9cb1398cc92a7bbbcc69d23b88708d088 --- wflow-py/wflow/wflow_sbm.py (.../wflow_sbm.py) (revision bdd2b962a90dd8a3b47f0a47506c29092438c0c8) +++ wflow-py/wflow/wflow_sbm.py (.../wflow_sbm.py) (revision ced522d9cb1398cc92a7bbbcc69d23b88708d088) @@ -1203,7 +1203,7 @@ self.InfiltExcessCubic = self.InfiltExcess * self.ToCubic self.ReinfiltCubic = -1.0 * self.reinfiltwater * self.ToCubic self.Inwater = self.Inwater + self.Inflow # Add abstractions/inflows in m^3/sec - + #self.Inwater = self.Inwater + self.Inflow # Add abstractions/inflows in m^3/sec ########################################################################## # Runoff calculation via Kinematic wave ################################## ##########################################################################