Update ugrid conventions #114
Assignees
Labels
enhancement
New feature or request
Comments
|
Thx @veenstrajelmer. Did you do this for OLDIO or new scribed IO? For the latter we have UGRID part added, so we can easily add those changes. |
|
Thanks for your reply. I realized after the creation of this issue that it does not have enough context. At the moment i cannot answer your questions, but I will get back to you after a meeting tomorrow. |
|
Hi @josephzhang8. After meeting with @BenjaminJacob86 I realized I used scribedio schism output that was already partly processed. In the below function I have defined the preprocessing that is needed in order for our ugrid reader to read and visualize the raw schism scribedio output. This is a bit more than I communicated before, but I hope this provides useful information to enrich the schism output dataset. import xarray as xr
import numpy as np
from netCDF4 import default_fillvals
import xugrid as xu # pip install xugrid, adds ugrid and grid accessors to xarray datasets
def preprocess_schism_scribedio(ds):
"""
This preprocessing function describes the minimal changes that would have to be made
in the SCHISM output in order for it to be read via ugrid conventions with xugrid.
It is probably not a complete overview.
"""
# set varnames
gridname = "SCHISM_hgrid"
fnc_varn = f"{gridname}_face_nodes"
enc_varn = f"{gridname}_edge_nodes"
node_x = f"{gridname}_node_x"
node_y = f"{gridname}_node_y"
# set topology attributes to empty variable
topo_attrs = {"cf_role": "mesh_topology",
"topology_dimension": 2, # has to be int, not str
"node_coordinates": f"{node_x} {node_y}",
"face_node_connectivity": fnc_varn,
"edge_node_connectivity": enc_varn,
}
ds[gridname] = xr.DataArray(np.array(default_fillvals["i4"], dtype=np.int32), attrs=topo_attrs)
# assign necessary attributes to connectivity variables
fnc_attrs = {"_FillValue":-1, "start_index":1}
ds[fnc_varn] = ds[fnc_varn].assign_attrs(fnc_attrs)
ds[enc_varn] = ds[enc_varn].assign_attrs(fnc_attrs)
# set node_x/node_y as coordinate variables instead of data_vars
ds = ds.set_coords([node_x,node_y])
# to prevent xugrid UserWarning, but this is hardcoded and it should be different for lat/lon models.
# "UserWarning: No standard_name of ('projection_x_coordinate', 'longitude', 'projection_y_coordinate', 'latitude') in
# ['SCHISM_hgrid_node_x', 'SCHISM_hgrid_node_y']. Using SCHISM_hgrid_node_x and SCHISM_hgrid_node_y as projected x and y coordinates."
projected = True
if projected:
ds[node_x] = ds[node_x].assign_attrs({"standard_name":"projection_x_coordinate"})
ds[node_y] = ds[node_y].assign_attrs({"standard_name":"projection_y_coordinate"})
else:
ds[node_x] = ds[node_x].assign_attrs({"standard_name":"longitude"})
ds[node_y] = ds[node_y].assign_attrs({"standard_name":"latitude"})
# add variable with coordinate system, optional but convenient for loading into QGIS and other tools
# not yet properly read/updated by xugrid: https://github.com/Deltares/xugrid/issues/42
if projected:
grid_mapping_name = 'Unknown projected'
crs_varn = 'projected_coordinate_system'
crs_num = 25832 #UTM Zone 32N from communication with BJ
else:
grid_mapping_name = 'latitude_longitude'
crs_varn = 'wgs84'
crs_num = 4326
crs_str = f'EPSG:{crs_num}'
crs_attrs = {'epsg': crs_num, # epsg or EPSG_code are required for correct interpretation by QGIS
'EPSG_code': crs_str, # epsg or EPSG_code are required for correct interpretation by QGIS
'grid_mapping_name': grid_mapping_name,
}
ds[crs_varn] = xr.DataArray(np.array(default_fillvals['i4'],dtype=int),dims=(),attrs=crs_attrs)
# mesh attribute is required for d-ecoimpact
# valueable other attrs are "location" (node/face/edge),
# "standard_name", "long_name", "units", "grid_mapping"
for varn in ds.data_vars:
ds[varn] = ds[varn].assign_attrs({'mesh': gridname})
# time requires "units" attribute to be converted by xarray and other tools
# refdate taken from params.nml
ds['time'] = ds['time'].assign_attrs({"units":"seconds since 2017-01-02 00:00:00"})
ds = xr.decode_cf(ds)
#TODO: set _FillValue attribute for data_vars, test dataset did not seem to have nans
return ds
file_nc_pat = r"c:\Users\veenstra\Downloads\example\1_raw_outputs\*_1.nc"
# open the file with xarray, add ugrid conventions with preprocessing function and convert to UgridDataset
# in case of a pre-merged file this oneliner also works: `uds = xu.open_dataset(ds, preprocess=preprocess_schism_scribedio)`
# that would also work if the topology variables would be present in all datasets
ds = xr.open_mfdataset(file_nc_pat)
ds = preprocess_schism_scribedio(ds)
uds = xu.UgridDataset(ds)
# test with a ugrid plot to show the salinity variable was properly connected to the ugrid topology
uds.salinity.isel(time=-1, nSCHISM_vgrid_layers=0).ugrid.plot()It mostly comes down to creating an empty variable containing the mesh topology and adding some attributes. In the comments some other suggestions or the reasoning behind it is given. Other suggestions would be:
|
|
Thanks @veenstrajelmer for following up on this. We @josephzhang8 should in any case add the relevant (empty variable) metadata to any output. @veenstrajelmer Could you attach a We had a discussion previously on your suggestion
Our opinions differ on this issue, as this would require more storage space. @josephzhang8 would you reconsider? |
|
Regarding the comments:
|
|
Thanks for the responses. I recommend to not add the topology variable if the actual topology is not added, it will result in a corrupt file according to the ugrid conventions since the variables where the topology attributes point to are not present. The current setup is in that case better, so I would recommend to leave that as is. However, enriching the dataset with the topology variables inside is still beneficial. I do not have netcdf/ncdump installed outside of python, but used this python workaround. I have pasted the entire output of both files below. Original schism scribedio output (out2d_1.nc): Enriched file: |
|
Yes I believe there are very simple tools to add this info, e.g. NCL/NCO etc. |
|
Sure, but I guess the output is comparable to what I send. Is there any additional information you need from my side? Do you think it would be a valuable update to add the ugrid metadata to schism output? |
|
I'll let @platipodium add the missing meta in SCHISM code. As I explained, we need to be frugal about adding mesh info in all outputs. |
|
Yes, understandable, I am not sure if my response was clear. But I think in that case the current setup is good. Plain 3D variables in separate files, that are dependent on the out2d file in the same folder that contains all the connectivity/coordinates (and in the future also the topology metadata). Anyway, happy to be able to contribute something. |
|
Great; thx @veenstrajelmer @platipodium : can u plz add the missing metadata. Thx |
|
@platipodium Was this enhancement by any chance already implemented? I am asking this since we are in the continuing process of simplifying our processing tools and specifically making https://github.com/Deltares/D-EcoImpact support different model outputs. This SCHISM enhancement would be very helpful to make this process more transparent and simplify workflows for all D-EcoImpact users that want to work with SCHISM output. |
Sorry @veenstrajelmer, this has not yet been added. Feel free to do so :=) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Recently I recieved SCHISM model output, converted to ugrid format. We could not directly read this data with our generic ugrid reader, since there were some attributes missing or incorrect. I have written a preprocess function with which most things seem to work. Could you consider adding the adjustments to your ugrid writer/converter? I am not sure if this issue is the most convenient way to communicate this, but it might be a start.
The text was updated successfully, but these errors were encountered: