"resample_nearest" GEO swath to eqc has tiny periodic "glitch" in results

[LSUESL]

Code Sample, a minimal, complete, and verifiable piece of code

The bulk of this code is just to set up the reprojection from CONUS GEO swath to "eqc" (Plate Carree). The resulting grid seems "flawed" (see plot).

Edit: Link to the data file; it's ~16MB.

from pyresample import get_area_def, kd_tree
from pyresample.geometry import SwathDefinition
from pyresample.kd_tree import resample_nearest

from netCDF4 import Dataset

import numpy as np
import math

import matplotlib.pyplot as plt

EARTH_RADIUS = 6371228.0
DEG2RAD = np.pi / 180.0
# For interpolation of reprojection (pyresample)
RADIUS_OF_INFLUENCE = 5000

def haversine_distance(pt1, pt2):
    lat1, lon1 = pt1
    lat2, lon2 = pt2
    dlat = DEG2RAD * (lat2-lat1)
    dlon = DEG2RAD * (lon2-lon1)
    haver1 = math.sin(dlat/2.0) * math.sin(dlat/2.0)
    haver2 = math.cos(DEG2RAD * lat1) * math.cos(DEG2RAD * lat2)
    haver3 = math.sin(dlon/2.0) * math.sin(dlon/2.0)
    a = haver1 + haver2 * haver3
    c = 2.0 * math.atan2(math.sqrt(a), math.sqrt(1.0 - a))
    return EARTH_RADIUS * c


# My data is just a composite of many CONUS scans from GOES-16
# Any GOES-16 CONUS netCDF that includes latitude/longitude variables could be substituted
ds = Dataset('20240301_J061_G16_C07_CONUS.nc')
lats = ds.variables['latitude']
lons = ds.variables['longitude']

# Specific details to GoM
row1 = 95
rown = row1 + 1500
col1 = 50
coln = col1 + 2000

W, E, S, N = (-99.0, -78.0, 15.0, 32.0)
geos_lon_0 = -75.0
resolution = 1000.

lat_0 = (S + N) / 2.0
lon_0 = (W + E) / 2.0

# Commented lines below result in values shown, not salient to issue
# total_x_extent = haversine_distance((lat_0, W), (lat_0, E))
# total_y_extent = haversine_distance((S, lon_0), (N, lon_0))
# x_extent = total_x_extent / (2.0 * math.cos(DEG2RAD * lat_0))
# y_extent = total_y_extent / 2.0
# area_extent = (-x_extent, -y_extent, x_extent, yes_extent)

area_extent = (-1166537.7983851556, -945190.700959653, 1166537.7983851556, 945190.700959653) 

# x_size = int(np.ceil(total_x_extent/resolution))
# y_size = int(np.ceil(total_y_extent/resolution))

x_size, y_size = 2140, 1891

# BEGIN HERE: setup area defs and resample

# Build EQC parameters
proj_dict = {'a': str(EARTH_RADIUS),
             'units': 'm',
             'proj': 'eqc',
             'lat_0': str(lat_0),
             'lon_0': str(lon_0),
             }

# Get area def for both projections
RectDef = get_area_def('eqc', 'eqc', 'eqc', proj_dict,
                       x_size, y_size, area_extent)
SwathDef = SwathDefinition(lats=lats[row1:rown, col1:coln],
                           lons=lons[row1:rown, col1:coln])

repro_lat = resample_nearest(SwathDef, lats[row1:rown, col1:coln], RectDef,
                                 radius_of_influence=RADIUS_OF_INFLUENCE,
                                 epsilon=1,
                                 )
# Plot any latitude "row"
plt.plot(repro_lat[500,:140])

Problem description

The plot of any (random) row of the resulting lats is shown.

The line oscillates over a small amplitude, but I would expect each row to be a constant latitude value. Larger values of "RADIUS_OF_INFLUENCE" and/or "epsilon" do not change the behavior, and its "periodic" nature seems "suspect"(?).

Or is this simply a limitation of the interpolation process, or do I misunderstand the "eqc" projection?

Despite being small values of variation, using the resulting lat/lon grid data gives perceivably "wobbly" gridline overlays, for example...

I appreciate your time and consideration!

Expected Output

Actual Result, Traceback if applicable

Versions of Python, package at hand and relevant dependencies

Python=3.9
pyresample=1.26.1
I tried Python 3.11, but it (in anaconda) pulls version 1.26.1 of pyresample.

[djhoese]

I'm very sorry for taking so long to answer this. I'm on paternity leave and had this issue open in a browser tab that got lost in the sea of tabs I've opened while checking emails. If I'm understanding your code and issue correctly, I think the main thing stems from your input longitude and latitudes. Have you checked them for consistency? I ask because it is a (well?) known issue that ABI L1b NetCDF files (C01-C16 for example) have 32-bit float scale factor and offset attributes for the X and Y coordinate variables. 32 bits is not enough precision for the distance covered by these coordinate variables. In Satpy (another pytroll project) I went to a lot of work to convert these scale and offset attributes to 64-bit floats and even round some of the decimal values when computing the geostationary X/Y metered coordinate arrays to get the theoretical/accurate values.

My guess (hope?) is that your lon/lats in your files were not generated with this level of care and are seeing the effects of the 32-bit scale factor attribute.