Visualizing Global Meshes #58
Comments
|
A pure-python implementation for dropping the elements that cross the IDL (i.e. Joseph's suggestion) is the following: def drop_elements_crossing_idl(ds: xr.Dataset, max_lon: float = 100) -> xr.Dataset:
"""
Drop triface elements crossing the International Date Line (IDL).
``max_lon`` is the maximum longitudinal "distance" in degrees for an element.
"""
if max_lon <= 0:
raise ValueError(f"Maximum longitudinal \"distance\" must be positive: {max_lon}")
a, b, c = ds.triface_nodes.T
indices_of_triface_nodes_crossing_idl = np.where(
(np.abs(ds.lon[a] - ds.lon[b]) > max_lon)
| (np.abs(ds.lon[a] - ds.lon[c]) > max_lon)
| (np.abs(ds.lon[b] - ds.lon[c]) > max_lon)
)[0]
ds = ds.drop_isel(triface=indices_of_triface_nodes_crossing_idl)
return dsIt can be used like this: from thalassa import api
ds = api.open_dataset("/path/to/output.nc")
ds = drop_elements_corssing_idl(ds=ds)
# use the normal api functions ...We could add this functionality in thalassa but is relatively slow for big meshes (30+ seconds on my laptop). I don't think it is something that we want to be using interactively. IMHV this belongs to the post-processing steps of a model. |
|
For the record, this is a faster, but also more memory intensive implementation (my laptop run out of memory when I run this on the STOFS output): def drop_elements_crossing_idl2(ds: xr.Dataset, max_lon: float = 100) -> xr.Dataset:
"""
Drop triface elements crossing the International Date Line (IDL).
``max_lon`` is the maximum longitudinal "distance" in degrees for an element.
"""
if max_lon <= 0:
raise ValueError(f"Maximum longitudinal \"distance\" must be positive: {max_lon}")
indices_of_triface_nodes_crossing_idl = np.where(
np.abs(ds.lon[ds.triface_nodes] - ds.lon[ds.triface_nodes.roll(three=1)]) > max_lon
)[0]
ds = ds.drop_isel(triface=indices_of_triface_nodes_crossing_idl)
return ds |
|
In the end, I added in utils.py a function that drops elements that cross the IDL: Lines 119 to 159 in 3dcc1b4 With the current implementation there are false positives, i.e. close to the poles, even elements that don't cross the IDL may be dropped. Nevertheless, it is quite fast and it should be good enough for a quick check of the results. Closing |
This is something that came up after investigating why the STOFS 2D Global model could not be visualized (#54). STOFS is an ADCIRC model, but the same problem should exist for SCHISM models, too, therefore I think that we should have a dedicated issue.
So, visualizing the STOFS model gives this output:
The problem is that there are elements crossing over the international time line and this messes up the visualization. The white line at ~65 latitude is because at that particular range the international line crosses the region North of Kamtchatka (not sure of the actual name).
This problem is something that we have encountered before while trying to visualize global meshes and pyposeidon does provide a fix for this and the docs do explain in a bit more detail. Applying the fix though takes significant time. So you need to do it a post-processing step. You can't do it dynamically in thalassa. At least not for these 12 million nodes files.
@brey can provide more details if necessary.
pinging @saeed-moghimi-noaa
The text was updated successfully, but these errors were encountered: