Merge pull request #225 from djhoese/feature-smart-radius

Add smarter default radius_of_influence to XArrayResamplerNN resampling

pyresample/geometry.py

@@ -475,6 +475,67 @@ def __str__(self):
  str(self.lons),
  str(self.lats))
 
+ def geocentric_resolution(self, ellps='WGS84', radius=None, nadir_factor=2):
+ """Calculate maximum geocentric pixel resolution.
+
+ If `lons` is a :class:`xarray.DataArray` object with a `resolution`
+ attribute, this will be used instead of loading the longitude and
+ latitude data. In this case the resolution attribute is assumed to
+ mean the nadir resolution of a swath and will be multiplied by the
+ `nadir_factor` to adjust for increases in the spatial resolution
+ towards the limb of the swath.
+
+ Args:
+ ellps (str): PROJ Ellipsoid for the Cartographic projection
+ used as the target geocentric coordinate reference system.
+ Default: 'WGS84'. Ignored if `radius` is provided.
+ radius (float): Spherical radius of the Earth to use instead of
+ the definitions in `ellps`.
+ nadir_factor (int): Number to multiply the nadir resolution
+ attribute by to reflect pixel size on the limb of the swath.
+
+ Returns: Estimated maximum pixel size in meters on a geocentric
+ coordinate system (X, Y, Z) representing the Earth.
+
+ Raises: RuntimeError if a simple search for valid longitude/latitude
+ data points found no valid data points.
+
+ """
+ if hasattr(self.lons, 'attrs') and 'resolution' in self.lons.attrs:
+ return self.lons.attrs['resolution'] * nadir_factor
+ if self.ndim == 1:
+ raise RuntimeError("Can't confidently determine geocentric "
+ "resolution for 1D swath.")
+ from pyproj import transform
+ rows = self.shape[0]
+ start_row = rows // 2 # middle row
+ src = Proj('+proj=latlong +datum=WGS84')
+ if radius:
+ dst = Proj("+proj=cart +a={} +b={}".format(radius, radius))
+ else:
+ dst = Proj("+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 hasattr(lons.data, 'compute'):
+ # dask arrays, compute them together
+ import dask.array as da
+ lons, lats = da.compute(lons, lats)
+ if hasattr(lons, 'values'):
+ # convert xarray to numpy array
+ lons = lons.values
+ lats = lats.values
+ lons = lons.ravel()
+ lats = lats.ravel()
+ alt = np.zeros_like(lons)
+
+ xyz = np.stack(transform(src, dst, lons, lats, alt), axis=1)
+ dist = np.linalg.norm(xyz[1] - xyz[0])
+ dist = dist[np.isfinite(dist)]
+ if not dist.size:
+ raise RuntimeError("Could not calculate geocentric resolution")
+ return dist[0]
+
 
 class GridDefinition(CoordinateDefinition):
  """Grid defined by lons and lats.
@@ -1936,6 +1997,35 @@ def __getitem__(self, key):
  self.crop_offset[1] + xslice.start)
  return new_area
 
+ def geocentric_resolution(self, ellps='WGS84', radius=None):
+ """Find best estimate for overall geocentric resolution."""
+ from pyproj import transform
+ rows, cols = self.shape
+ mid_row = rows // 2
+ mid_col = cols // 2
+ x, y = self.get_proj_vectors()
+ mid_col_x = np.repeat(x[mid_col], y.size)
+ mid_row_y = np.repeat(y[mid_row], x.size)
+ src = Proj(getattr(self, 'crs', self.proj_dict))
+ if radius:
+ dst = Proj("+proj=cart +a={} +b={}".format(radius, radius))
+ else:
+ dst = Proj("+proj=cart +ellps={}".format(ellps))
+ alt_x = np.zeros(x.size)
+ alt_y = np.zeros(y.size)
+ hor_xyz = np.stack(transform(src, dst, x, mid_row_y, alt_x), axis=1)
+ vert_xyz = np.stack(transform(src, dst, mid_col_x, y, alt_y), axis=1)
+ hor_dist = np.linalg.norm(np.diff(hor_xyz, axis=0), axis=1)
+ vert_dist = np.linalg.norm(np.diff(vert_xyz, axis=0), axis=1)
+ dist = np.concatenate((hor_dist, vert_dist))
+ dist = dist[np.isfinite(dist)]
+ if not dist.size:
+ raise RuntimeError("Could not calculate geocentric resolution")
+ # return np.max(dist) # alternative to histogram
+ # use the average of the largest histogram bin to avoid
+ # outliers and really large values
+ return np.mean(np.histogram_bin_edges(dist)[:2])
+
 
 def _make_slice_divisible(sli, max_size, factor=2):
  """Make the given slice even in size."""

pyresample/kd_tree.py

@@ -964,7 +964,7 @@ class XArrayResamplerNN(object):
  def __init__(self,
  source_geo_def,
  target_geo_def,
- radius_of_influence,
+ radius_of_influence=None,
  neighbours=1,
  epsilon=0):
  """
@@ -975,8 +975,10 @@ def __init__(self,
  Geometry definition of source
  target_geo_def : object
  Geometry definition of target
- radius_of_influence : float
- Cut off distance in meters
+ radius_of_influence : float, optional
+ Cut off distance in geocentric meters.
+ If not provided this will be estimated based on the source
+ and target geometry definition.
  neighbours : int, optional
  The number of neigbours to consider for each grid point.
  Default 1. Currently 1 is the only supported number.
@@ -997,10 +999,33 @@ def __init__(self,
  self.epsilon = epsilon
  self.source_geo_def = source_geo_def
  self.target_geo_def = target_geo_def
+ if radius_of_influence is None:
+ radius_of_influence = self._compute_radius_of_influence()
  self.radius_of_influence = radius_of_influence
  assert (self.target_geo_def.ndim == 2), \
  "Target area definition must be 2 dimensions"
 
+ def _compute_radius_of_influence(self):
+ """Estimate a good default radius_of_influence."""
+ try:
+ src_res = self.source_geo_def.geocentric_resolution()
+ except RuntimeError:
+ logger.warning("Could not calculate source definition resolution")
+ src_res = np.nan
+ try:
+ dst_res = self.target_geo_def.geocentric_resolution()
+ except RuntimeError:
+ logger.warning("Could not calculate destination definition "
+ "resolution")
+ dst_res = np.nan
+ radius_of_influence = np.nanmax([src_res, dst_res])
+ if np.isnan(radius_of_influence):
+ logger.warning("Could not calculate radius_of_influence, falling "
+ "back to 10000 meters. This may produce lower "
+ "quality results than expected.")
+ radius_of_influence = 10000
+ return radius_of_influence
+
  def _create_resample_kdtree(self, chunks=CHUNK_SIZE):
  """Set up kd tree on input"""
  source_lons, source_lats = self.source_geo_def.get_lonlats(

pyresample/test/test_geometry.py

@@ -25,8 +25,7 @@
 
 
 class Test(unittest.TestCase):
-
- """Unit testing the geometry and geo_filter modules"""
+ """Unit testing the geometry and geo_filter modules."""
 
  def test_lonlat_precomp(self):
  area_def = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD',
@@ -1222,6 +1221,30 @@ def test_get_lonlats_options(self):
  self.assertEqual(lon.dtype, np.dtype("f8"))
  self.assertIsInstance(lon, dask_array)
 
+ def test_area_def_geocentric_resolution(self):
+ """Test the AreaDefinition.geocentric_resolution method."""
+ from pyresample import get_area_def
+ area_extent = (-5570248.477339745, -5561247.267842293, 5567248.074173927, 5570248.477339745)
+ proj_dict = {'a': 6378169.0, 'b': 6356583.8, 'h': 35785831.0,
+ 'lon_0': 0.0, 'proj': 'geos', 'units': 'm'}
+ # metered projection
+ area_def = get_area_def('orig', 'Test area', 'test',
+ proj_dict,
+ 3712, 3712,
+ area_extent)
+ geo_res = area_def.geocentric_resolution()
+ np.testing.assert_allclose(10646.562531, geo_res)
+
+ # lon/lat
+ proj_dict = {'a': 6378169.0, 'b': 6356583.8, 'proj': 'latlong'}
+ area_def = get_area_def('orig', 'Test area', 'test',
+ proj_dict,
+ 3712, 3712,
+ [-130, 30, -120, 40],
+ area_extent)
+ geo_res = area_def.geocentric_resolution()
+ np.testing.assert_allclose(248.594116, geo_res)
+
 
 class TestMakeSliceDivisible(unittest.TestCase):
 
@@ -1264,7 +1287,6 @@ def assert_np_dict_allclose(dict1, dict2):
 
 
 class TestSwathDefinition(unittest.TestCase):
-
  """Test the SwathDefinition."""
 
  def test_swath(self):
@@ -1572,6 +1594,27 @@ def test_striding(self):
  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):
+ """Test the SwathDefinition.geocentric_resolution method."""
+ import dask.array as da
+ import xarray as xr
+ import numpy as np
+ from pyresample.geometry import SwathDefinition
+ 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 = SwathDefinition(xlons, xlats)
+ 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)
+
+ # 1D
+ xlats = xr.DataArray(da.from_array(lats.ravel(), chunks=2), dims=['y'])
+ xlons = xr.DataArray(da.from_array(lons.ravel(), chunks=2), dims=['y'])
+ sd = SwathDefinition(xlons, xlats)
+ self.assertRaises(RuntimeError, sd.geocentric_resolution)
 
 
 class TestStackedAreaDefinition(unittest.TestCase):

pyresample/test/test_kd_tree.py

@@ -9,10 +9,13 @@
 from pyresample import geometry, kd_tree, utils
 from pyresample.test.utils import catch_warnings
 
-if sys.version_info < (2, 7):
- import unittest2 as unittest
-else:
- import unittest
+import unittest
+
+try:
+ from unittest import mock
+except ImportError:
+ # python 2.7
+ import mock
 
 
 class Test(unittest.TestCase):
@@ -887,6 +890,47 @@ def test_nearest_area_2d_to_area_1n(self):
  expected = 27706753.0
  self.assertEqual(cross_sum, expected)
 
+ def test_nearest_area_2d_to_area_1n_no_roi(self):
+ """Test 2D area definition to 2D area definition; 1 neighbor, no radius of influence."""
+ from pyresample.kd_tree import XArrayResamplerNN
+ import xarray as xr
+ import dask.array as da
+ data = self.data_2d
+ resampler = XArrayResamplerNN(self.src_area_2d, self.area_def,
+ neighbours=1)
+ ninfo = resampler.get_neighbour_info()
+ for val in ninfo[:3]:
+ # vii, ia, voi
+ self.assertIsInstance(val, da.Array)
+ self.assertRaises(AssertionError,
+ resampler.get_sample_from_neighbour_info, data)
+
+ # rename data dimensions to match the expected area dimensions
+ data = data.rename({'my_dim_y': 'y', 'my_dim_x': 'x'})
+ res = resampler.get_sample_from_neighbour_info(data)
+ self.assertIsInstance(res, xr.DataArray)
+ self.assertIsInstance(res.data, da.Array)
+ res = res.values
+ cross_sum = np.nansum(res)
+ expected = 32114793.0
+ self.assertEqual(cross_sum, expected)
+
+ # pretend the resolutions can't be determined
+ with mock.patch.object(self.src_area_2d, 'geocentric_resolution') as sgr, \
+ mock.patch.object(self.area_def, 'geocentric_resolution') as dgr:
+ sgr.side_effect = RuntimeError
+ dgr.side_effect = RuntimeError
+ resampler = XArrayResamplerNN(self.src_area_2d, self.area_def,
+ neighbours=1)
+ resampler.get_neighbour_info()
+ res = resampler.get_sample_from_neighbour_info(data)
+ self.assertIsInstance(res, xr.DataArray)
+ self.assertIsInstance(res.data, da.Array)
+ res = res.values
+ cross_sum = np.nansum(res)
+ expected = 1855928.0
+ self.assertEqual(cross_sum, expected)
+
  def test_nearest_area_2d_to_area_1n_3d_data(self):
  """Test 2D area definition to 2D area definition; 1 neighbor, 3d data."""
  from pyresample.kd_tree import XArrayResamplerNN