Gradient search resampling swath data gives transposed results

[mraspaud]

          Yikes, the resulting image has transposed!

Originally posted by @pnuu in #618 (comment)

[pnuu]

So, the issue is that resampling data that is a SwathDefinition with gradient search resample and specific projecktions results in images that are transposed. I've tested this with AVHRR AAPP L1b data via Satpy. Replacing the resampler with traditional nearest neighbour resampler gives images with correct orientation.

To reproduce:

import glob
from pyresample.geometry import AreaDefinition
from satpy import Scene

adef = AreaDefinition("EPSG_3035_4km", "EPSG:3035", "EPSG_3035_4km", "EPSG:3035", 967, 980, [2426378.0132, 1528101.2618, 6293974.6215, 5446513.5222])
# Use some day-time data that covers Europe
fnames = glob.glob("/data/hrpt_noaa19_20240826_1002_80154.l1b")
lcl = glbl.resample(adef, resampler='gradient_search')
# Optional for comparing with nn
# lcl = glbl.resample(adef)
lcl.save_datasets(base_dir='/tmp', tiled=True, blockxsize=512, blockysize=512, overviews=[], driver="COG")

Comparing the output of gdalinfo overview_20240826_100301.tif between the two resamplers shows they are identical, only the data are transposed with gradient search.

Switch to the built-in euro4 area and gradient search works. So something in the usage of Proj library, perhaps?

[pnuu]

So I thought it could be due to the projection (laea) is handled somehow weirdly, but the built-in ease_nh has the data in proper orientation. So why is my own area def behaving weirdly 🙈

[djhoese]

So a couple things that come to mind:

1. Does it have anything to do with the ascending/descending nature of the input swath data?
2. A while ago PROJ made it so axes in projection are passed in and returned in the order of the projection. For example, EPSG:4326 is actually defined as (lat, lon). Pyresample generally assumes it is getting (x, y) to/from pyproj. There is an always_xy=True keyword argument that can be provided to Transformer and Proj object creation that would fix this. I could imagine that if we are getting swapped axes (y, x) instead of (x, y) from PROJ that it could look like the data is being flipped, but they are actually in the wrong location. Just a guess.

[pnuu]

I need to check the ascending/descending orbit case. The first two I tested with behaved the same. But then I'd assume also the ease_nh area would show the data transposed as my "custom" area and it has the same projection (laea).