diff --git a/dfm_tools/interpolate_grid2bnd.py b/dfm_tools/interpolate_grid2bnd.py index 5bcaf44b1..0f8dad876 100644 --- a/dfm_tools/interpolate_grid2bnd.py +++ b/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') diff --git a/dfm_tools/modelbuilder.py b/dfm_tools/modelbuilder.py index 1bafdfd53..18cb19a4d 100644 --- a/dfm_tools/modelbuilder.py +++ b/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): diff --git a/docs/whats-new.md b/docs/whats-new.md index ab0f840ed..75a4d012d 100644 --- a/docs/whats-new.md +++ b/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) diff --git a/tests/examples/preprocess_ini_cmems_to_nc.py b/tests/examples/preprocess_ini_cmems_to_nc.py deleted file mode 100644 index e871fc884..000000000 --- a/tests/examples/preprocess_ini_cmems_to_nc.py +++ /dev/null @@ -1,39 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Created on Tue Sep 27 12:14:14 2022 - -@author: veenstra -""" - -import os -import glob -import datetime as dt -import xarray as xr - -#TODO: merge with other ini script and make generic for getting an inifield out of CMEMS/etc regulargrid Dataset or a 2D/3D FM map/rst Dataset - -tSimStart = dt.datetime(1998,11,1) -dir_data = r'p:\i1000668-tcoms\03_newModel\01_input\02_bnd\data_opendap' #folder containing CMEMS so and thetao netcdf files - -dir_out = '.' - -file_nc_list_so = glob.glob(f'{dir_data}\\cmems_so_1998-1*.nc') -file_nc_list_thetao = glob.glob(f'{dir_data}\\cmems_thetao_1998-1*.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) -data_xr_times_pd = data_xr.time.to_series() - -interp = False -if interp: #this would be the proper way to do it, but FM needs two timesteps for some reason - print('ds.interp()') - data_xr_ontime = data_xr.interp(time=[tSimStart],kwargs=dict(bounds_error=True)) #bounds_error makes sure, outofbounds time results in "ValueError: A value in x_new is below the interpolation range." -else: - print('ds.sel()') - 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) - diff --git a/tests/examples_workinprogress/workinprogress_modelbuilder.py b/tests/examples_workinprogress/workinprogress_modelbuilder.py index 5e894afb4..07f9d89b6 100644 --- a/tests/examples_workinprogress/workinprogress_modelbuilder.py +++ b/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,6 +143,7 @@ # 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, @@ -151,7 +151,7 @@ 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) @@ -159,20 +159,20 @@ #%% 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 diff --git a/tests/test_external_packages.py b/tests/test_external_packages.py index 622d8596a..dff1beb7a 100644 --- a/tests/test_external_packages.py +++ b/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) + diff --git a/tests/test_interpolate_grid2bnd.py b/tests/test_interpolate_grid2bnd.py index b82fa77c4..001f2a781 100644 --- a/tests/test_interpolate_grid2bnd.py +++ b/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(): diff --git a/tests/test_modelbuilder.py b/tests/test_modelbuilder.py new file mode 100644 index 000000000..f5bf98af0 --- /dev/null +++ b/tests/test_modelbuilder.py @@ -0,0 +1,70 @@ +# -*- coding: utf-8 -*- +""" +Created on Thu Oct 26 16:08:26 2023 + +@author: veenstra +""" + +import os +import pytest +import dfm_tools as dfmt +import hydrolib.core.dflowfm as hcdfm +import xarray as xr +from dfm_tools.hydrolib_helpers import get_ncbnd_construct + + +@pytest.mark.systemtest +def test_cmems_nc_to_ini(): + + # TODO: create fixture + from tests.test_interpolate_grid2bnd import cmems_dataset_4times + ds1 = cmems_dataset_4times().isel(time=slice(None,2)) + ds2 = cmems_dataset_4times().isel(time=slice(2,None)) + + dir_pattern = "./temp_cmems_2day_*.nc" + file_nc1 = dir_pattern.replace("*","sal_p1") + file_nc2 = dir_pattern.replace("*","sal_p2") + file_nc3 = dir_pattern.replace("*","tem_p1") + file_nc4 = dir_pattern.replace("*","tem_p2") + ds1.to_netcdf(file_nc1) + ds2.to_netcdf(file_nc2) + ds1.rename({"so":"thetao"}).to_netcdf(file_nc3) + ds2.rename({"so":"thetao"}).to_netcdf(file_nc4) + + ext_old = hcdfm.ExtOldModel() + + ext_old = dfmt.cmems_nc_to_ini(ext_old=ext_old, + dir_output='.', + list_quantities=["salinitybnd"], + tstart="2020-01-01", + dir_pattern=dir_pattern) + + file_expected = "./nudge_salinity_temperature_2020-01-01_00-00-00.nc" + + times_expected = ['2019-12-31 12:00:00', '2020-01-01 12:00:00'] + + assert os.path.exists(file_expected) + ds_out = xr.open_dataset(file_expected) + + times_actual = ds_out.time.to_pandas().dt.strftime("%Y-%m-%d %H:%M:%S").tolist() + assert times_expected == times_actual + + assert "so" in ds_out.data_vars + assert ds_out.so.isnull().sum().load() == 0 + + ncbnd_construct = get_ncbnd_construct() + varn_depth = ncbnd_construct['varn_depth'] + assert varn_depth in ds_out.coords + # the below is inconsistent since depth is actually defined positive up, but FM requires this for inifields: https://issuetracker.deltares.nl/browse/UNST-7455 + assert ds_out[varn_depth].attrs['positive'] == 'down' + + # cleanup + del ds_out + del ds1 + del ds2 + os.remove(file_expected) + os.remove(file_nc1) + os.remove(file_nc2) + os.remove(file_nc3) + os.remove(file_nc4) + \ No newline at end of file