Error setting up cold start ICs (Start+count exceeds dimension bound)

[NicWayand]

SUMMA v3.2.0 Latest, gfortran, ubuntu:xenial

I am trying to create a cold start IC file (summa_zInitialCond.nc) and getting the below errror. I suspect its because I have something wrong in the dimensions of the file. I am starting with the example file settings/syntheticTestCases/wigmosta1999/summa_zInitialCond.nc Then simply extending it to 1000 hrus for my test run. I also remove what I believe are variables not required for initialization. So my summa_zInitialCond.nc looks like:

The error:

trim(nf90_strerror(err) = NetCDF: Start+count exceeds dimension bound
 read_icond/[NetCDF: Start+count exceeds dimension bound]


FATAL ERROR: summa_readRestart/read_icond/[NetCDF: Start+count exceeds dimension bound]: problem getting the data for variable dt_init

Note: The following floating-point exceptions are signalling: IEEE_INVALID_FLAG IEEE_DIVIDE_BY_ZERO

Are there any existing scripts or best practices for creating cold start IC files?

[andywood]

You're probably missing the dt_init variable. Anyway, I have a little script to generate them, using another file that lists the hruIds (eg the trialParam file) as a template. I still use this but there could be changes in SUMMA since I last used it. It doesn't initialize GRU routing params because summa will initialize those if they aren't in the state file. Note, there's a spot where layer depths are hardwired in this. Andy ========= #!/usr/bin/env python ''' Create a vector cold state file for SUMMA from constant values''' # Author: Andy Wood, Feb 2017 # Notes: quick n dirty to generate constant initial states across a domain # all values hardwired, just gets HRU index from an existing parameter file # improvements: could read another cold state file to get list of variables # to populate; or read a metadata dictionary (names, types, etc) # no mapping required here -- but one could map another resolution vals, similar # to the param mapping scripts # # Requirements: run with a python (eg miniconda) 3.x that includes netCDF4 # ========================================================================= import sys, os, time, getopt import numpy as np import netCDF4 as nc4 #import xarray as xr ######################################################################## # Subroutines # ######################################################################## def getNetCDFData(fn, varname): """Read <varname> variables available to be mapped from NetCDF <fn> """ f = nc4.Dataset(fn,'r') data = f.variables[varname][:] f.close() # ds = xr.open_dataset(fn) # data = ds[varname] return data def getOutputPolyIDs(nc_file): outPolyIDs = getNetCDFData(nc_file, 'hruId') print("read output outPolyIds ('hruId') from example domain file") return outPolyIDs # write variables to netcdf output file def writeNC_state_vars(nc_out, newVarName, newVarDim, newVarType, newVarVals, fillVal): """ Write <vars>[hru] array in netCDF4 file,<fn> and variable of <varname> """ print("adding data") ncvar = nc_out.createVariable(newVarName, newVarType, (newVarDim, 'hru',),fill_value=fillVal) ncvar[:] = newVarVals # store data in netcdf file # write dimensions and dimension variables to netcdf output file def writeNC_dims(fn, scalarv, midSoil, midToto, ifcToto, hrus, hru_type): """ Write <vars>[hru] array in netCDF4 file,<fn> and variable of <varname> """ print("writing output file") nc_out = nc4.Dataset(fn, 'w', format='NETCDF4') # Create dimensions dim_hru = nc_out.createDimension('hru', len(hrus)) dim_scalarv = nc_out.createDimension('scalarv', scalarv) dim_midSoil = nc_out.createDimension('midSoil', midSoil) dim_midToto = nc_out.createDimension('midToto', midToto) dim_ifcToto = nc_out.createDimension('ifcToto', ifcToto) # --- Create HRU ID variable (can be either int or string) if hru_type == 'str': # string HRU (need to add string length) max_strlen = 20 # EC dim_str = nc_out.createDimension('strlen', max_strlen) hruId = nc_out.createVariable('hruId', 'S1', ('hru', 'strlen'),fill_value='-999') print("writing hru_type", hru_type) hruId[:] = nc4.stringtochar(np.asarray(hrus, dtype='S{}'.format(max_strlen))) elif hru_type == 'int64': # integer HRU hruId = nc_out.createVariable('hruId', 'i8', ('hru', ), fill_value='-999') print("writing hru_type", hru_type) hruId[:] = hrus #hruId[:] = np.asarray(hrus, dtype='int64') elif hru_type == 'int': # integer HRU hruId = nc_out.createVariable('hruId', 'i4', ('hru', ), fill_value='-999') print("writing hru_type", hru_type) hruId[:] = hrus #hruId[:] = np.asarray(hrus, dtype='int') else: # not recognized sys.exit("ERROR, hru_type not recognized: must be str, int64, or int") # add attribute hruId.long_name = 'USGS HUC12 ID' return nc_out ############################################ # Main # ############################################ use = ''' Usage: %s -[h] <existing_inputfile_with_hruId <output_netCDF> <hru_type (int|int64|str)> -h help ''' if __name__ == '__main__': def usage(): sys.stderr.write(use % sys.argv[0]) sys.exit(1) try: (opts, args) = getopt.getopt(sys.argv[1:], 'h') except getopt.error: usage() verbose = False grid_info = False proj_info = True for (opt,val) in opts: if opt == '-h': usage() elif opt == '-v': verbose = True else: raise OptionError(opt) usage() if len(args) == 3: nc_example_name = args[0] # template file (other params) nc_out_name = args[1] # output cold-state file hru_type = args[2] # 'int', 'int64' or 'string' # nc_in_name = args[4] # existing cold-state file # hardwired to forcing formats (hru index rather than grid) outPolyIDs=getOutputPolyIDs(nc_example_name) nOutPolygons = len(outPolyIDs) # === now start to create the cold state variables using the variable template === # settings (hardwire for now) # -- 3 layer scalarv = 1 midSoil = 3 midToto = 3 ifcToto = midToto+1 dT = 10800 # timestep of forcings in seconds #dT = 3600 # timestep of forcings in seconds #lyrDepth = [0.1, 0.5, 0.9] #lyrHeight = [0.0, 0.1, 0.6, 1.5] lyrDepth = [0.1, 0.3, 0.6] lyrHeight = [0.0, 0.1, 0.4, 1.0] # initialize netcdf file by storing dimensions and hru variable nc_out = writeNC_dims(nc_out_name, scalarv, midSoil, midToto, ifcToto, outPolyIDs, hru_type) # === now loop through variables and write # this could be done by looping through the input state file and xferring values # could also read vars and intFillVal = -999 realFillVal = -999.0 # nSoil, nSnow newVarVals = np.full((1,nOutPolygons), midSoil, dtype='int') writeNC_state_vars(nc_out, 'nSoil', 'scalarv', 'i4', newVarVals, intFillVal) newVarVals = np.zeros((1,nOutPolygons), dtype='int') writeNC_state_vars(nc_out, 'nSnow', 'scalarv', 'i4', newVarVals, intFillVal) # dT newVarVals = np.full((1,nOutPolygons), dT) writeNC_state_vars(nc_out, 'dt_init', 'scalarv', 'f8', newVarVals, realFillVal) # SWE, SnowDepth, SfcMeltPond, SnowAlbedo, CanopyLiq, CanopyIce newVarVals = np.zeros((1,nOutPolygons)) writeNC_state_vars(nc_out, 'scalarSWE', 'scalarv', 'f8', newVarVals, realFillVal) writeNC_state_vars(nc_out, 'scalarSnowDepth', 'scalarv', 'f8', newVarVals, realFillVal) writeNC_state_vars(nc_out, 'scalarSfcMeltPond', 'scalarv', 'f8', newVarVals, realFillVal) writeNC_state_vars(nc_out, 'scalarSnowAlbedo', 'scalarv', 'f8', newVarVals, realFillVal) writeNC_state_vars(nc_out, 'scalarCanopyLiq', 'scalarv', 'f8', newVarVals, realFillVal) writeNC_state_vars(nc_out, 'scalarCanopyIce', 'scalarv', 'f8', newVarVals, realFillVal) # CanairTemp, CanopyTemp newVarVals = np.full((1,nOutPolygons), 283.16) writeNC_state_vars(nc_out, 'scalarCanairTemp', 'scalarv', 'f8', newVarVals, realFillVal) writeNC_state_vars(nc_out, 'scalarCanopyTemp', 'scalarv', 'f8', newVarVals, realFillVal) # AquiferStorage newVarVals = np.full((1,nOutPolygons), 1.0) writeNC_state_vars(nc_out, 'scalarAquiferStorage', 'scalarv', 'f8', newVarVals, realFillVal) # layer MatricHead newVarVals = np.full((midSoil,nOutPolygons), -1.0) writeNC_state_vars(nc_out, 'mLayerMatricHead', 'midSoil', 'f8', newVarVals, realFillVal) # layer Temp newVarVals = np.full((midToto,nOutPolygons), 283.16) writeNC_state_vars(nc_out, 'mLayerTemp', 'midToto', 'f8', newVarVals, realFillVal) # layer VolFracLiq newVarVals = np.full((midToto,nOutPolygons), 0.2) writeNC_state_vars(nc_out, 'mLayerVolFracLiq', 'midToto', 'f8', newVarVals, realFillVal) # layer VolFracIce newVarVals = np.full((midToto,nOutPolygons), 0.0) writeNC_state_vars(nc_out, 'mLayerVolFracIce', 'midToto', 'f8', newVarVals, realFillVal) # layer Depth, Height newVarVals = np.full((nOutPolygons,midToto), lyrDepth).transpose() writeNC_state_vars(nc_out, 'mLayerDepth', 'midToto', 'f8', newVarVals, realFillVal) newVarVals = np.full((nOutPolygons,ifcToto), lyrHeight).transpose() writeNC_state_vars(nc_out, 'iLayerHeight', 'ifcToto', 'f8', newVarVals, realFillVal) nc_out.close() else: usage()

…

On Thu, Aug 3, 2023 at 10:06 AM Nic Wayand ***@***.***> wrote: SUMMA v3.2.0 Latest, gfortran, ubuntu:xenial I am trying to create a cold start IC file (summa_zInitialCond.nc) and getting the below errror. I suspect its because I have something wrong in the dimensions of the file. I am starting with the example file ettings/syntheticTestCases/wigmosta1999/summa_zInitialCond.nc Then simply extending it to 1000 hrus for my test run. I also remove what I believe are variables not required for initialization. So my summa_zInitialCond.nc looks like: [image: Screenshot 2023-08-03 at 9 06 02 AM] <https://user-images.githubusercontent.com/1117224/258170946-670d91fa-79dd-4942-80d2-fa616215369f.png> The error: trim(nf90_strerror(err) = NetCDF: Start+count exceeds dimension bound read_icond/[NetCDF: Start+count exceeds dimension bound] FATAL ERROR: summa_readRestart/read_icond/[NetCDF: Start+count exceeds dimension bound]: problem getting the data for variable dt_init Note: The following floating-point exceptions are signalling: IEEE_INVALID_FLAG IEEE_DIVIDE_BY_ZERO Are there any existing scripts or best practices for creating cold start IC files? — Reply to this email directly, view it on GitHub <#544>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABIKARPAZESWI5IR732NPW3XTPEAXANCNFSM6AAAAAA3DABEJQ> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[NicWayand]

Thanks Andy for sharing your script, I will test it out.

I do have dt_init defined in the IC file though, does it need to be set anywhere else? Also, it wasn't clear from the definition if dt_init should be set at the input weather forcing timestep (1 hour = 3600 s) or at the best guess of sub-timestep (60s is what I saw in the example).

[NicWayand]

Ok I found the issue, due to the way I was modifying the IC xr.Dataset() in python, the order of dims for each variable was not the same.

bad and get error:

mLayerVolFracIce (midToto, hru)
scalarCanairTemp(hru, scalarv)

Good no error:

mLayerVolFracIce (midToto, hru)
scalarCanairTemp(scalarv, hru)

[wknoben]

I've seen this happen with the forcing files as well (Python dims needing to be swapped to make it work in Fortran), so just a heads-up in case you run into that. There's more SUMMA setup code examples here in case that helps: https://github.com/CH-Earth/CWARHM