EWA resampler fails (with METimage/VII data) with non-contiguous lonlat arrays

[ameraner]

The EWA resampler currently fails when trying to resample a METimage/VII swath.
To reproduce:

scn = Scene(filenames=filenames, reader='vii_l1b_nc')
scn.load(["vii_10690"])
scn_r = scn.resample("EPSG_4326_720x360", resampler='ewa', rows_per_scan=24)
scn_r.save_datasets(writer='simple_image')

results in

Traceback (most recent call last):
  File "/tcenas/home/andream/anaconda3/envs/devenv/lib/python3.11/site-packages/dask/base.py", line 662, in compute
    results = schedule(dsk, keys, **kwargs)
              ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/tcenas/home/andream/code/satpy_latest/pyresample/pyresample/ewa/dask_ewa.py", line 60, in _call_ll2cr
    swath_points_in_grid, cols, rows = ll2cr(new_src, target_geo_def)
                                       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/tcenas/home/andream/code/satpy_latest/pyresample/pyresample/ewa/ewa.py", line 78, in ll2cr
    swath_points_in_grid = _ll2cr.ll2cr_static(lons, lats, fill,
                           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "pyresample/ewa/_ll2cr.pyx", line 242, in pyresample.ewa._ll2cr.ll2cr_static
  File "<stringsource>", line 664, in View.MemoryView.memoryview_cwrapper
  File "<stringsource>", line 352, in View.MemoryView.memoryview.__cinit__
ValueError: ndarray is not C-contiguous

A simple fix seems to be adding the following lines after this block (where the lonlats are already computed ndarrays)

pyresample/pyresample/ewa/ewa.py

Lines 59 to 68 in 5d91b3a

	lons, lats = swath_def.get_lonlats()
	# ll2cr requires 64-bit floats due to pyproj limitations
	# also need a copy of lons, lats since they are written to in-place
	try:
	lons = lons.astype(np.float64, copy=copy)
	lats = lats.astype(np.float64, copy=copy)
	except TypeError:
	lons = lons.astype(np.float64)
	lats = lats.astype(np.float64)

    lons = np.ascontiguousarray(lons)
    lats = np.ascontiguousarray(lats)

Is this a good idea? Or is there a better way?

Fixing it at the source without breaking lazyness is probably not trivial, due to the lon-lats of METimage being computed with rather complex geotiepoints interpolations and coordinate transforms...

cc @cesclc

[djhoese]

That should probably be fine or the astype could maybe be changed with np.ascontiguousarray with the dtype passed to it as a kwarg. I'm not sure how that works. Oh...but yeah we do want it writable so maybe np.require should be used: https://numpy.org/doc/stable/reference/generated/numpy.require.html#numpy.require

Regardless, I'm not sure I understand why this is a thing. The dask array feeding into this should be rechunked to be aligned with the scan size. I've always assumed that the creation of the numpy array underneath would involve a np.concatenate meaning a new contiguous array. I'm not sure how it is possible that this isn't doing that here. Could I get access to the MetImage data you're using in this test case?

Side note: Can rows_per_scan be added to the .attrs of the data? Then it doesn't need to be specified in the resample call.

[ameraner]

Thanks for looking into this! I've sent you an example file.

Regarding the rows_per_scan: I'm currently adding it in my reader renaming PR: pytroll/satpy#3395

[djhoese]

Do I need a non-merged version of the reader to read this file? I'm getting a "No supported files found" error.

[djhoese]

Nevermind, copying/pasting the filename included the \ to escape the commas. Man these filenames.

[djhoese]

@ameraner I'm going to make a fix PR for this in a bit, but I'm noticing that numpy considers the lons/lats coming from the reader as being fortran (F) contiguous. Is that expected? That seems like something wrong in the reader.

[ameraner]

Thanks for the fix! I'm not sure, I don't know that part of code well - probably the issue is in geotiepoints somewhere (https://github.com/pytroll/python-geotiepoints/blob/main/geotiepoints/viiinterpolator.py), I see that other interpolators force contiguous arrays. Is it a problem or possibly just a performance penalty?

[djhoese]

Yeah likely a performance penalty both in generating the geolocation and using it. It would be better if that interpolation code pre-allocated the entire geolocation array and then filled it in rather than concatenating/appending. But it looks like the final step uses xarray to do the interp and has a comment saying it is faster than interpolating a 2D array so maybe this is xarray's fault...but actually I wonder if the two interp calls were swapped if the result would then be C contiguous.

Edit: Oh the way it is being called xarray is just calling numpy.interp.