[pyresample] 05/11: Imported Upstream version 1.1.0
Antonio Valentino
a_valentino-guest at moszumanska.debian.org
Sat Jun 21 09:18:19 UTC 2014
This is an automated email from the git hooks/post-receive script.
a_valentino-guest pushed a commit to branch master
in repository pyresample.
commit a2064955ae7be4c42373e1942b994a0862c20921
Author: Antonio Valentino <antonio.valentino at tiscali.it>
Date: Sat Jun 21 10:53:05 2014 +0200
Imported Upstream version 1.1.0
---
PKG-INFO | 2 +-
docs/source/_static/images/uncert_conc_nh.png | Bin 0 -> 193508 bytes
docs/source/_static/images/uncert_count_nh.png | Bin 0 -> 217710 bytes
docs/source/_static/images/uncert_stddev_nh.png | Bin 0 -> 201293 bytes
docs/source/conf.py | 25 +++
docs/source/plot.rst | 4 +-
docs/source/swath.rst | 37 +++-
pyresample.egg-info/PKG-INFO | 2 +-
pyresample.egg-info/SOURCES.txt | 3 +
pyresample/kd_tree.py | 230 ++++++++++++++++--------
pyresample/plot.py | 42 +----
pyresample/version.py | 2 +-
test/test_kd_tree.py | 144 ++++++++++++---
test/test_plot.py | 18 +-
14 files changed, 364 insertions(+), 145 deletions(-)
diff --git a/PKG-INFO b/PKG-INFO
index a8e3315..9e9ceba 100644
--- a/PKG-INFO
+++ b/PKG-INFO
@@ -1,6 +1,6 @@
Metadata-Version: 1.1
Name: pyresample
-Version: 1.0.0
+Version: 1.1.0
Summary: Resampling of remote sensing data in Python
Home-page: UNKNOWN
Author: Esben S. Nielsen
diff --git a/docs/source/_static/images/uncert_conc_nh.png b/docs/source/_static/images/uncert_conc_nh.png
new file mode 100644
index 0000000..56437d4
Binary files /dev/null and b/docs/source/_static/images/uncert_conc_nh.png differ
diff --git a/docs/source/_static/images/uncert_count_nh.png b/docs/source/_static/images/uncert_count_nh.png
new file mode 100644
index 0000000..90f505f
Binary files /dev/null and b/docs/source/_static/images/uncert_count_nh.png differ
diff --git a/docs/source/_static/images/uncert_stddev_nh.png b/docs/source/_static/images/uncert_stddev_nh.png
new file mode 100644
index 0000000..0b8f0d5
Binary files /dev/null and b/docs/source/_static/images/uncert_stddev_nh.png differ
diff --git a/docs/source/conf.py b/docs/source/conf.py
index d258b7f..8076a20 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -13,6 +13,31 @@
import sys, os
+class Mock(object):
+ def __init__(self, *args, **kwargs):
+ pass
+
+ def __call__(self, *args, **kwargs):
+ return Mock()
+
+ @classmethod
+ def __getattr__(cls, name):
+ if name in ('__file__', '__path__'):
+ return '/dev/null'
+ elif name[0] == name[0].upper():
+ mockType = type(name, (), {})
+ mockType.__module__ = __name__
+ return mockType
+ elif name == "inf":
+ return 0
+ else:
+ return Mock()
+
+MOCK_MODULES = ['numpy', 'pykdtree', 'configobj', 'pyproj',
+ 'scipy', 'scipy.spatial']
+for mod_name in MOCK_MODULES:
+ sys.modules[mod_name] = Mock()
+
# If extensions (or modules to document with autodoc) are in another directory,
# add these directories to sys.path here. If the directory is relative to the
# documentation root, use os.path.abspath to make it absolute, like shown here.
diff --git a/docs/source/plot.rst b/docs/source/plot.rst
index e808d2d..836998d 100644
--- a/docs/source/plot.rst
+++ b/docs/source/plot.rst
@@ -1,3 +1,5 @@
+.. _plot:
+
Plotting with pyresample and Basemap
====================================
Pyresample supports basic integration with Basemap (http://matplotlib.sourceforge.net/basemap).
@@ -133,7 +135,7 @@ Currently only the following set of Proj.4 arguments can be interpreted in the c
{'proj', 'a', 'b', 'ellps', 'lon_0', 'lat_0', 'lon_1', 'lat_1', 'lon_2', 'lat_2', 'lat_ts'}
Any other Proj.4 parameters will be ignored.
-If the ellipsoid is not defined in terms of 'ellps', 'a' or ('a', 'b') an exception will be raised.
+If the ellipsoid is not defined in terms of 'ellps', 'a' or ('a', 'b') it will default to WGS84.
The xsize and ysize in an AreaDefinition will only be used during resampling when the image data for use in
**basemap.imshow** is created. The actual size and shape of the final plot is handled by matplotlib.
diff --git a/docs/source/swath.rst b/docs/source/swath.rst
index cfa4bca..8359bce 100644
--- a/docs/source/swath.rst
+++ b/docs/source/swath.rst
@@ -29,10 +29,10 @@ The ImageContainerNearest class can be used for nearest neighbour resampling of
>>> area_con = swath_con.resample(area_def)
>>> result = area_con.image_data
-For other resampling types or splitting the process in two steps use the functions in **pyresample.swath** described below.
+For other resampling types or splitting the process in two steps use the functions in **pyresample.kd_tree** described below.
-pyresample.swath
-----------------
+pyresample.kd_tree
+------------------
This module contains several functions for resampling swath data.
@@ -164,6 +164,37 @@ If more channels are present in **data** the keyword argument **weight_funcs** m
If **data** is a masked array any pixel in the result data that has been "contaminated" by weighting of a masked pixel is masked.
+Uncertainty estimates
+*********************
+
+Uncertainty estimates in the form of weighted standard deviation can be obtained from the **resample_custom** and **resample_gauss** functions.
+By default the functions return the result of the resampling as a single numpy array. If the functions are given the keyword argument **with_uncert=True**
+then the following list of numpy arrays will be returned instead: **(result, stddev, count)**. **result** is the usual result. **stddev** is the weighted standard deviation for each element in the result. **count** is the number of data values used in the weighting for each element in the result.
+
+The principle is to view the calculated value for each element in the result as a weighted average of values sampled from a statistical variable.
+An estimate of the standard deviation of the distribution is calculated using the unbiased weighted estimator given as
+**stddev = sqrt((V1 / (V1 ** 2 + V2)) * sum(wi * (xi - result) ** 2))** where **result** is the result of the resampling. **xi** is the value of a contributing neighbour
+and **wi** is the corresponding weight. The coefficients are given as **V1 = sum(wi)** and **V2 = sum(wi ** 2)**. The standard deviation is only calculated for elements in
+the result where more than one neighbour has contributed to the weighting. The **count** numpy array can be used for filtering at a higher number of contributing neigbours.
+
+Usage only differs in the number of return values from **resample_gauss** and **resample_custom**. E.g.:
+
+ >>> result, stddev, count = pr.kd_tree.resample_gauss(swath_def, ice_conc, area_def,
+ ... radius_of_influence=20000,
+ ... sigmas=pr.utils.fwhm2sigma(35000),
+ ... fill_value=None, with_uncert=True)
+
+Below is shown a plot of the result of the resampling using a real data set:
+ .. image:: _static/images/uncert_conc_nh.png
+
+The corresponding standard deviations:
+ .. image:: _static/images/uncert_stddev_nh.png
+
+And the number of contributing neighbours for each element:
+ .. image:: _static/images/uncert_count_nh.png
+
+Notice the standard deviation is only calculated where there are more than one contributing neighbour.
+
Resampling from neighbour info
******************************
The resampling can be split in two steps:
diff --git a/pyresample.egg-info/PKG-INFO b/pyresample.egg-info/PKG-INFO
index a8e3315..9e9ceba 100644
--- a/pyresample.egg-info/PKG-INFO
+++ b/pyresample.egg-info/PKG-INFO
@@ -1,6 +1,6 @@
Metadata-Version: 1.1
Name: pyresample
-Version: 1.0.0
+Version: 1.1.0
Summary: Resampling of remote sensing data in Python
Home-page: UNKNOWN
Author: Esben S. Nielsen
diff --git a/pyresample.egg-info/SOURCES.txt b/pyresample.egg-info/SOURCES.txt
index 8f240ad..99a3d5a 100644
--- a/pyresample.egg-info/SOURCES.txt
+++ b/pyresample.egg-info/SOURCES.txt
@@ -20,6 +20,9 @@ docs/source/_static/images/tb37v_bmng.png
docs/source/_static/images/tb37v_ortho.png
docs/source/_static/images/tb37v_pc.png
docs/source/_static/images/tb37v_quick.png
+docs/source/_static/images/uncert_conc_nh.png
+docs/source/_static/images/uncert_count_nh.png
+docs/source/_static/images/uncert_stddev_nh.png
pyresample/__init__.py
pyresample/_multi_proc.py
pyresample/_spatial_mp.py
diff --git a/pyresample/kd_tree.py b/pyresample/kd_tree.py
index fc1027d..4302a52 100644
--- a/pyresample/kd_tree.py
+++ b/pyresample/kd_tree.py
@@ -89,15 +89,15 @@ def resample_nearest(source_geo_def, data, target_geo_def,
def resample_gauss(source_geo_def, data, target_geo_def,
radius_of_influence, sigmas, neighbours=8, epsilon=0,
- fill_value=0, reduce_data=True, nprocs=1, segments=None):
+ fill_value=0, reduce_data=True, nprocs=1, segments=None, with_uncert=False):
"""Resamples data using kd-tree gaussian weighting neighbour approach
:Parameters:
source_geo_def : object
Geometry definition of source
data : numpy array
- 1d array of single channel data points or
- (source_size, k) array of k channels of datapoints
+ Array of single channel data points or
+ (source_geo_def.shape, k) array of k channels of datapoints
target_geo_def : object
Geometry definition of target
radius_of_influence : float
@@ -123,14 +123,21 @@ def resample_gauss(source_geo_def, data, target_geo_def,
segments : {int, None}
Number of segments to use when resampling.
If set to None an estimate will be calculated
+ with_uncert : bool, optional
+ Calculate uncertainty estimates
- :Returns:
- data : numpy array
+ :Returns:
+ data : numpy array (default)
Source data resampled to target geometry
+
+ data, stddev, counts : numpy array, numpy array, numpy array (if with_uncert == True)
+ Source data resampled to target geometry.
+ Weighted standard devaition for all pixels having more than one source value
+ Counts of number of source values used in weighting per pixel
"""
-
+
def gauss(sigma):
- #Return gauss functino object
+ #Return gauss function object
return lambda r: np.exp(-r**2 / float(sigma)**2)
#Build correct sigma argument
@@ -156,20 +163,20 @@ def resample_gauss(source_geo_def, data, target_geo_def,
return _resample(source_geo_def, data, target_geo_def, 'custom',
radius_of_influence, neighbours=neighbours,
epsilon=epsilon, weight_funcs=weight_funcs, fill_value=fill_value,
- reduce_data=reduce_data, nprocs=nprocs, segments=segments)
+ reduce_data=reduce_data, nprocs=nprocs, segments=segments, with_uncert=with_uncert)
def resample_custom(source_geo_def, data, target_geo_def,
radius_of_influence, weight_funcs, neighbours=8,
epsilon=0, fill_value=0, reduce_data=True, nprocs=1,
- segments=None):
+ segments=None, with_uncert=False):
"""Resamples data using kd-tree custom radial weighting neighbour approach
:Parameters:
source_geo_def : object
Geometry definition of source
data : numpy array
- 1d array of single channel data points or
- (source_size, k) array of k channels of datapoints
+ Array of single channel data points or
+ (source_geo_def.shape, k) array of k channels of datapoints
target_geo_def : object
Geometry definition of target
radius_of_influence : float
@@ -197,10 +204,16 @@ def resample_custom(source_geo_def, data, target_geo_def,
Number of segments to use when resampling.
If set to None an estimate will be calculated
- :Returns:
- data : numpy array
+ :Returns:
+ data : numpy array (default)
Source data resampled to target geometry
+
+ data, stddev, counts : numpy array, numpy array, numpy array (if with_uncert == True)
+ Source data resampled to target geometry.
+ Weighted standard devaition for all pixels having more than one source value
+ Counts of number of source values used in weighting per pixel
"""
+
try:
for weight_func in weight_funcs:
if not isinstance(weight_func, types.FunctionType):
@@ -213,11 +226,11 @@ def resample_custom(source_geo_def, data, target_geo_def,
radius_of_influence, neighbours=neighbours,
epsilon=epsilon, weight_funcs=weight_funcs,
fill_value=fill_value, reduce_data=reduce_data,
- nprocs=nprocs, segments=segments)
+ nprocs=nprocs, segments=segments, with_uncert=with_uncert)
def _resample(source_geo_def, data, target_geo_def, resample_type,
radius_of_influence, neighbours=8, epsilon=0, weight_funcs=None,
- fill_value=0, reduce_data=True, nprocs=1, segments=None):
+ fill_value=0, reduce_data=True, nprocs=1, segments=None, with_uncert=False):
"""Resamples swath using kd-tree approach"""
valid_input_index, valid_output_index, index_array, distance_array = \
@@ -230,13 +243,16 @@ def _resample(source_geo_def, data, target_geo_def, resample_type,
nprocs=nprocs,
segments=segments)
+
return get_sample_from_neighbour_info(resample_type,
- target_geo_def.shape,
- data, valid_input_index,
- valid_output_index, index_array,
- distance_array=distance_array,
- weight_funcs=weight_funcs,
- fill_value=fill_value)
+ target_geo_def.shape,
+ data, valid_input_index,
+ valid_output_index,
+ index_array,
+ distance_array=distance_array,
+ weight_funcs=weight_funcs,
+ fill_value=fill_value,
+ with_uncert=with_uncert)
def get_neighbour_info(source_geo_def, target_geo_def, radius_of_influence,
neighbours=8, epsilon=0, reduce_data=True, nprocs=1, segments=None):
@@ -528,11 +544,11 @@ def _create_empty_info(source_geo_def, target_geo_def, neighbours):
return valid_output_index, index_array, distance_array
-
def get_sample_from_neighbour_info(resample_type, output_shape, data,
valid_input_index, valid_output_index,
index_array, distance_array=None,
- weight_funcs=None, fill_value=0):
+ weight_funcs=None, fill_value=0,
+ with_uncert=False):
"""Resamples swath based on neighbour info
:Parameters:
@@ -541,6 +557,8 @@ def get_sample_from_neighbour_info(resample_type, output_shape, data,
'custom': Resample based on weight_funcs
output_shape : (int, int)
Shape of output as (rows, cols)
+ data : numpy array
+ Source data
valid_input_index : numpy array
valid_input_index from get_neighbour_info
valid_output_index : numpy array
@@ -562,10 +580,10 @@ def get_sample_from_neighbour_info(resample_type, output_shape, data,
with undetermined pixels masked
:Returns:
- data : numpy array
+ result : numpy array
Source data resampled to target geometry
"""
-
+
if data.ndim > 2 and data.shape[0] * data.shape[1] == valid_input_index.size:
data = data.reshape(data.shape[0] * data.shape[1], data.shape[2])
elif data.shape[0] != valid_input_index.size:
@@ -575,16 +593,15 @@ def get_sample_from_neighbour_info(resample_type, output_shape, data,
raise ValueError('Mismatch between geometry and dataset')
is_multi_channel = (data.ndim > 1)
-
valid_input_size = valid_input_index.sum()
valid_output_size = valid_output_index.sum()
- if valid_input_size == 0 or valid_output_size == 0:
+ #Handle empty result set
+ if valid_input_size == 0 or valid_output_size == 0:
if is_multi_channel:
output_shape = list(output_shape)
output_shape.append(data.shape[1])
- #Handle empty result set
if fill_value is None:
#Use masked array for fill values
return np.ma.array(np.zeros(output_shape, data.dtype),
@@ -592,7 +609,7 @@ def get_sample_from_neighbour_info(resample_type, output_shape, data,
else:
#Return fill vaues for all pixels
return np.ones(output_shape, dtype=data.dtype) * fill_value
-
+
#Get size of output and reduced input
input_size = valid_input_size
if len(output_shape) > 1:
@@ -644,6 +661,11 @@ def get_sample_from_neighbour_info(resample_type, output_shape, data,
#Add the mask as channels to the dataset
is_masked_data = True
new_data = np.column_stack((new_data.data, new_data.mask))
+
+ if new_data.ndim > 1: # Multiple channels or masked input
+ output_shape = list(output_shape)
+ output_shape.append(new_data.shape[1])
+
#Prepare weight_funcs argument for handeling mask data
if weight_funcs is not None and is_masked_data:
@@ -666,76 +688,131 @@ def get_sample_from_neighbour_info(resample_type, output_shape, data,
result = new_data[new_index_array]
result[index_mask] = fill_value
else:
- #Calculate result using weighting
+ # Calculate result using weighting.
+ # Note: the code below has low readability in order
+ # to avoid looping over numpy arrays
#Get neighbours and masks of valid indices
ch_neighbour_list = []
index_mask_list = []
- for i in range(neighbours):
+ for i in range(neighbours): # Iterate over number of neighbours
+ # Make working copy neighbour index and
+ # set out of bounds indices to zero
index_ni = index_array[:, i].copy()
index_mask_ni = (index_ni == input_size)
index_ni[index_mask_ni] = 0
+
+ # Get channel data for the corresponing indices
ch_ni = new_data[index_ni]
ch_neighbour_list.append(ch_ni)
index_mask_list.append(index_mask_ni)
#Calculate weights
weight_list = []
- for i in range(neighbours):
- #Set out of bounds distance to 1 in order to avoid numerical Inf
+ for i in range(neighbours): # Iterate over number of neighbours
+ # Make working copy of neighbour distances and
+ # set out of bounds distance to 1 in order to avoid numerical Inf
distance = distance_array[:, i].copy()
distance[index_mask_list[i]] = 1
- if new_data.ndim > 1:
- #Calculate weights for each channel
- num_weights = valid_output_index.sum()
+ if new_data.ndim > 1: # More than one channel in data set.
+ # Calculate weights for each channel
weights = []
- for j in range(new_data.shape[1]):
+ num_weights = valid_output_index.sum()
+ num_channels = new_data.shape[1]
+ for j in range(num_channels):
calc_weight = weight_funcs[j](distance)
- #Use broadcasting to account for constant weight
+ # Turn a scalar weight into a numpy array
+ # (no effect if calc_weight already is an array)
expanded_calc_weight = np.ones(num_weights) * calc_weight
weights.append(expanded_calc_weight)
+
+ # Collect weights for all channels for neighbour number
weight_list.append(np.column_stack(weights))
- else:
+ else: # Only one channel
weights = weight_funcs(distance)
weight_list.append(weights)
result = 0
norm = 0
-
- #Calculate result
- for i in range(neighbours):
- #Find invalid indices to be masked of from calculation
- if new_data.ndim > 1:
+ count = 0
+ norm_sqr = 0
+ stddev = 0
+
+ # Calculate result
+ for i in range(neighbours): # Iterate over number of neighbours
+ # Find invalid indices to be masked of from calculation
+ if new_data.ndim > 1: # More than one channel in data set.
inv_index_mask = np.expand_dims(np.invert(index_mask_list[i]), axis=1)
- else:
+ else: # Only one channel
inv_index_mask = np.invert(index_mask_list[i])
#Aggregate result and norm
- result += inv_index_mask * ch_neighbour_list[i] * weight_list[i]
- norm += inv_index_mask * weight_list[i]
-
+ weights_tmp = inv_index_mask * weight_list[i]
+ result += weights_tmp * ch_neighbour_list[i]
+ norm += weights_tmp
+
#Normalize result and set fillvalue
- new_valid_index = (norm > 0)
- result[new_valid_index] /= norm[new_valid_index]
- result[np.invert(new_valid_index)] = fill_value
-
- #Add fill values
- if new_data.ndim > 1:
- full_result = np.ones((output_size, new_data.shape[1])) * fill_value
- else:
- full_result = np.ones(output_size) * fill_value
+ result_valid_index = (norm > 0)
+ result[result_valid_index] /= norm[result_valid_index]
+
+ if with_uncert: # Calculate uncertainties
+ # 2. pass to calculate standard deviation
+ for i in range(neighbours): # Iterate over number of neighbours
+ # Find invalid indices to be masked of from calculation
+ if new_data.ndim > 1: # More than one channel in data set.
+ inv_index_mask = np.expand_dims(np.invert(index_mask_list[i]), axis=1)
+ else: # Only one channel
+ inv_index_mask = np.invert(index_mask_list[i])
+
+ # Aggregate stddev information
+ weights_tmp = inv_index_mask * weight_list[i]
+ count += inv_index_mask
+ norm_sqr += weights_tmp ** 2
+ values = inv_index_mask * ch_neighbour_list[i]
+ stddev += weights_tmp * (values - result) ** 2
+
+ # Calculate final stddev
+ new_valid_index = (count > 1)
+ v1 = norm[new_valid_index]
+ v2 = norm_sqr[new_valid_index]
+ stddev[new_valid_index] = np.sqrt((v1 / (v1 ** 2 - v2)) * stddev[new_valid_index])
+ stddev[~new_valid_index] = np.NaN
+
+ #Add fill values
+ result[np.invert(result_valid_index)] = fill_value
+
+ # Create full result
+ if new_data.ndim > 1: # More than one channel
+ output_raw_shape = ((output_size, new_data.shape[1]))
+ else: # One channel
+ output_raw_shape = output_size
+
+ full_result = np.ones(output_raw_shape) * fill_value
full_result[valid_output_index] = result
result = full_result
-
- #Calculte correct output shape
- if new_data.ndim > 1:
- output_shape = list(output_shape)
- output_shape.append(new_data.shape[1])
-
+
+ if with_uncert: # Add fill values for uncertainty
+ full_stddev = np.ones(output_raw_shape) * np.nan
+ full_count = np.zeros(output_raw_shape)
+ full_stddev[valid_output_index] = stddev
+ full_count[valid_output_index] = count
+ stddev = full_stddev
+ count = full_count
+
+ stddev = stddev.reshape(output_shape)
+ count = count.reshape(output_shape)
+
+ if is_masked_data: # Ignore uncert computation of masks
+ stddev = _remask_data(stddev, is_to_be_masked=False)
+ count = _remask_data(count, is_to_be_masked=False)
+
+ # Set masks for invalid stddev
+ stddev = np.ma.array(stddev, mask=np.isnan(stddev))
+
#Reshape resampled data to correct shape
result = result.reshape(output_shape)
-
+
#Remap mask channels to create masked output
if is_masked_data:
result = _remask_data(result)
@@ -746,8 +823,15 @@ def get_sample_from_neighbour_info(resample_type, output_shape, data,
#Set output data type to input data type if relevant
if conserve_input_data_type:
- result = result.astype(input_data_type)
- return result
+ result = result.astype(input_data_type)
+
+ if with_uncert:
+ if np.ma.isMA(result):
+ stddev = np.ma.array(stddev, mask=(result.mask | stddev.mask))
+ count = np.ma.array(count, mask=result.mask)
+ return result, stddev, count
+ else:
+ return result
def _get_fill_mask_value(data_dtype):
"""Returns the maximum value of dtype"""
@@ -761,14 +845,18 @@ def _get_fill_mask_value(data_dtype):
data_dtype.type)
return fill_value
-def _remask_data(data):
+def _remask_data(data, is_to_be_masked=True):
"""Interprets half the array as mask for the other half"""
channels = data.shape[-1]
- mask = data[..., (channels // 2):]
- #All pixels affected by masked pixels are masked out
- mask = (mask != 0)
- data = np.ma.array(data[..., :(channels // 2)], mask=mask)
+ if is_to_be_masked:
+ mask = data[..., (channels // 2):]
+ #All pixels affected by masked pixels are masked out
+ mask = (mask != 0)
+ data = np.ma.array(data[..., :(channels // 2)], mask=mask)
+ else:
+ data = data[..., :(channels // 2)]
+
if data.shape[-1] == 1:
data = data.reshape(data.shape[:-1])
return data
diff --git a/pyresample/plot.py b/pyresample/plot.py
index f034015..86e1eed 100644
--- a/pyresample/plot.py
+++ b/pyresample/plot.py
@@ -114,40 +114,14 @@ def area_def2basemap(area_def, **kwargs):
basemap_args['projection'] = 'cyl'
else:
basemap_args['projection'] = area_def.proj_dict['proj']
- try:
- basemap_args['lon_0'] = float(area_def.proj_dict['lon_0'])
- except KeyError:
- pass
-
- try:
- basemap_args['lat_0'] = float(area_def.proj_dict['lat_0'])
- except KeyError:
- pass
-
- try:
- basemap_args['lon_1'] = float(area_def.proj_dict['lon_1'])
- except KeyError:
- pass
-
- try:
- basemap_args['lat_1'] = float(area_def.proj_dict['lat_1'])
- except KeyError:
- pass
-
- try:
- basemap_args['lon_2'] = float(area_def.proj_dict['lon_2'])
- except KeyError:
- pass
-
- try:
- basemap_args['lat_2'] = float(area_def.proj_dict['lat_2'])
- except KeyError:
- pass
-
- try:
- basemap_args['lat_ts'] = float(area_def.proj_dict['lat_ts'])
- except KeyError:
- pass
+
+ # Try adding potentially remaining args
+ for key in ('lon_0', 'lat_0', 'lon_1', 'lat_1', 'lon_2', 'lat_2',
+ 'lat_ts'):
+ try:
+ basemap_args[key] = float(area_def.proj_dict[key])
+ except KeyError:
+ pass
return Basemap(**basemap_args)
diff --git a/pyresample/version.py b/pyresample/version.py
index c77622e..a278cb1 100644
--- a/pyresample/version.py
+++ b/pyresample/version.py
@@ -15,4 +15,4 @@
#You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
-__version__ = '1.0.0'
+__version__ = '1.1.0'
diff --git a/test/test_kd_tree.py b/test/test_kd_tree.py
index 247b5ed..cdd5ef4 100644
--- a/test/test_kd_tree.py
+++ b/test/test_kd_tree.py
@@ -3,14 +3,16 @@ from __future__ import with_statement
import os
import sys
import unittest
+
import warnings
if sys.version_info < (2, 6):
warnings.simplefilter("ignore")
-else:
+else:
warnings.simplefilter("always")
+
import numpy
-from pyresample import kd_tree, utils, geometry, grid, data_reduce
+from pyresample import kd_tree, utils, geometry, data_reduce
def mp(f):
@@ -45,18 +47,17 @@ class Test(unittest.TestCase):
tlons = numpy.array([11.280789, 12.649354, 12.080402])
tlats = numpy.array([56.011037, 55.629675, 55.641535])
tswath = geometry.SwathDefinition(lons=tlons, lats=tlats)
- #grid = numpy.ones((1, 1, 2))
- #grid[0, 0, 0] = 12.562036
- #grid[0, 0, 1] = 55.715613
tgrid = geometry.CoordinateDefinition(lons=numpy.array([12.562036]),
lats=numpy.array([55.715613]))
-
+
+
def test_nearest_base(self):
res = kd_tree.resample_nearest(self.tswath,\
self.tdata.ravel(), self.tgrid,\
100000, reduce_data=False, segments=1)
self.assertTrue(res[0] == 2, 'Failed to calculate nearest neighbour')
+ @tmp
def test_gauss_base(self):
if sys.version_info < (2, 6):
res = kd_tree.resample_gauss(self.tswath, \
@@ -71,7 +72,8 @@ class Test(unittest.TestCase):
self.failIf(('Searching' not in str(w[0].message)), 'Failed to create correct neighbour warning')
self.assertAlmostEqual(res[0], 2.2020729, 5, \
'Failed to calculate gaussian weighting')
-
+
+
def test_custom_base(self):
def wf(dist):
return 1 - dist/100000.0
@@ -81,7 +83,7 @@ class Test(unittest.TestCase):
self.tdata.ravel(), self.tgrid,\
50000, wf, reduce_data=False, segments=1)
else:
- with warnings.catch_warnings(record=True) as w:
+ with warnings.catch_warnings(record=True) as w:
res = kd_tree.resample_custom(self.tswath,\
self.tdata.ravel(), self.tgrid,\
50000, wf, reduce_data=False, segments=1)
@@ -89,6 +91,56 @@ class Test(unittest.TestCase):
self.failIf(('Searching' not in str(w[0].message)), 'Failed to create correct neighbour warning')
self.assertAlmostEqual(res[0], 2.4356757, 5,\
'Failed to calculate custom weighting')
+
+ @tmp
+ def test_gauss_uncert(self):
+ sigma = utils.fwhm2sigma(41627.730557884883)
+ if sys.version_info < (2, 6):
+ res, stddev, count = kd_tree.resample_gauss(self.tswath, self.tdata,
+ self.tgrid, 100000, sigma,
+ with_uncert=True)
+ else:
+ with warnings.catch_warnings(record=True) as w:
+ res, stddev, count = kd_tree.resample_gauss(self.tswath, self.tdata,
+ self.tgrid, 100000, sigma,
+ with_uncert=True)
+ self.failIf(len(w) != 1, 'Failed to create neighbour warning')
+ self.failIf(('Searching' not in str(w[0].message)), 'Failed to create correct neighbour warning')
+
+ expected_res = 2.20206560694
+ expected_stddev = 0.707115076173
+ expected_count = 3
+ self.assertAlmostEqual(res[0], expected_res, 5, \
+ 'Failed to calculate gaussian weighting with uncertainty')
+ self.assertAlmostEqual(stddev[0], expected_stddev, 5, \
+ 'Failed to calculate uncertainty for gaussian weighting')
+ self.assertEqual(count[0], expected_count, 'Wrong data point count for gaussian weighting with uncertainty')
+
+ @tmp
+ def test_custom_uncert(self):
+ def wf(dist):
+ return 1 - dist/100000.0
+
+ if sys.version_info < (2, 6):
+ res, stddev, counts = kd_tree.resample_custom(self.tswath,
+ self.tdata, self.tgrid,
+ 100000, wf, with_uncert=True)
+ else:
+ with warnings.catch_warnings(record=True) as w:
+ res, stddev, counts = kd_tree.resample_custom(self.tswath,
+ self.tdata, self.tgrid,
+ 100000, wf, with_uncert=True)
+ self.failIf(len(w) != 1, 'Failed to create neighbour warning')
+ self.failIf(('Searching' not in str(w[0].message)), 'Failed to create correct neighbour warning')
+
+ self.assertAlmostEqual(res[0], 2.32193149, 5, \
+ 'Failed to calculate custom weighting with uncertainty')
+ self.assertAlmostEqual(stddev[0], 0.81817972, 5, \
+ 'Failed to calculate custom for gaussian weighting')
+ self.assertEqual(counts[0], 3, 'Wrong data point count for custom weighting with uncertainty')
+
+
+
def test_nearest(self):
data = numpy.fromfunction(lambda y, x: y*x, (50, 10))
lons = numpy.fromfunction(lambda y, x: 3 + x, (50, 10))
@@ -100,6 +152,7 @@ class Test(unittest.TestCase):
expected = 15874591.0
self.assertEqual(cross_sum, expected,\
msg='Swath resampling nearest failed')
+
def test_nearest_1d(self):
data = numpy.fromfunction(lambda x, y: x * y, (800, 800))
@@ -115,6 +168,7 @@ class Test(unittest.TestCase):
self.assertEqual(cross_sum, expected,
msg='Swath resampling nearest 1d failed')
+
def test_nearest_empty(self):
data = numpy.fromfunction(lambda y, x: y*x, (50, 10))
lons = numpy.fromfunction(lambda y, x: 165 + x, (50, 10))
@@ -127,6 +181,7 @@ class Test(unittest.TestCase):
self.assertEqual(cross_sum, expected,\
msg='Swath resampling nearest empty failed')
+
def test_nearest_empty_multi(self):
data = numpy.fromfunction(lambda y, x: y*x, (50, 10))
lons = numpy.fromfunction(lambda y, x: 165 + x, (50, 10))
@@ -139,6 +194,7 @@ class Test(unittest.TestCase):
self.assertEqual(res.shape, (800, 800, 3),\
msg='Swath resampling nearest empty multi failed')
+
def test_nearest_empty_multi_masked(self):
data = numpy.fromfunction(lambda y, x: y*x, (50, 10))
lons = numpy.fromfunction(lambda y, x: 165 + x, (50, 10))
@@ -151,7 +207,8 @@ class Test(unittest.TestCase):
fill_value=None)
self.assertEqual(res.shape, (800, 800, 3),
msg='Swath resampling nearest empty multi masked failed')
-
+
+
def test_nearest_empty_masked(self):
data = numpy.fromfunction(lambda y, x: y*x, (50, 10))
lons = numpy.fromfunction(lambda y, x: 165 + x, (50, 10))
@@ -165,6 +222,7 @@ class Test(unittest.TestCase):
self.assertTrue(cross_sum == expected,
msg='Swath resampling nearest empty masked failed')
+
def test_nearest_segments(self):
data = numpy.fromfunction(lambda y, x: y*x, (50, 10))
lons = numpy.fromfunction(lambda y, x: 3 + x, (50, 10))
@@ -177,6 +235,7 @@ class Test(unittest.TestCase):
self.assertEqual(cross_sum, expected,\
msg='Swath resampling nearest segments failed')
+
def test_nearest_remap(self):
data = numpy.fromfunction(lambda y, x: y*x, (50, 10))
lons = numpy.fromfunction(lambda y, x: 3 + x, (50, 10))
@@ -191,7 +250,7 @@ class Test(unittest.TestCase):
self.assertEqual(cross_sum, expected,\
msg='Grid remapping nearest failed')
- @mp
+
def test_nearest_mp(self):
data = numpy.fromfunction(lambda y, x: y*x, (50, 10))
lons = numpy.fromfunction(lambda y, x: 3 + x, (50, 10))
@@ -203,7 +262,8 @@ class Test(unittest.TestCase):
expected = 15874591.0
self.assertEqual(cross_sum, expected,\
msg='Swath resampling mp nearest failed')
-
+
+
def test_nearest_multi(self):
data = numpy.fromfunction(lambda y, x: y*x, (50, 10))
lons = numpy.fromfunction(lambda y, x: 3 + x, (50, 10))
@@ -217,7 +277,8 @@ class Test(unittest.TestCase):
expected = 3 * 15874591.0
self.assertEqual(cross_sum, expected,\
msg='Swath multi channel resampling nearest failed')
-
+
+
def test_nearest_multi_unraveled(self):
data = numpy.fromfunction(lambda y, x: y*x, (50, 10))
lons = numpy.fromfunction(lambda y, x: 3 + x, (50, 10))
@@ -230,7 +291,8 @@ class Test(unittest.TestCase):
expected = 3 * 15874591.0
self.assertEqual(cross_sum, expected,\
msg='Swath multi channel resampling nearest failed')
-
+
+
def test_gauss_sparse(self):
data = numpy.fromfunction(lambda y, x: y*x, (50, 10))
lons = numpy.fromfunction(lambda y, x: 3 + x, (50, 10))
@@ -242,7 +304,8 @@ class Test(unittest.TestCase):
expected = 15387753.9852
self.assertAlmostEqual(cross_sum, expected, places=3,\
msg='Swath gauss sparse nearest failed')
-
+
+
def test_gauss(self):
data = numpy.fromfunction(lambda y, x: (y + x)*10**-5, (5000, 100))
lons = numpy.fromfunction(lambda y, x: 3 + (10.0/100)*x, (5000, 100))
@@ -262,7 +325,7 @@ class Test(unittest.TestCase):
self.assertAlmostEqual(cross_sum, expected,\
msg='Swath resampling gauss failed')
- @tmp
+
def test_gauss_fwhm(self):
data = numpy.fromfunction(lambda y, x: (y + x)*10**-5, (5000, 100))
lons = numpy.fromfunction(lambda y, x: 3 + (10.0/100)*x, (5000, 100))
@@ -281,7 +344,8 @@ class Test(unittest.TestCase):
expected = 4872.81050892
self.assertAlmostEqual(cross_sum, expected,\
msg='Swath resampling gauss failed')
-
+
+
def test_gauss_multi(self):
data = numpy.fromfunction(lambda y, x: (y + x)*10**-6, (5000, 100))
lons = numpy.fromfunction(lambda y, x: 3 + (10.0/100)*x, (5000, 100))
@@ -302,6 +366,40 @@ class Test(unittest.TestCase):
expected = 1461.84313918
self.assertAlmostEqual(cross_sum, expected,\
msg='Swath multi channel resampling gauss failed')
+
+ @tmp
+ def test_gauss_multi_uncert(self):
+ data = numpy.fromfunction(lambda y, x: (y + x)*10**-6, (5000, 100))
+ lons = numpy.fromfunction(lambda y, x: 3 + (10.0/100)*x, (5000, 100))
+ lats = numpy.fromfunction(lambda y, x: 75 - (50.0/5000)*y, (5000, 100))
+ swath_def = geometry.SwathDefinition(lons=lons, lats=lats)
+ data_multi = numpy.column_stack((data.ravel(), data.ravel(),\
+ data.ravel()))
+ if sys.version_info < (2, 6):
+ res, stddev, counts = kd_tree.resample_gauss(swath_def, data_multi,\
+ self.area_def, 50000, [25000, 15000, 10000],
+ segments=1, with_uncert=True)
+ else:
+ with warnings.catch_warnings(record=True) as w:
+ res, stddev, counts = kd_tree.resample_gauss(swath_def, data_multi,\
+ self.area_def, 50000, [25000, 15000, 10000],
+ segments=1, with_uncert=True)
+ self.failIf(len(w) != 1, 'Failed to create neighbour radius warning')
+ self.failIf(('Possible more' not in str(w[0].message)), 'Failed to create correct neighbour radius warning')
+ cross_sum = res.sum()
+ cross_sum_stddev = stddev.sum()
+ cross_sum_counts = counts.sum()
+ expected = 1461.84313918
+ expected_stddev = 0.446204424799
+ expected_counts = 4934802.0
+ self.assertTrue(res.shape == stddev.shape and stddev.shape == counts.shape and counts.shape == (800, 800, 3))
+ self.assertAlmostEqual(cross_sum, expected,
+ msg='Swath multi channel resampling gauss failed on data')
+ self.assertAlmostEqual(cross_sum_stddev, expected_stddev,
+ msg='Swath multi channel resampling gauss failed on stddev')
+ self.assertAlmostEqual(cross_sum_counts, expected_counts,
+ msg='Swath multi channel resampling gauss failed on counts')
+
def test_gauss_multi_mp(self):
data = numpy.fromfunction(lambda y, x: (y + x)*10**-6, (5000, 100))
@@ -325,7 +423,8 @@ class Test(unittest.TestCase):
expected = 1461.84313918
self.assertAlmostEqual(cross_sum, expected,\
msg='Swath multi channel resampling gauss failed')
-
+
+
def test_gauss_multi_mp_segments(self):
data = numpy.fromfunction(lambda y, x: (y + x)*10**-6, (5000, 100))
lons = numpy.fromfunction(lambda y, x: 3 + (10.0/100)*x, (5000, 100))
@@ -348,7 +447,8 @@ class Test(unittest.TestCase):
expected = 1461.84313918
self.assertAlmostEqual(cross_sum, expected,\
msg='Swath multi channel segments resampling gauss failed')
-
+
+
def test_gauss_multi_mp_segments_empty(self):
data = numpy.fromfunction(lambda y, x: (y + x)*10**-6, (5000, 100))
lons = numpy.fromfunction(lambda y, x: 165 + (10.0/100)*x, (5000, 100))
@@ -364,6 +464,7 @@ class Test(unittest.TestCase):
msg=('Swath multi channel segments empty '
'resampling gauss failed'))
+
def test_custom(self):
def wf(dist):
return 1 - dist/100000.0
@@ -385,7 +486,8 @@ class Test(unittest.TestCase):
expected = 4872.81050729
self.assertAlmostEqual(cross_sum, expected,\
msg='Swath custom resampling failed')
-
+
+
def test_custom_multi(self):
def wf1(dist):
return 1 - dist/100000.0
@@ -642,6 +744,7 @@ class Test(unittest.TestCase):
self.assertEqual(cross_sum, expected,\
msg='Swath resampling from neighbour info nearest failed')
+
def test_custom_multi_from_sample(self):
def wf1(dist):
return 1 - dist/100000.0
@@ -712,9 +815,6 @@ class Test(unittest.TestCase):
lons = numpy.fromfunction(lambda y, x: 3 + x, (50, 10))
lats = numpy.fromfunction(lambda y, x: 75 - y, (50, 10))
swath_def = geometry.SwathDefinition(lons=lons, lats=lats)
-# res = swath.resample_nearest(lons.ravel(), lats.ravel(),
-# masked_data, self.area_def, 50000,
-# fill_value=None)
valid_input_index, valid_output_index, index_array, distance_array = \
kd_tree.get_neighbour_info(swath_def,
self.area_def,
diff --git a/test/test_plot.py b/test/test_plot.py
index a1580a5..d513f6e 100644
--- a/test/test_plot.py
+++ b/test/test_plot.py
@@ -5,6 +5,12 @@ import numpy as np
import pyresample as pr
+try:
+ import matplotlib
+ matplotlib.use('Agg')
+except ImportError:
+ pass # Postpone fail to individual tests
+
def tmp(f):
f.tmp = True
return f
@@ -17,8 +23,7 @@ class Test(unittest.TestCase):
lons = data[:, 0].astype(np.float64)
lats = data[:, 1].astype(np.float64)
tb37v = data[:, 2].astype(np.float64)
-
-
+
def test_ellps2axis(self):
a, b = pr.plot.ellps2axis('WGS84')
self.assertAlmostEqual(a, 6378137.0,
@@ -26,7 +31,6 @@ class Test(unittest.TestCase):
self.assertAlmostEqual(b, 6356752.3142451793,
msg='Failed to get semi-minor axis of ellipsis')
- @tmp
def test_area_def2basemap(self):
area_def = pr.utils.parse_area_file(os.path.join(os.path.dirname(__file__),
'test_files', 'areas.cfg'), 'ease_sh')[0]
@@ -35,10 +39,7 @@ class Test(unittest.TestCase):
bmap.rmajor == 6371228.0,
'Failed to create Basemap object')
-
def test_plate_carreeplot(self):
- import matplotlib
- matplotlib.use('Agg')
area_def = pr.utils.parse_area_file(os.path.join(os.path.dirname(__file__),
'test_files', 'areas.cfg'), 'pc_world')[0]
swath_def = pr.geometry.SwathDefinition(self.lons, self.lats)
@@ -47,10 +48,7 @@ class Test(unittest.TestCase):
fill_value=None)
plt = pr.plot._get_quicklook(area_def, result, num_meridians=0,
num_parallels=0)
-
def test_easeplot(self):
- import matplotlib
- matplotlib.use('Agg')
area_def = pr.utils.parse_area_file(os.path.join(os.path.dirname(__file__),
'test_files', 'areas.cfg'), 'ease_sh')[0]
swath_def = pr.geometry.SwathDefinition(self.lons, self.lats)
@@ -60,8 +58,6 @@ class Test(unittest.TestCase):
plt = pr.plot._get_quicklook(area_def, result)
def test_orthoplot(self):
- import matplotlib
- matplotlib.use('Agg')
area_def = pr.utils.parse_area_file(os.path.join(os.path.dirname(__file__),
'test_files', 'areas.cfg'), 'ortho')[0]
swath_def = pr.geometry.SwathDefinition(self.lons, self.lats)
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-grass/pyresample.git
More information about the Pkg-grass-devel
mailing list