618 add initial fields for tracers (#619)

* initial values in extfile, two timesteps, z conversion

* added testcase

* updated whatsnew

* auto remove temporary pytest nc files

* remove outdated example file

dfm_tools/interpolate_grid2bnd.py

@@ -248,6 +248,8 @@ def open_dataset_extra(dir_pattern, quantity, tstart, tstop, conversion_dict=Non
 
  if quantity=='uxuyadvectionvelocitybnd': #T3Dvector
  quantity_list = ['ux','uy']
+ elif isinstance(quantity, list):
+ quantity_list = quantity
  else:
  quantity_list = [quantity]
 
@@ -336,10 +338,11 @@ def open_dataset_extra(dir_pattern, quantity, tstart, tstop, conversion_dict=Non
  data_xr_vars.time.encoding['units'] = refdate_str
 
  if varn_depth in data_xr_vars.coords:
- #make negative down
+ #make positive up
  if 'positive' in data_xr_vars[varn_depth].attrs.keys():
  if data_xr_vars[varn_depth].attrs['positive'] == 'down': #attribute appears in CMEMS, GFDL and CMCC, save to assume presence?
  data_xr_vars[varn_depth] = -data_xr_vars[varn_depth]
+ data_xr_vars[varn_depth].attrs['positive'] = 'up'
  #optional reversing depth dimension for comparison to coastserv
  if reverse_depth:
  print('> reversing depth dimension')

dfm_tools/modelbuilder.py

@@ -9,14 +9,13 @@
 """
 
 import os
-import xarray as xr
 import pandas as pd
 import dfm_tools as dfmt
 import hydrolib.core.dflowfm as hcdfm
 import datetime as dt
-import glob
 from hydrolib.core.dimr.models import DIMR, FMComponent, Start
 import warnings
+from dfm_tools.hydrolib_helpers import get_ncbnd_construct
 
 
 def cmems_nc_to_bc(ext_bnd, list_quantities, tstart, tstop, file_pli, dir_pattern, dir_output, refdate_str):
@@ -50,46 +49,79 @@ def cmems_nc_to_bc(ext_bnd, list_quantities, tstart, tstop, file_pli, dir_patter
 
  return ext_bnd
 
+
+def preprocess_ini_cmems_to_nc(**kwargs):
+ raise DeprecationWarning("`dfmt.preprocess_ini_cmems_to_nc()` was deprecated, use `cmems_nc_to_ini()` instead")
+
+
+def cmems_nc_to_ini(ext_old, dir_output, list_quantities, tstart, dir_pattern, conversion_dict=None):
 
-def preprocess_ini_cmems_to_nc(ext_old, tstart, dir_data, dir_out):
-
- file_nc_list_so = glob.glob(f'{dir_data}\\cmems_so_*.nc')
- file_nc_list_thetao = glob.glob(f'{dir_data}\\cmems_thetao_*.nc')
- file_nc_list = file_nc_list_so + file_nc_list_thetao
-
- print(f'opening {len(file_nc_list)} datasets')
- data_xr = xr.open_mfdataset(file_nc_list)
-
- # fill nans
- # start with lat/lon to avoid values from shallow coastal areas in deep layers
- # first interpolate nans to get smooth filling of e.g. islands, this cannot fill nans at the edge of the dataset
- data_xr = data_xr.interpolate_na(dim='latitude').interpolate_na(dim='longitude')
-
- # then use bfill/ffill to fill nans at the edge for lat/lon/depth
- data_xr = data_xr.ffill(dim='latitude').bfill(dim='latitude')
- data_xr = data_xr.ffill(dim='longitude').bfill(dim='longitude')
- data_xr = data_xr.ffill(dim='depth').bfill(dim='depth')
-
- tSimStart = pd.Timestamp(tstart)
- data_xr_ontime = data_xr.sel(time=slice(tSimStart-dt.timedelta(days=1),tSimStart+dt.timedelta(days=1)))
-
- print('writing file')
- outFile = os.path.join(dir_out,f'InitialField_{tSimStart.strftime("%Y-%m-%d_%H-%M-%S")}.nc')
- data_xr_ontime.to_netcdf(outFile,format="NETCDF4_CLASSIC") #TODO: why the format setting?
+ if conversion_dict is None:
+ conversion_dict = dfmt.get_conversion_dict()
 
- #append forcings to ext
- # 3D initialsalinity/initialtemperature fields are silently ignored, initial 3D conditions are only possible via nudging 1st timestep via quantity=nudge_salinity_temperature
- forcing_saltem = hcdfm.ExtOldForcing(quantity='nudge_salinity_temperature',
- filename=outFile,
- filetype=hcdfm.ExtOldFileType.NetCDFGridData,
- method=hcdfm.ExtOldMethod.InterpolateTimeAndSpaceSaveWeights, #3
- operand=hcdfm.Operand.override, #O
- )
- ext_old.forcing.append(forcing_saltem)
+ tstart_pd = pd.Timestamp(tstart)
+ tstart_str = tstart_pd.strftime("%Y-%m-%d_%H-%M-%S")
+ tstart_round, tstop_round = dfmt.round_timestamp_to_outer_noon(tstart,tstart)
+ for quan_bnd in list_quantities:
+
+ if quan_bnd in ["temperaturebnd","uxuyadvectionvelocitybnd"]:
+ # silently skipped, temperature is handled with salinity, uxuy not supported
+ continue
+ elif quan_bnd=="salinitybnd":
+ # 3D initialsalinity/initialtemperature fields are silently ignored
+ # initial 3D conditions are only possible via nudging 1st timestep via quantity=nudge_salinity_temperature
+ data_xr = dfmt.open_dataset_extra(dir_pattern=dir_pattern, quantity=["salinitybnd","temperaturebnd"],
+ tstart=tstart_round, tstop=tstop_round,
+ conversion_dict=conversion_dict)
+ data_xr = data_xr.rename_vars({"salinitybnd":"so", "temperaturebnd":"thetao"})
+ quantity = "nudge_salinity_temperature"
+ varname = None
+ else:
+ data_xr = dfmt.open_dataset_extra(dir_pattern=dir_pattern, quantity=quan_bnd,
+ tstart=tstart_round, tstop=tstop_round,
+ conversion_dict=conversion_dict)
+ quantity = f'initial{quan_bnd.replace("bnd","")}'
+ varname = quantity
+ data_xr = data_xr.rename_vars({quan_bnd:quantity})
+
+ # open_dataset_extra converted depths from positive down to positive up, including update of the "positive" attribute
+ # TODO: this correctly updated attr negatively impacts model results when using netcdf inifields, so we revert it here
+ # https://issuetracker.deltares.nl/browse/UNST-7455
+ ncbnd_construct = get_ncbnd_construct()
+ varn_depth = ncbnd_construct['varn_depth']
+ if varn_depth in data_xr.coords:
+ data_xr[varn_depth].attrs['positive'] = 'down'
+
+ # subset two times. interp to tstart would be the proper way to do it, 
+ # but FM needs two timesteps for nudge_salinity_temperature and initial waq vars
+ data_xr = data_xr.sel(time=slice(tstart_round, tstop_round))
+
+ # fill nans, start with lat/lon to avoid values from shallow coastal areas in deep layers
+ # first interpolate nans to get smooth filling of e.g. islands, this cannot fill nans at the edge of the dataset
+ data_xr = data_xr.interpolate_na(dim='latitude').interpolate_na(dim='longitude')
+
+ # then use bfill/ffill to fill nans at the edge for lat/lon/depth
+ data_xr = data_xr.ffill(dim='latitude').bfill(dim='latitude')
+ data_xr = data_xr.ffill(dim='longitude').bfill(dim='longitude')
+ data_xr = data_xr.ffill(dim='depth').bfill(dim='depth')
+
+
+ print('writing file')
+ file_output = os.path.join(dir_output,f"{quantity}_{tstart_str}.nc")
+ data_xr.to_netcdf(file_output)
+
+ #append forcings to ext
+ forcing_saltem = hcdfm.ExtOldForcing(quantity=quantity,
+ varname=varname,
+ filename=file_output,
+ filetype=hcdfm.ExtOldFileType.NetCDFGridData,
+ method=hcdfm.ExtOldMethod.InterpolateTimeAndSpaceSaveWeights, #3
+ operand=hcdfm.Operand.override) #O
+ ext_old.forcing.append(forcing_saltem)
 
  return ext_old
- 
- 
+
+
 def preprocess_merge_meteofiles_era5(ext_old, varkey_list, dir_data, dir_output, time_slice):
 
  if isinstance(varkey_list[0], list):

docs/whats-new.md

@@ -1,6 +1,10 @@
 ## UNRELEASED
 
+### Feat
+- support for initial fields for variables other than salinity/temperature with `dfmt.cmems_nc_to_ini()` by [@veenstrajelmer](https://github.com/veenstrajelmer) in [#619](https://github.com/Deltares/dfm_tools/pull/619)
+
 ### Fix
+- fix initial fields for salinity/temperature with `dfmt.cmems_nc_to_ini()` by [@veenstrajelmer](https://github.com/veenstrajelmer) in [#619](https://github.com/Deltares/dfm_tools/pull/619)
 - increased buffer in `dfmt.download_ERA5()` by [@veenstrajelmer](https://github.com/veenstrajelmer) in [#612](https://github.com/Deltares/dfm_tools/pull/612)
 - support for Polygon geometries in `dfmt.geodataframe_to_PolyFile()` by [@veenstrajelmer](https://github.com/veenstrajelmer) in [#610](https://github.com/Deltares/dfm_tools/pull/610)
 - fill nan-values in initial salinity/temperature netcdf dataset in `dfmt.preprocess_ini_cmems_to_nc()` by [@veenstrajelmer](https://github.com/veenstrajelmer) in [#617](https://github.com/Deltares/dfm_tools/pull/617)

tests/examples_workinprogress/workinprogress_modelbuilder.py

@@ -23,7 +23,6 @@
 is_geographic = True
 crs = 'EPSG:4326'
 
-inisaltem = True #initialsalinity/initialtemperature gives 33.8ppt uniform and sal instabilities right from the start of the model run. Proper way seems to be with iniwithnudge=2 and nudge_salinity_temperature, which gives ini sal/tem indeed but also instable. Providing nudge_salinity_temperature and iniwithnudge=0 gives more stable model but also inisal is 33.8 (not spatially varying) (is same as False maybe?)
 #TODO: salinity instable, also waterlevel and velocity magnitude are instable at northeast side of island (latter is with incorrect ordering/selection in extfile)
 """
 ** INFO : Min. salinity limited, number of cells Limmin = 20
@@ -144,35 +143,36 @@
 # or else "ERROR : update_ghostboundvals: not all ghost boundary flowlinks are being updated" is raised (https://issuetracker.deltares.nl/browse/UNST-7011).
 # Two waterlevelbnds need to share same physical plifile in order to be appended (https://issuetracker.deltares.nl/browse/UNST-5320).
 list_quantities = ['waterlevelbnd','salinitybnd','temperaturebnd','uxuyadvectionvelocitybnd','tracerbndNO3','tracerbndPON1']
+dir_pattern = os.path.join(dir_output_data_cmems,'cmems_{ncvarname}_*.nc')
 ext_new = dfmt.cmems_nc_to_bc(ext_bnd=ext_new,
  refdate_str=f'minutes since {ref_date} 00:00:00 +00:00',
  dir_output=dir_output,
  list_quantities=list_quantities,
  tstart=date_min,
  tstop=date_max, 
  file_pli=poly_file,
- dir_pattern=os.path.join(dir_output_data_cmems,'cmems_{ncvarname}_*.nc'))
+ dir_pattern=dir_pattern)
 
 #save new ext file
 ext_new.save(filepath=ext_file_new,path_style=path_style)
 
 
 #%% old ext
 
-# CMEMS - initial condition file
+# CMEMS - initial conditions
 ext_file_old = os.path.join(dir_output, f'{model_name}_old.ext')
 ext_old = hcdfm.ExtOldModel()
 
-if inisaltem:
- ext_old = dfmt.preprocess_ini_cmems_to_nc(ext_old=ext_old,
-  tstart=date_min,
-  dir_data=dir_output_data_cmems,
-  dir_out=dir_output)
+ext_old = dfmt.cmems_nc_to_ini(ext_old=ext_old,
+  dir_output=dir_output,
+ list_quantities=list_quantities,
+ tstart=date_min,
+ dir_pattern=dir_pattern)
 
 # ERA5 - download
 dir_output_data_era5 = os.path.join(dir_output_data,'ERA5')
 os.makedirs(dir_output_data_era5, exist_ok=True)
- 
+
 if ERA5_meteo_option == 1:
  varlist_list = [['msl','u10n','v10n','chnk']]
 elif ERA5_meteo_option == 2:
@@ -230,7 +230,10 @@
 mdu.physics.secchidepth = 4
 mdu.physics.salimax = 50
 mdu.physics.tempmax = 50
-if inisaltem:
+if 'nudge_salinity_temperature' in [x.quantity for x in ext_old.forcing]:
+ # initialsalinity/initialtemperature gives 33.8ppt uniform (so is not read)
+ # the only possible way is with iniwithnudge=2 and nudge_salinity_temperature in old extfile.
+ # for other variables we can use initialtracerbndNO3 for instance
  mdu.physics.iniwithnudge = 2 #TODO: commented in oldextfile in reference run, initial sal/tem profiles from deep layer were used instead (not yet derived, but 3D inifields also do not have an effect)
 
 mdu.wind.icdtyp = 4

tests/test_external_packages.py

@@ -5,6 +5,7 @@
 @author: veenstra
 """
 
+import os
 import pytest
 import xarray as xr
 import xugrid as xu
@@ -146,7 +147,6 @@ def test_xarray_decode_default_fillvals():
  #write file
  file_out = 'temp_fnc_default_fillvals_map.nc'
  ds.to_netcdf(file_out)
- ds.close()
 
  #open dataset with decode_fillvals
  try:
@@ -161,7 +161,15 @@ def test_xarray_decode_default_fillvals():
  "version is set as requirement")
 
  #this should be successful
- uds = dfmt.open_partitioned_dataset(file_out,decode_fillvals=True)
- fnc_new = uds.grid.face_node_connectivity
+ uds2 = dfmt.open_partitioned_dataset(file_out,decode_fillvals=True)
+ fnc_new = uds2.grid.face_node_connectivity
 
  assert fill_value_default in fnc_new
+
+ # cleanup
+ # del ds
+ # del uds2
+ # del fnc_new
+ # PermissionError: [WinError 32] The process cannot access the file because it is being used by another process: 'temp_fnc_default_fillvals_map.nc'
+ # os.remove(file_out)
+

tests/test_interpolate_grid2bnd.py

@@ -10,7 +10,6 @@
 import dfm_tools as dfmt
 import numpy as np
 import datetime as dt
-import pandas as pd
 import xarray as xr
 import shapely
 import geopandas as gpd
@@ -183,6 +182,10 @@ def test_open_dataset_extra_correctdepths():
 
  assert (np.abs(depth_actual - depth_expected) < 1e-9).all()
  assert len(ds_moretime_import.time) == 2
+
+ # cleanup
+ del ds_moretime_import
+ os.remove(file_nc)
 
 
 @pytest.mark.unittest
@@ -200,8 +203,10 @@ def test_open_dataset_extra_slightly_different_latlons():
  ds_lon[1] += 1e-8
  ds2['longitude'] = xr.DataArray(ds_lon,dims='longitude')
 
- ds1.to_netcdf('temp_cmems_2day_p1.nc')
- ds2.to_netcdf('temp_cmems_2day_p2.nc')
+ file_nc1 = 'temp_cmems_2day_p1.nc'
+ file_nc2 = 'temp_cmems_2day_p2.nc'
+ ds1.to_netcdf(file_nc1)
+ ds2.to_netcdf(file_nc2)
 
  try:
  ds = dfmt.open_dataset_extra('temp_cmems_2day_*.nc', quantity='salinitybnd', tstart='2020-01-01', tstop='2020-01-03')
@@ -212,6 +217,12 @@ def test_open_dataset_extra_slightly_different_latlons():
  # ValueError: cannot align objects with join='exact' where index/labels/sizes are not equal along these coordinates (dimensions): 'longitude' ('longitude',)
  pass # this is expected, so pass
 
+ # cleanup
+ del ds1
+ del ds2
+ os.remove(file_nc1)
+ os.remove(file_nc2)
+
 
 @pytest.mark.unittest
 def test_interp_regularnc_to_plipointsDataset():