Skip to content

Commit

Permalink
Make code still work if lon,lat arrays are Numpy arrays
Browse files Browse the repository at this point in the history
Signed-off-by: Adam.Dybbroe <[email protected]>
  • Loading branch information
Adam.Dybbroe committed Nov 20, 2023
1 parent 05d7e37 commit 17aabb9
Show file tree
Hide file tree
Showing 2 changed files with 38 additions and 6 deletions.
19 changes: 16 additions & 3 deletions pyresample/geometry.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -694,16 +695,28 @@ 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:
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.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
Expand Down
25 changes: 22 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 All @@ -15,6 +15,7 @@
"""Test AreaDefinition objects."""
import contextlib

import logging
import dask.array as da
import numpy as np
import pytest
Expand Down Expand Up @@ -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
Expand Down

0 comments on commit 17aabb9

Please sign in to comment.