From 05d7e37623628c30ef3940fcc7a3cf067fd860a8 Mon Sep 17 00:00:00 2001 From: "Adam.Dybbroe" Date: Fri, 17 Nov 2023 16:24:15 +0100 Subject: [PATCH 1/3] Make use of the fact that longitudes are an xarray data array, and use 'y' dimension for rows/scanlines This makes the code resillient towards any other (awkward) data layout than the standard (where first dimension is usually the rows. Signed-off-by: Adam.Dybbroe --- pyresample/geometry.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/pyresample/geometry.py b/pyresample/geometry.py index 7586a06d8..5b8250d83 100644 --- a/pyresample/geometry.py +++ b/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 @@ -693,7 +693,8 @@ 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] + + rows = self.lons['y'].shape[0] start_row = rows // 2 # middle row src = CRS('+proj=latlong +datum=WGS84') if radius: @@ -701,8 +702,8 @@ def geocentric_resolution(self, ellps='WGS84', radius=None, nadir_factor=2): 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] + lons = self.lons.sel(y=start_row)[:2] + lats = self.lats.sel(y=start_row)[:2] if hasattr(lons.data, 'compute'): # dask arrays, compute them together import dask.array as da From 17aabb94cf0bc2c4398c66fb4b21c0b3f5ce59e6 Mon Sep 17 00:00:00 2001 From: "Adam.Dybbroe" Date: Mon, 20 Nov 2023 15:04:21 +0100 Subject: [PATCH 2/3] Make code still work if lon,lat arrays are Numpy arrays Signed-off-by: Adam.Dybbroe --- pyresample/geometry.py | 19 +++++++++++++--- pyresample/test/test_geometry/test_swath.py | 25 ++++++++++++++++++--- 2 files changed, 38 insertions(+), 6 deletions(-) diff --git a/pyresample/geometry.py b/pyresample/geometry.py index 5b8250d83..285f405c4 100644 --- a/pyresample/geometry.py +++ b/pyresample/geometry.py @@ -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 @@ -694,7 +695,14 @@ def geocentric_resolution(self, ellps='WGS84', radius=None, nadir_factor=2): raise RuntimeError("Can't confidently determine geocentric " "resolution for 1D swath.") - rows = self.lons['y'].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: @@ -702,8 +710,13 @@ def geocentric_resolution(self, ellps='WGS84', radius=None, nadir_factor=2): else: dst = CRS("+proj=cart +ellps={}".format(ellps)) # simply take the first two columns of the middle of the swath - lons = self.lons.sel(y=start_row)[:2] - lats = self.lats.sel(y=start_row)[: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 diff --git a/pyresample/test/test_geometry/test_swath.py b/pyresample/test/test_geometry/test_swath.py index 47a4a81ed..876923848 100644 --- a/pyresample/test/test_geometry/test_swath.py +++ b/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 From a7650f4c5725ccb3557ea381b3e8eb396d149259 Mon Sep 17 00:00:00 2001 From: "Adam.Dybbroe" Date: Tue, 21 Nov 2023 10:40:26 +0100 Subject: [PATCH 3/3] Make code work evenb without xarray available, and skip verbose warning Signed-off-by: Adam.Dybbroe --- pyresample/geometry.py | 10 ++++------ pyresample/test/test_geometry/test_swath.py | 11 ++--------- 2 files changed, 6 insertions(+), 15 deletions(-) diff --git a/pyresample/geometry.py b/pyresample/geometry.py index 285f405c4..c7b22618b 100644 --- a/pyresample/geometry.py +++ b/pyresample/geometry.py @@ -28,7 +28,6 @@ from typing import Optional, Sequence, Union import numpy as np -import xarray as xr import pyproj import yaml from pyproj import Geod, Proj @@ -695,12 +694,11 @@ def geocentric_resolution(self, ellps='WGS84', radius=None, nadir_factor=2): raise RuntimeError("Can't confidently determine geocentric " "resolution for 1D swath.") - if isinstance(self.lons, xr.DataArray): + if isinstance(self.lons, 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)') + # Data have no information on dimensions, so we assume first + # dimension (the rows) is the y-axis (the satellite scans): rows = self.shape[0] start_row = rows // 2 # middle row @@ -710,7 +708,7 @@ def geocentric_resolution(self, ellps='WGS84', radius=None, nadir_factor=2): else: dst = CRS("+proj=cart +ellps={}".format(ellps)) # simply take the first two columns of the middle of the swath - if isinstance(self.lons, xr.DataArray): + if isinstance(self.lons, DataArray): lons = self.lons.sel(y=start_row)[:2] lats = self.lats.sel(y=start_row)[:2] else: diff --git a/pyresample/test/test_geometry/test_swath.py b/pyresample/test/test_geometry/test_swath.py index 876923848..3a5014314 100644 --- a/pyresample/test/test_geometry/test_swath.py +++ b/pyresample/test/test_geometry/test_swath.py @@ -15,7 +15,6 @@ """Test AreaDefinition objects.""" import contextlib -import logging import dask.array as da import numpy as np import pytest @@ -513,18 +512,12 @@ 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_numpy(self, caplog, create_test_swath): + def test_swath_def_geocentric_resolution_numpy(self, 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 + geo_res = sd.geocentric_resolution() # google says 1 degrees of longitude is about ~111.321km # so this seems good