Fix derivation of correct radius of influence when data layout is not standard

pyresample/geometry.py

@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
 #
-# Copyright (C) 2010-2020 Pyresample developers
+# Copyright (C) 2010-2023 Pyresample developers
 #
 # This program is free software: you can redistribute it and/or modify it under
 # the terms of the GNU Lesser General Public License as published by the Free
@@ -28,6 +28,7 @@
 from typing import Optional, Sequence, Union
 
 import numpy as np
+import xarray as xr
 import pyproj
 import yaml
 from pyproj import Geod, Proj
@@ -693,16 +694,29 @@ def geocentric_resolution(self, ellps='WGS84', radius=None, nadir_factor=2):
         if self.ndim == 1:
             raise RuntimeError("Can't confidently determine geocentric "
                                "resolution for 1D swath.")
-        rows = self.shape[0]
+
+        if isinstance(self.lons, xr.DataArray):
+            rows = self.lons['y'].shape[0]
+        else:
+            # Data have no information on dimensions, so we assume first dimension (the rows) is the y-axis:
+            logger.warning('As Numpy data arrays carry no information on the data layout we here ' +
+                           'assume the first dimension (the rows) is the y-axis (the satellite scans)')
+            rows = self.shape[0]
+
         start_row = rows // 2  # middle row
         src = CRS('+proj=latlong +datum=WGS84')
         if radius:
             dst = CRS("+proj=cart +a={} +b={}".format(radius, radius))
         else:
             dst = CRS("+proj=cart +ellps={}".format(ellps))
         # simply take the first two columns of the middle of the swath
-        lons = self.lons[start_row: start_row + 1, :2]
-        lats = self.lats[start_row: start_row + 1, :2]
+        if isinstance(self.lons, xr.DataArray):
+            lons = self.lons.sel(y=start_row)[:2]
+            lats = self.lats.sel(y=start_row)[:2]
+        else:
+            lons = self.lons[start_row: start_row + 1, :2]
+            lats = self.lats[start_row: start_row + 1, :2]
+
         if hasattr(lons.data, 'compute'):
             # dask arrays, compute them together
             import dask.array as da

pyresample/test/test_geometry/test_swath.py

@@ -1,4 +1,4 @@
-# Copyright (C) 2010-2022 Pyresample developers
+# Copyright (C) 2010-2023 Pyresample developers
 #
 # This program is free software: you can redistribute it and/or modify it under
 # the terms of the GNU Lesser General Public License as published by the Free
@@ -15,6 +15,7 @@
 """Test AreaDefinition objects."""
 import contextlib
 
+import logging
 import dask.array as da
 import numpy as np
 import pytest
@@ -512,13 +513,31 @@ def test_striding(self, create_test_swath):
         np.testing.assert_allclose(res.lons, [[178.5, -179.5]])
         np.testing.assert_allclose(res.lats, [[0, 0]], atol=2e-5)
 
-    def test_swath_def_geocentric_resolution(self, create_test_swath):
-        """Test the SwathDefinition.geocentric_resolution method."""
+    def test_swath_def_geocentric_resolution_numpy(self, caplog, create_test_swath):
+        """Test the SwathDefinition.geocentric_resolution method - lon/lat are a numpy arrays."""
+        lats = np.array([[0, 0, 0, 0], [1, 1, 1, 1.0]])
+        lons = np.array([[178.5, 179.5, -179.5, -178.5], [178.5, 179.5, -179.5, -178.5]])
+        sd = create_test_swath(lons, lats)
+
+        with caplog.at_level(logging.DEBUG):
+            geo_res = sd.geocentric_resolution()
+
+        log_output = ('As Numpy data arrays carry no information on the data layout we here assume ' +
+                      'the first dimension (the rows) is the y-axis (the satellite scans)')
+        assert log_output in caplog.text
+
+        # google says 1 degrees of longitude is about ~111.321km
+        # so this seems good
+        np.testing.assert_allclose(111301.237078, geo_res)
+
+    def test_swath_def_geocentric_resolution_xarray(self, create_test_swath):
+        """Test the SwathDefinition.geocentric_resolution method lon/lat are Xarray data arrays."""
         lats = np.array([[0, 0, 0, 0], [1, 1, 1, 1.0]])
         lons = np.array([[178.5, 179.5, -179.5, -178.5], [178.5, 179.5, -179.5, -178.5]])
         xlats = xr.DataArray(da.from_array(lats, chunks=2), dims=['y', 'x'])
         xlons = xr.DataArray(da.from_array(lons, chunks=2), dims=['y', 'x'])
         sd = create_test_swath(xlons, xlats)
+
         geo_res = sd.geocentric_resolution()
         # google says 1 degrees of longitude is about ~111.321km
         # so this seems good