Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Open
wants to merge 3 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 16 additions & 4 deletions pyresample/geometry.py
Original file line number Diff line number Diff line change
@@ -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
Expand Down Expand Up @@ -693,16 +693,28 @@ 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, 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 (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, 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
Expand Down
18 changes: 15 additions & 3 deletions pyresample/test/test_geometry/test_swath.py
Original file line number Diff line number Diff line change
@@ -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
Expand Down Expand Up @@ -512,13 +512,25 @@ 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, 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)
geo_res = sd.geocentric_resolution()

# 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
Expand Down
Loading