Interpolation error with meridional/zonal sections

[malmans2]

Hi Mattia,

I am extracting data along sections and it seems to be working for most sections but it fails when I try to take a section along a line of constant longitude. I get the following error:

20200226 141720.440 ERROR PET0 ESMF_FieldRegrid.F90:2570 ESMF_FieldRegridGetIwts Invalid argument - - can't currently regrid a grid that contains a DE of width less than 2

But if I try the same section but make the longitudes of the end points a little different it works. Is it really not possible to extract a section along a constant longitude? What am I missing? The bit of code is below.

lons = [2, 2]
lats = [66, 72]
section_4 = od.subsample.survey_stations(
     Xsurv = lons,
     Ysurv = lats,
     delta = delta_sec1,
     ZRange = [0, -700],
     timeRange = ["2018-01-01", "2018-02-01"],
     timeFreq = "6H",
     varList = ["Temp", "S", "U", "V", "drC", "drF", "HFacC", "HFacW",
"HFacS"],
     )
sec4_orth = section_4.compute.survey_aligned_velocities()
sec4_orth.to_netcdf('section4.nc',)

[malmans2]

To make things faster OceanSpy internally runs the cutout function before interpolating. However, the interpolation fails if the cutout grid does not have at least 2 lons and lats points. Because of that, OceanSpy guesses a pad to add to the grid. I'm not sure why in your case the OceanSpy guess is to small (something I should fix in the code at some point!).

Anyways, you can specify the pad to add to the cutout grid. It's the add_Hbdr argument. So add_Hbdr=1 will add/subtract 1 deg of longitude and latitude to the cutout. Keep it as small as possible. It's still running, but this is working for me:

lons = [2, 2]
lats = [66, 72]
section_4 = od.subsample.survey_stations(
     Xsurv = lons,
     Ysurv = lats,
     add_Hbdr=0.1,
     delta = 2,
     ZRange = [0, -700],
     timeRange = ["2018-01-01", "2018-02-01"],
     timeFreq = "6H",
     varList = ["Temp", "S", "U", "V", "drC", "drF", "HFacC", "HFacW",
"HFacS"],
     )

[malmans2]

Currently the horizontal pad for survey_stations is decided using add_Hbdr=True, which is an argument of cutout. The pad is based on the mean grid spacing. We should either rethink the way the pad is computed in cutout, or compute in a different way the pad in survey_stations: https://github.com/malmans2/oceanspy/blob/07391a8cdcceef3836926c28fc1f20d1bd6f2eb8/oceanspy/subsample.py#L191-L203

We always want to keep the pad as small as possible, but at the same time we don't want to trigger expensive computations to decide the pad.

[ThomasHaine]

@malmans2, what's the fix for this issue?

[ThomasHaine]

If the pad isn't big enough, expand it automatically?

[MaceKuailv]

I am testing it on ecco, there is a different issue coming up. Somehow HFacS and HFacW after the cutout become an array with both X and Xp1 in their dimension. @Mikejmnez. Could you take a quick look?

[MaceKuailv]

I have three proposals after looking into this.

1. Nicer error message. add a try...except structure when creating the xesmf regridder

2. give the pad a 1.5 factor. there is no guarantee that you can find two gridpoints in $2\Delta x$, but it is guaranteed when you change it to $2\Delta x$. Of course if the grid spacing is inhomogeneous, there is no guarantee you will find a solution.

3. use CKD tree to find the nearest indexes of the corner of the cutouts. If the corner index is at least 3 apart, then it would be enough for xesmf. The con is that building the tree for large dataset still takes quite a while (mostly because of reading XC and YC into memory). The code is also going to be a bit more complicated.

I suggest doing a combination of 1 and 2.

[ThomasHaine]

Sounds good to me. Any concerns with @MaceKuailv's proposal?

[Mikejmnez]

I am testing it on ecco, there is a different issue coming up. Somehow HFacS and HFacW after the cutout become an array with both X and Xp1 in their dimension. @Mikejmnez. Could you take a quick look?

Yes, I also saw that. At some point within llc_rearrange I defined a list of metric variables twice in two different parts of the code, one of them excluding the HFac variables which cause the weird behavior... It is being fixed in the new PR.

[MaceKuailv]

closed with PR #278