Improve interpolation of uds/ugrid to plipoints

[veenstrajelmer]

After fixing #912 and #911 the conversion_dict can be applied to a uds (ugriddataset) also. Furthermore, there are some additional improvements to the bc files (like depth correction and quantity renaming) via the conventions/conversion public functions:

import dfm_tools as dfmt
import hydrolib.core.dflowfm as hcdfm
import matplotlib.pyplot as plt
plt.close("all")
import xarray as xr

# read plifile as gdf
file_pli = r'P:\11209133-ultfarms\02-Simulations\North_Sea_runs_Loreta\pre-processing\offshore_boundaries\FuMA_scripts\DCSM-FM_OB_all_20181108.pli' # polygon containing open boundary forcing points
pli_obj = hcdfm.PolyFile(file_pli)
pli_gdf = dfmt.PolyFile_to_geodataframe_points(pli_obj)

# only take first points to speed up for testing
pli_gdf_test = pli_gdf.iloc[:6]
time_slice = slice(None,5) # use None to use all timesteps

scen = 'historical' # test
# yr = 2005#, 2006, 2007, 2008] # test

conversion_dict = dfmt.get_conversion_dict()
conversion_dict = {
    # conversion is phyc in mol/m3 to newvar in g/m3
    'tracerbndOXY'        : {'ncvarname': 'o2',          'unit': 'g/m3', 'conversion': 32.0 },
    'tracerbndNO3'        : {'ncvarname': 'no3',         'unit': 'g/m3', 'conversion': 14.0 },
    'tracerbndPO4'        : {'ncvarname': 'po4',         'unit': 'g/m3', 'conversion': 30.97 },
    'tracerbndSi'         : {'ncvarname': 'si',          'unit': 'g/m3', 'conversion': 28.08},
    'tracerbndPON1'       : {'ncvarname': 'phyc',        'unit': 'g/m3', 'conversion': 14.0}, # Caution: this empirical relation might not be applicable to your use case
    'tracerbndPOP1'       : {'ncvarname': 'phyc',        'unit': 'g/m3', 'conversion': 30.97}, # Caution: this empirical relation might not be applicable to your use case
    'tracerbndPOC1'       : {'ncvarname': 'phyc',        'unit': 'g/m3', 'conversion': 14.0 * (106 / 16)}, # Caution: this empirical relation might not be applicable to your use case
    'salinitybnd'         : {'ncvarname': 'so'},         #'1e-3'
    'temperaturebnd'      : {'ncvarname': 'thetao'},         #'degC'
    'ux'                  : {'ncvarname': 'uo'},          #'m/s'
    'uy'                  : {'ncvarname': 'vo'},          #'m/s'
    'waterlevelbnd'       : {'ncvarname': 'zos'},         #'m' #steric
    'tide'                : {'ncvarname': ''},            #'m' #tide (dummy entry)
}

ext_new = hcdfm.ExtModel()

quan_bnd_list = ['uxuyadvectionvelocitybnd','salinitybnd','temperaturebnd','waterlevelbnd']
# quan_bnd_list = ['tracerbndNO3'] # fails on "KeyError: 'units'" (will be fixed in dfm_tools 0.26.0)

