Invalid Pole Cells for ne30x8 SCRIP Grid

[lantaosun]

Version

v2025.03.0

How did you install UXarray?

Conda

What happened?

I’m testing the Uxarray zonal mean function using an Arctic-refined SE grid (ne30x8). After calling the zonal mean function, the output contains 19 values corresponding to latitudes from -90° to 90° at 5-degree intervals. However, the zonal mean values for surface air temperature at both the North and South Poles are zero, which seems odd. At other latitudes, the values are all above 250 K, which are reasonable for a temperature field.

What did you expect to happen?

I was expecting values above 250 K. However, I understand that since the grid does not extend all the way to 90°N and 90°S, there may be no data at the poles, which could explain the zero values. One way to address this is to manually specify the latitude range—for example, using .zonal_mean(lat=(-85, 85, 5))—to avoid including the poles. That said, it might be helpful if the default behavior started slightly below 90°N/S to avoid this issue altogether.

Can you provide a MCVE to repoduce the bug?

# Recognition of unstructured grids and data handling
import uxarray as ux

# General Plotting
import cartopy.crs as ccrs

# Plotting with HoloViz
import holoviews as hv

# import hvplot.pandas
import geoviews.feature as gf

import numpy as np

file_dir = "/glade/derecho/scratch/lantao/Uxarray/"

grid_dir = "/glade/derecho/scratch/lantao/Uxarray/"

grid_file = "PAMIP_ne30x8_np4_SCRIP.nc"

data_file = "f.e22.F2000climo.ne30x8_mt12.preindustrial.100.cam.h1.0002-01.nc"

ds = ux.open_dataset(grid_dir + grid_file, file_dir + data_file)

zonal_mean_trefht = ds["TREFHT"][0].zonal_mean()

(zonal_mean_trefht.plot.line() * zonal_mean_trefht.plot.scatter(color="red")).opts(
    title="Zonal Average Plot (Default)", xticks=np.arange(-90, 100, 20), xlim=(-95, 95)
)

[rajeeja]

Thanks for sharing this bug. Is it possible to share actual files referenced in your bug report?

[philipc2]

Thanks for sharing this bug. Is it possible to share actual files referenced in your bug report?

The file is on glade, so I should be able to check it out if the permissions are open.

[lantaosun]

Thanks! I’ve moved the data to /glade/derecho/scratch/lantao/Uxarray/ and updated the script accordingly. Please let me know if you run into any issues with the data or the script.

[philipc2]

Hi @lantaosun

Looking at the grid structure, this is what it looks like over the poles.

9 total cells with one of them exactly centered about the pole (the stitching is a visualization artificial on how we handle periodic elements).

So for the case when the latitude is exactly +-90, we aren't correctly identifying the intersection. Looking further, this is what the latitude bounds look like for that cell.

[89.95580607, 89.96875017]

When it should be:

[89.95580607, 90.0]

@hongyuchen1030

I think somewhere along the way with our optimizations, we may have overlooked this edge case. I'm thinking that we should be able to modify our bounds implementation with the following to correctly handle this:

• If a face contains a pole point, set its maximum/minimum latitude to $\pm \frac{\pi}{2}$

[hongyuchen1030]

@hongyuchen1030

I think somewhere along the way with our optimizations, we may have overlooked this edge case. I'm thinking that we should be able to modify our bounds implementation with the following to correctly handle this:

• If a face contains a pole point, set its maximum/minimum latitude to
  ±

  π
  2
  

It should already be set this way in the latlon bound implementation

    if has_north_pole or has_south_pole:
        # Initial assumption that the pole point is inside the face
        is_center_pole = True

        pole_point_xyz = np.zeros(3, dtype=np.float64)
        pole_point_lonlat = np.zeros(2, dtype=np.float64)
        new_pt_latlon = np.zeros(2, dtype=np.float64)

        if has_north_pole:
            pole_point_xyz[2] = 1.0  # [0.0, 0.0, 1.0]
            pole_point_lonlat[1] = math.pi / 2  # [0.0, pi/2]
            new_pt_latlon[0] = math.pi / 2  # [pi/2, INT_FILL_VALUE]
            new_pt_latlon[1] = np.nan
        else:
            pole_point_xyz[2] = -1.0  # [0.0, 0.0, -1.0]
            pole_point_lonlat[1] = -math.pi / 2  # [0.0, -pi/2]
            new_pt_latlon[0] = -math.pi / 2  # [-pi/2, INT_FILL_VALUE]
            new_pt_latlon[1] = np.nan

So what actually happened is that this face doesn't actually contain the pole point. but is pretty close to the pole.

