Index: wflow-py/UnitTests/Testwflow_sbm_nc.py =================================================================== diff -u -r048e9fa4ca4b2c6a7c14169de0b192ea2ff9b983 -re71ad2de03777dca02510a47f65dd6f7075a876b --- wflow-py/UnitTests/Testwflow_sbm_nc.py (.../Testwflow_sbm_nc.py) (revision 048e9fa4ca4b2c6a7c14169de0b192ea2ff9b983) +++ wflow-py/UnitTests/Testwflow_sbm_nc.py (.../Testwflow_sbm_nc.py) (revision e71ad2de03777dca02510a47f65dd6f7075a876b) @@ -22,7 +22,7 @@ myModel = wf.WflowModel(wflow_cloneMap, caseName,runId,configfile) # initialise the framework - dynModelFw = wf.wf_DynamicFramework(myModel, stopTime,startTime) + dynModelFw = wf.wf_DynamicFramework(myModel, 0,0) # Load model config from files and check directory structure dynModelFw.createRunId(NoOverWrite=False,level=wf.logging.DEBUG) @@ -55,7 +55,7 @@ self.assertAlmostEquals(3.574188167654313e-09,my_data[:,2].sum(),places=4) my_data = wf.genfromtxt(os.path.join(caseName,runId,"wbsoil.csv"), delimiter=',') print("Checking soil water budget ....") - self.assertAlmostEquals(0.0006305182757557759,my_data[:,2].sum(),places=4) + self.assertAlmostEquals(0.002721056049267645,my_data[:,2].sum(),places=4) print("Checking precip sum ....") my_data = wf.genfromtxt(os.path.join(caseName,runId,"P.csv"), delimiter=',') self.assertAlmostEquals(sump,my_data[:,2].sum()) Index: wflow-py/UnitTests/wflow_sbm/instates.nc =================================================================== diff -u Binary files differ Index: wflow-py/UnitTests/wflow_sbm/wflow_sbm_nc.ini =================================================================== diff -u -r56eaa92ab85d71e87bc31ff0f3015dc890deb253 -re71ad2de03777dca02510a47f65dd6f7075a876b --- wflow-py/UnitTests/wflow_sbm/wflow_sbm_nc.ini (.../wflow_sbm_nc.ini) (revision 56eaa92ab85d71e87bc31ff0f3015dc890deb253) +++ wflow-py/UnitTests/wflow_sbm/wflow_sbm_nc.ini (.../wflow_sbm_nc.ini) (revision e71ad2de03777dca02510a47f65dd6f7075a876b) @@ -44,8 +44,11 @@ # required, base timestep of the model timestepsecs = 3600 #start model with cold state -reinit=1 -runlengthdetermination=intervals +reinit=0 +#runlengthdetermination=intervals +runlengthdetermination=steps +starttime=1990-01-30 01:00:00 +endtime=1990-01-31 06:00:00 # Model parameters and settings [model] @@ -63,9 +66,9 @@ UpdMaxDist=300000.0 #SubCatchFlowOnly = 1 origTopogLateral = 1 -reinit=0 + #DynamicVegetation=1 [misc] @@ -80,8 +83,8 @@ # netcdfoutput requires also outputformat = 1 (default) and additionally the name of the file netcdfoutput = outmaps.nc netcdfstaticoutput = staticoutmaps.nc -netcdfstateoutput = states.nc -#netcdfstatesinput = instates.nc +netcdfstatesoutput = states.nc +netcdfstatesinput = instates.nc netcdfwritebuffer=100 [layout] Index: wflow-py/wflow/wf_DynamicFramework.py =================================================================== diff -u -r9e189ac4f4cf4cc959992318e9121098809a3573 -re71ad2de03777dca02510a47f65dd6f7075a876b --- wflow-py/wflow/wf_DynamicFramework.py (.../wf_DynamicFramework.py) (revision 9e189ac4f4cf4cc959992318e9121098809a3573) +++ wflow-py/wflow/wf_DynamicFramework.py (.../wf_DynamicFramework.py) (revision e71ad2de03777dca02510a47f65dd6f7075a876b) @@ -1065,10 +1065,10 @@ 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.ncinfilestates = 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.ncoutfilestate = configget(self._userModel().config, 'framework', 'netcdfstatesoutput', "None") + self.ncoutfilestates = 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") @@ -1131,11 +1131,11 @@ meta.update(metafrom_config) - if self.ncfilestates != "None": + if self.ncinfilestates != "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) + self.NcInputStates = netcdfinputstates(os.path.join(caseName, self.ncinfilestates), self.logger, maps) if self.ncfilestatic != "None": @@ -1157,8 +1157,8 @@ maxbuf=1, metadata=meta, EPSG=self.EPSG,Format=self.ncfileformat, zlib=self.ncfilecompression,least_significant_digit=self.ncfiledigits) - if self.ncoutfilestate != 'None': # Ncoutput - self.NcOutputState = netcdfoutputstatic(os.path.join(caseName, runId, self.ncoutfilestate), + if self.ncoutfilestates != 'None': # Ncoutput + self.NcOutputState = netcdfoutputstatic(os.path.join(caseName, runId, self.ncoutfilestates), self.logger, self.DT.runEndTime,1,timestepsecs=self.DT.timeStepSecs, maxbuf=1, metadata=meta, EPSG=self.EPSG,Format=self.ncfileformat, zlib=self.ncfilecompression,least_significant_digit=self.ncfiledigits) @@ -1425,19 +1425,29 @@ while stop == 0: name = os.path.join(directory, var + "_" + str(nr) + ".map").replace("\\", "/") - if os.path.exists(name): + try: + tvar = self.wf_readmap(name, 0.0, ncfilesource=self.ncinfilestates,fail=True,silent=True) if nr == 0: exec "self._userModel()." + var + "= []" - execstr = "self._userModel()." + var + ".append(readmap(\"" + name + "\"))" + execstr = "self._userModel()." + var + ".append(tvar)" exec execstr nr = nr + 1 - else: + except: stop = 1 + + #if os.path.exists(name): + # if nr == 0: + # exec "self._userModel()." + var + "= []" + # execstr = "self._userModel()." + var + ".append(readmap(\"" + name + "\"))" + # exec execstr + # nr = nr + 1 + #else: + # stop = 1 if nr == 0: try: mpath = os.path.join(directory, var + ".map").replace("\\", "/") - tvar = self.wf_readmap(mpath,0.0,ncfilesource=self.ncfilestates) - wf_readmtvar = self.wf_readmap(mpath,0.0,ncfilesource=self.ncfilestates,fail=True) + tvar = self.wf_readmap(mpath, 0.0, ncfilesource=self.ncinfilestates) + #wf_readmtvar = self.wf_readmap(mpath,0.0,ncfilesource=self.ncinfilestates,fail=True) setattr(self._userModel(), var,tvar) except: self.logger.error( @@ -2413,7 +2423,7 @@ - def wf_readmap(self, name, default, verbose=True,fail=False,ncfilesource="not set"): + 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 how the data is read:: @@ -2500,7 +2510,7 @@ self.logger.debug("Input data (" + os.path.abspath(path) + ") for timestep not present, returning " + str(default)) if fail: self.logger.error("Required map: " + os.path.abspath(path) + " not found, exiting..") - sys.exit(1) + raise ValueError('Input map not found') return cover(scalar(default)) elif self._userModel()._inInitial(): @@ -2510,8 +2520,9 @@ return retval else: if fail: - self.logger.error("Required map: " + os.path.abspath(path) + " not found in " + self.ncfilestatic + " exiting..") - sys.exit(1) + if not silent: + self.logger.error("Required map: " + os.path.abspath(path) + " not found in " + self.ncfilestatic + " exiting..") + raise ValueError('Input static variable not found in netcdf') else: return self.TheClone + default @@ -2522,16 +2533,22 @@ if verbose: self.logger.debug("Static 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) + if not silent: + self.logger.error("Required map: " + os.path.abspath(path) + " not found, exiting..") + raise ValueError('Input static variable not found') return self.TheClone + default elif self._inResume(): - if ncfilesource == self.ncfilestates and ncfilesource not in 'None': - retval, succ = self.NcInputStates.gettimestep(1, self.logger, var=varname) + if ncfilesource == self.ncinfilestates and ncfilesource not in 'None': + retval, succ = self.NcInputStates.gettimestep(1, self.logger, var=varname,tsdatetime=self.DT.runStateTime) if succ: return retval else: + if fail: + if not silent: + self.logger.error("Required map: " + os.path.abspath(path) + " not found, exiting..") + raise ValueError('Input state variable not found') + return self.TheClone + default if os.path.isfile(path): @@ -2541,8 +2558,9 @@ 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) + if not silent: + self.logger.error("Required map: " + os.path.abspath(path) + " not found, exiting..") + raise ValueError('Input state variable not found') return cover(scalar(default)) else: # Assuming we are in pre-or post loop within the framwork if "None" not in self.ncfilestatic: @@ -2551,8 +2569,9 @@ return retval else: if fail: - self.logger.error("Required map: " + os.path.abspath(path) + " not found in " + self.ncfilestatic + " exiting..") - sys.exit(1) + if not silent: + self.logger.error("Required map: " + os.path.abspath(path) + " not found in " + self.ncfilestatic + " exiting..") + raise ValueError('Input variable not found in netcdf') else: return self.TheClone + default @@ -2563,8 +2582,9 @@ if verbose: self.logger.debug("Static 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) + if not silent: + self.logger.error("Required map: " + os.path.abspath(path) + " not found, exiting..") + raise ValueError('Input variable not found') return self.TheClone + default Index: wflow-py/wflow/wf_netcdfio.py =================================================================== diff -u -rfb4c96a491f37f94570d811301e74bd91819d8a3 -re71ad2de03777dca02510a47f65dd6f7075a876b --- wflow-py/wflow/wf_netcdfio.py (.../wf_netcdfio.py) (revision fb4c96a491f37f94570d811301e74bd91819d8a3) +++ wflow-py/wflow/wf_netcdfio.py (.../wf_netcdfio.py) (revision e71ad2de03777dca02510a47f65dd6f7075a876b) @@ -606,6 +606,7 @@ vars: list of variables to get from file """ + self.fname = netcdffile if os.path.exists(netcdffile): self.dataset = netCDF4.Dataset(netcdffile, mode='r') else: @@ -622,31 +623,76 @@ self.fstep = 0 self.lstep = self.fstep + self.maxsteps + self.datetime = self.dataset.variables['time'][:] + if hasattr(self.dataset.variables['time'],'units'): + self.timeunits=self.dataset.variables['time'].units + else: + self.timeunits ='Seconds since 1970-01-01 00:00:00' + if hasattr(self.dataset.variables['time'], 'calendar'): + self.calendar= self.dataset.variables['time'].calendar + else: + self.calendar ='gregorian' + self.datetimelist=netCDF4.num2date(self.datetime,self.timeunits, calendar=self.calendar) + try: self.x = self.dataset.variables['x'][:] except: self.x = self.dataset.variables['lon'][:] + # Now check Y values to see if we must flip the data try: self.y = self.dataset.variables['y'][:] except: self.y = self.dataset.variables['lat'][:] + # test if 1D or 2D array + if len(self.y.shape) == 1: + if self.y[0] > self.y[-1]: + self.flip = False + else: + self.flip = True + else: # not sure if this works + self.y = self.y[:][0] + if self.y[0] > self.y[-1]: + self.flip = False + else: + self.flip = True + x = _pcrut.pcr2numpy(_pcrut.xcoordinate(_pcrut.boolean(_pcrut.cover(1.0))), NaN)[0, :] y = _pcrut.pcr2numpy(_pcrut.ycoordinate(_pcrut.boolean(_pcrut.cover(1.0))), NaN)[:, 0] - (self.latidx,) = logical_and(self.x >= x.min(), self.x < x.max()).nonzero() - (self.lonidx,) = logical_and(self.y >= x.min(), self.y < y.max()).nonzero() + #Get average cell size + acc = diff(x).mean() * 0.25 # non-exact match needed becuase of possible rounding problems + if self.flip: + (self.latidx,) = logical_and(self.y[::-1] +acc >= y.min(), self.y[::-1] <= y.max() + acc).nonzero() + (self.lonidx,) = logical_and(self.x + acc >= x.min(), self.x <= x.max() + acc).nonzero() + else: + (self.latidx,) = logical_and(self.y +acc >= y.min(), self.y <= y.max() + acc).nonzero() + (self.lonidx,) = logical_and(self.x +acc >= x.min(), self.x <= x.max() + acc).nonzero() + if len(self.lonidx) != len(x): + logging.error("error in determining X coordinates in netcdf...") + logging.error("model expects: " + str(x.min()) + " to " + str(x.max())) + logging.error("got coordinates netcdf: " + str(self.x.min()) + " to " + str(self.x.max())) + logging.error("got len from netcdf x: " + str(len(x)) + " expected " + str(len(self.lonidx))) + raise ValueError("X coordinates in netcdf do not match model") + + if len(self.latidx) != len(y): + logging.error("error in determining Y coordinates in netcdf...") + logging.error("model expects: " + str(y.min()) + " to " + str(y.max())) + logging.error("got from netcdf: " + str(self.y.min()) + " to " + str(self.y.max())) + logging.error("got len from netcdf y: " + str(len(y)) + " expected " + str(len(self.latidx))) + raise ValueError("Y coordinates in netcdf do not match model") + 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'): + def gettimestep(self, timestep, logging, var='P', tsdatetime=None): """ Gets a map for a single timestep. reads data in blocks assuming sequential access @@ -655,19 +701,19 @@ 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, self.latidx.min():self.latidx.max() + 1, - self.lonidx.min():self.lonidx.max() + 1] + + if var in self.dataset.variables: + if tsdatetime != None: + if tsdatetime.replace(tzinfo=None) != self.datetimelist[ncindex].replace(tzinfo=None): + logging.warn("Date/time of state (" + var + " in " + self.fname + ")does not match. Wanted " + str(tsdatetime) + " got " + str(self.datetimelist[ncindex])) + + np_step = self.dataset.variables[var][ncindex, self.latidx.min():self.latidx.max() + 1, + self.lonidx.min():self.lonidx.max() + 1] + miss = float(self.dataset.variables[var]._FillValue) return numpy2pcr(Scalar, np_step, miss), True else: - logging.debug("Var (" + var + ") not found returning map with 0.0") + #logging.debug("Var (" + var + ") not found returning map with 0.0") return cover(scalar(0.0)), False Index: wflow-py/wflow/wflow_sbm.py =================================================================== diff -u -rab8caa7e1395fd178532feafe2fdae06f875f1b1 -re71ad2de03777dca02510a47f65dd6f7075a876b --- wflow-py/wflow/wflow_sbm.py (.../wflow_sbm.py) (revision ab8caa7e1395fd178532feafe2fdae06f875f1b1) +++ wflow-py/wflow/wflow_sbm.py (.../wflow_sbm.py) (revision e71ad2de03777dca02510a47f65dd6f7075a876b) @@ -349,8 +349,9 @@ if hasattr(self,'ReserVoirLocs'): states.append('ReservoirVolume') - if self.nrpaddyirri > 0: - states.append('PondingDepth') + if hasattr(self,'nrpaddyirri'): + if self.nrpaddyirri > 0: + states.append('PondingDepth') return states def supplyCurrentTime(self):