for quantity in quan_bnd_list:
    # quantity_list = ['ux','uy']
    # TODO: consider avoiding the need to use private functions
    quantity_list = dfmt.modelbuilder.get_quantity_list(quantity)
    
    ds_points_list = []
    for quantity_one in quantity_list:
        ncvarname = conversion_dict[quantity_one]['ncvarname']
        
        dir_pattern = fr'P:\archivedprojects\11206304-futuremares-rawdata-preps\data\CMIP6_BC\CMCC-ESM2\{ncvarname}_Omon_CMCC-ESM2_{scen}_r1i1p1f1_gn_*.nc'
        
        uds = dfmt.open_dataset_curvilinear(dir_pattern, convert_360to180=True)
        
        # update refdate
        uds.time.encoding['units'] = 'days since 2004-01-01'
    
        uds = uds.ugrid.to_nonperiodic(xmax=180)
        
        # subset in time to speed up for testing
        uds_sel = uds.isel(time=time_slice)
        
        # # plot (first subset in space for visualization purposes)
        # buffer = 5
        # xmin = pli_gdf_test.geometry.x.min() - buffer
        # xmax = pli_gdf_test.geometry.x.max() + buffer
        # ymin = pli_gdf_test.geometry.y.min() - buffer
        # ymax = pli_gdf_test.geometry.y.max() + buffer
        # fig, ax = plt.subplots()
        # uds.ugrid.sel(x=slice(xmin,xmax), y=slice(ymin,ymax))[ncvarname].isel(time=0, lev=0).ugrid.plot(ax=ax)
        # dfmt.plot_coastlines(ax=ax, res='i')
        # pli_gdf_test.plot(ax=ax, color='r')
    
        # apply conventions and conversions
        # TODO: consider avoiding the need to use private functions
        uds_conv = dfmt.interpolate_grid2bnd.ds_apply_conventions(uds_sel)
        uds_conv = dfmt.interpolate_grid2bnd.ds_apply_conversion_dict(uds_conv, conversion_dict=conversion_dict, quantity=quantity_one)
        
        # interpolate to points of polyfile
        ds_points = dfmt.interp_uds_to_plipoints(uds=uds_conv, gdf=pli_gdf_test)
        
        # keep only relevant variable
        ds_points = ds_points[quantity_one]
        
        ds_points_list.append(ds_points)
    
    ds_points_onevar = xr.merge(ds_points_list)
    
    # convert to bc/forcing hydrolib-core object
    ForcingModel_object = dfmt.plipointsDataset_to_ForcingModel(ds_points_onevar)
    # save to bc file
    ForcingModel_object.save(f"test_{quantity}_{scen}.bc")
    
    # append to extfile
    boundary_object = hcdfm.Boundary(quantity=quantity,
                                     locationfile=file_pli, #placeholder, will be replaced later on
                                     forcingfile=ForcingModel_object)
    
    # TODO: consider doing this per polyline (required if there is more than one in the file)
    # with ext_add_boundary_object_per_polyline
    ext_new.boundary.append(boundary_object)

ext_new.save(f"test_{scen}.ext")

Todo:

• It does not yet work for multiple datasets in time
• drop 0-sized cells in dfmt.open_dataset_curvilinear() #926
• fix slow performance of dfmt.open_dataset_curvilinear()` due to incorrectly added time dimension #928
• It does not yet work for ux/uy boundaries >> it being a vector quantity in FM is inconvenient, but the CMCC datasets also have ux/uy on different locations, so these would have to be interpolated separately and then merged. In other functions we combine the raw data (having the same coords) and interpolate the entire dataset to points.
• move drop_vars to dfmt.open_dataset_curvilinear(), but dfmt.interpolate_grid2bnd.ds_apply_conventions() will raise "KeyError: "No variable named 'longitude'. Variables on the dataset include...". >> Cleanup dataset returned by dfmt.open_dataset_curvilinear() #930
• refactor open_dataset_extra() #945
• fix incorrect number of returned points by dfmt.interp_uds_to_plipoints() with outofbounds plipoints #927
• instead of manual reconstruction, try xu.from_structured_multicoord() instead >> "ValueError: The input coordinate is not monotonic"
• also add the forcings to an extmodel
• units attribute is not added during conversion of CMCC #964
• currently requires using some newly implemented private functions and user code is still quite complex in case of ux/uy-vectors. This can be resolved by making a wrapper function similar to dfmt.open_dataset_extra(), but make sure to align functionalities across similar dfm_tools functions and avoid duplicate code. >> This would become quite specific given the different latlon locs of ux/uy and the sortof curvigrid, so maybe not useful to add a (dedicated) wrapper function
• open_dataset_curvilinear() drops faces with some nan nodes or some duplicate nodes (also triangles), consider fixing this. Is used for CMCC datasets. There were still 0-size cells which required drop 0-sized cells in dfmt.open_dataset_curvilinear() #926. Why is this the case? If the cell size is zero, all nodes should be equal which should be catched. >> size is also zero if the length of one of the sides of a square is zero, which won't have all-equal nodes.