Incorrect SpatialAverager weights with big regions

[RondeauG]

When trying to use SpatialAverager to perform a global average, me and @aulemahal found some strange behaviour in the weight given to pixels.

These are the weights generated when trying to average over a 40° x 20° rectangle. This seems adequate.

Something funky happens for a 160° x 80° rectangle however, with curves appearing:

For a 320° x 160° region, the rectangle appears to "flip" and the weights are generated only with a slice of latitudes:

For the globe (360° x 180°), the weights are generated along a singular latitude.

Sample code:

eps = 1e-1
for (x, y) in [(20, 10), (40, 20), (80, 40), (160, 80), (180, 90)]:
    glob = Polygon([(-x + eps, -y + eps), (x - eps, -y + eps), (x - eps, y - eps), (-x + eps, y - eps)])

    savg = xe.SpatialAverager(ds, [glob])

    w = ds.tg_mean.isel(time=0).copy(data=savg.weights.data.todense().reshape((64, 128)))
    fig, ax = plt.subplots()
    w.plot(ax=ax)
    ax.set_title(f"{x} x {y}")

[huard]

I suspect this has to do with how you draw a line from one coordinate to the other. The default ESMF option is to use a great circle, see https://earthsystemmodeling.org/esmpy_doc/release/ESMF_8_2_0/html/LineType.html?highlight=great%20circle

I suspect this would go away using the option for cartesian lines, but this option is not exposed in xesmf for now as far as I know.

I suspect that's a fairly straightforward PR if you're interested to see this work out of the box. Otherwise, you can probably work around this by using a polygon with more segments along the contour you want to follow.

[huard]

See also: https://earthsystemmodeling.org/esmpy_doc/release/ESMF_8_2_0/html/api.html?highlight=180#great-circle-cells

[aulemahal]

AFAIU, "cartesian" lines are in the cartesian 3D space, not in a 2D latitude-longitude space (that wouldn't be cartesian anyway). It seems like ESMF does not offer an option to work in a purely lat-lon space.

I did test switching the line_type to the SpatialAverager regrid call and nothing changed. I'm not sure if this is my switch was ineffective or if the result is the same for both options. The latter seems fishy, it would mean I do not understand what's going on...

My guess is that we need to "densify" the polygons before converting them to meshes. I'm not sure of only increasing the number of nodes in the outer boundary is enough or if we also need to divide large polygons in pieces...

[huard]

Increasing the density of the polygon contour seems to do the trick. The figure below shows the weights for input grids of 1, 5 and 15 degree resolution and an uniform input field of value one.

The polygon for the top row is defined from the four corners: p1 = Polygon([(-80, -40), (80, -40), (80, 40), (-80, 40)])

The polygon for the bottom row is an interpolated version.

distance = np.linspace(0, 1, 100)
p2 = Polygon([p1.boundary.interpolate(d, normalized=True) for d in distance])

The white line is the polygon contour and the red line the great circle between polygon vertices.

I must say that I'm baffled by the behavior on the top row, the dependence on the resolution of the input grid is unexpected.

[RondeauG]

Very interesting find!

[huard]

v0.8 has a function to densify the polygon segments.