ocean-eddy-cpt / gcm-filters

Diffusion-based Spatial Filtering of Gridded Data
https://gcm-filters.readthedocs.io/
Other
37 stars 24 forks source link

ValueError while trying to filter with fixed factor, using MOM5 #162

Closed marimpacheco closed 3 months ago

marimpacheco commented 1 year ago

Hi,

I want to use gcm-filters to filter a lat/lon(3S-46S/132W-68W) section of the GFDL CM2.6 model (MOM5, 0.1°) using a filter with fixed factor (spatially-varying and anisotropic filter), and with a 15° longitude, 5° latitude cut-off wavelength. I want to filter tracers (temp, o2, ph).

Since MOM5 is NE-oriented, I selected the dxu, dxt, dyu and dyt accordingly so dxu would be located southwest of dxt (not sure if that is necessary? or a wrong way to proceed?).

Then I followed all the steps in the example and choosing as filter scale:

dx_max = max(dxu.max(),dyu.max(),dxt.max(),dyt.max())
dx_max = dx_max.values

filter_scale = 150 * dx_max # 15 degrees in longitude; each cell is 0.1 degree
kappa_w = dxt * dxt / (dx_max * dx_max)
kappa_s =  (1/900) * dyt * dyt / (dx_max * dx_max) # 5 degrees in latitude

filter_fixed_factor = gcm_filters.Filter(
    filter_scale=filter_scale,
    dx_min=dx_min,
    filter_shape=gcm_filters.FilterShape.GAUSSIAN,
    grid_type=gcm_filters.GridType.IRREGULAR_WITH_LAND,
    grid_vars={
        'wet_mask': wet_mask, 
        'dxw': dxt, 'dyw': dyt, 'dxs': dxu, 'dys': dyu, 'area': area, 
        'kappa_w': kappa_w, 'kappa_s': kappa_s
    }
)
filter_fixed_factor
Filter(filter_scale=1665563.4834155159, dx_min=array(7591.0808314), filter_shape=<FilterShape.GAUSSIAN: 1>, transition_width=3.141592653589793, ndim=2, n_steps=242, grid_type=<GridType.IRREGULAR_WITH_LAND: 5>)

However, when I try to run the filtering, I get the error: ValueError: At least one place in the domain must have either kappa_w = 1.0 or kappa_s = 1.Otherwise the filter's scale will not be equal to filter_scale anywhere in the domain.

But:

kappa_w.max()
1

If I add the (1/900) * to the kappa_s and remove it from kappa_w, it works.

How should I solve this? Thanks!

NoraLoose commented 1 year ago

Hi @marimpacheco, thanks for reaching out! I think I came across the same problem a while ago, which is why I opened PR #161. Let's merge the PR and see if this solves your problem.

marimpacheco commented 12 months ago

Hi Nora, I reinstalled gcm-filters, but unfortunately the problem remains

iangrooms commented 12 months ago

@marimpacheco did you install directly from the website? The changes have not been pushed to conda or pip yet, so if you re-installed from there you wouldn't have gotten the updates

marimpacheco commented 11 months ago

Ah right, from the commit it works, thanks!

marimpacheco commented 9 months ago

Hi, I have a question on how to adjust the dxw/dyw/dxs/dys values using MOM5 with the grid_type=gcm_filters.GridType.IRREGULAR_WITH_LAND, since I think I was taking the wrong approach.

In issue #24 @NoraLoose mentioned we should set: dxw = 0.5 * (np.roll(dxt, 1, axis=1) + np.roll(np.roll(dxt, 1, axis=1), 1, axis=0)), and similarly for dyw, dxs, dys (assuming that your coordinates have the order (j,i)). I assume the same equation applies to dyw, dyt, but I am not sure for dxs and dys. Is that something like: dxs = 0.5 * ( dxu + np.roll(dxu, -1, axis=0)) ? Thanks!

NoraLoose commented 9 months ago

Hi @marimpacheco!

If you have MOM5 data, I recommend using the grid types gcm_filters.GridType.MOM5U and gcm_filters.GridType.MOM5T. That way, you don't have to compute any of the grid data; instead you can just input the native grid information that comes with the model output.

Let me know if that is helpful, or if you still want to use the gcm_filters.GridType.IRREGULAR_WITH_LAND grid instead.

marimpacheco commented 9 months ago

I need to use the IRREGULAR_WITH_LAND in order to use a fixed factor filter. gcm_filters.GridType.MOM5T and gcm_filters.GridType.MOM5U do not have kappa_w and kappa_s. Or is there another way?

NoraLoose commented 9 months ago

True, in that case the IRREGULAR_WITH_LAND grid type is the one to go.

I recommend to use the python package xgcm to compute the necessary grid information.

import xgcm

coords = {
    "X": {"center": "xt_ocean", "right": "xu_ocean"},
    "Y": {"center": "yt_ocean", "right": "yu_ocean"}
}

metrics = {
    ('X',): ['dxu', 'dxt'], # X distances
    ('Y',): ['dyu', 'dyt'], # Y distances
}

grid = Grid(ds_grid, coords=coords, metrics=metrics, periodic=['X', 'Y'])

where ds_grid contains the grid variables of your MOM5 configuration.

Next, you can usegrid.interp to create dxw from dxt etc.

If you run into more issues, feel free to share a minimal example notebook. Hope this helps!

NoraLoose commented 9 months ago

Another (easier) option is to use a simple fixed factor filter with REGULAR_WITH_LAND_AREA_WEIGHTED. This grid type only requires grid variables that should be directly available from the MOM5 output.

https://gcm-filters.readthedocs.io/en/latest/examples/example_filter_types.html#simple-fixed-factor-filter

marimpacheco commented 9 months ago

Ah ok, thank you very much! I'll try with xgcm :)

iangrooms commented 4 months ago

@marimpacheco did xgcm resolve your problem? Can we close this issue?

marimpacheco commented 3 months ago

Yes, it did, thanks :)