And when calculating with the zonal average, some floating point degeneracy occurs.

[hongyuchen1030]

@lantaosun

Can you also share which accuracy mode you're running with the zonal average? Are you using the robust weight calculation or not? Thanks

[lantaosun]

@lantaosun

Can you also share which accuracy mode you're running with the zonal average? Are you using the robust weight calculation or not? Thanks

Thanks for the question. Unfortunately, I’m not sure which accuracy mode is being used—my script (shared in the earlier post) uses the default settings. That being said, I’d be happy to run additional tests if you have any suggestions.

[philipc2]

@lantaosun

Can you also share which accuracy mode you're running with the zonal average? Are you using the robust weight calculation or not? Thanks

These are not exposed to the user

[hongyuchen1030]

@lantaosun
Can you also share which accuracy mode you're running with the zonal average? Are you using the robust weight calculation or not? Thanks

These are not exposed to the user

@philipc2
So my point is: this doesn't look like an implementation bug to me, but more like a floating point calculation issues. Since he's running the zonal_avg, do you know that's the default mode here?

[philipc2]

Some further digging into the face at the pole

# all for face index 209980

# face_node_connectivity
array([[6, 5, 7, 0, 3, 2, 1, 4, 4, 4]])

# node_lon 
array([-4.50000000e+01, -1.13686838e-13,  4.50000000e+01,  9.00000000e+01,
        1.35000000e+02, -1.80000000e+02, -1.35000000e+02, -9.00000000e+01,
       -9.00000000e+01, -9.00000000e+01])

# node_lat
array([89.95580607, 89.96875017, 89.95580607, 89.96875017, 89.95580607,
       89.96875017, 89.95580607, 89.96875017, 89.96875017, 89.96875017])


# double checking with the scrip "corner_lat" values
array([89.95580607, 89.96875017, 89.95580607, 89.96875017, 89.95580607,
       89.96875017, 89.95580607, 89.96875017, 89.96875017, 89.96875017])

# center is exactly at the pole 
array(90.)

Here is the single face projected over the north pole with an orthographic projection.

@hongyuchen1030

The geometry checks out here, so it must be either a precision issue or something that was lost in the translation from your previous implementation to the optimized one. Looking into it today.

[hongyuchen1030]

Some further digging into the face at the pole

all for face index 209980

face_node_connectivity

array([[6, 5, 7, 0, 3, 2, 1, 4, 4, 4]])

node_lon

array([-4.50000000e+01, -1.13686838e-13, 4.50000000e+01, 9.00000000e+01,
1.35000000e+02, -1.80000000e+02, -1.35000000e+02, -9.00000000e+01,
-9.00000000e+01, -9.00000000e+01])

node_lat

array([89.95580607, 89.96875017, 89.95580607, 89.96875017, 89.95580607,
89.96875017, 89.95580607, 89.96875017, 89.96875017, 89.96875017])

double checking with the scrip "corner_lat" values

array([89.95580607, 89.96875017, 89.95580607, 89.96875017, 89.95580607,
89.96875017, 89.95580607, 89.96875017, 89.96875017, 89.96875017])

center is exactly at the pole

array(90.)
Here is the single face projected over the north pole with an orthographic projection.

@hongyuchen1030

The geometry checks out here, so it must be either a precision issue or something that was lost in the translation from your previous implementation to the optimized one. Looking into it today.

So does the optimized one refer to the current implementation on the main branch, or to PR #1180? Because #1180 has not been merged yet.

Strangely, a similar case should already be covered by the existing test cases as well. I will also look into it today.

[hongyuchen1030]

Hi @philipc2

I just looked into this case and want to confirm a few points:

1. Repeated Node (Triple 4)
   In the face_node_connectivity array [[6, 5, 7, 0, 3, 2, 1, 4, 4, 4]], why is node index 4 repeated three times? Is this intentional or a data artifact? It seems unusual for a single face to reference the same node multiple times.

2. Convexity/Validity
   When visualizing the face by looking down from the north pole, the shape doesn’t appear to be strictly convex—possibly not even a typical “valid” face. Could you confirm whether this is expected? It might be that the face is meant to wrap around the pole, but it’s worth verifying.

3. Precision near the Pole
   As discussed earlier, querying points in the high-latitude range (89,90) often leads to precision issues. Because these nodes are only (10^{-4}) away from the north pole in Cartesian coordinates, tiny numerical differences can lead to large errors or unexpected geometry artifacts.

Let me know your thoughts on these points—especially the repeated node index and whether the face is indeed expected to be valid in its current form.