You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@climate.apache.org by go...@apache.org on 2013/07/25 00:45:01 UTC
svn commit: r1506769 -
/incubator/climate/branches/RefactorInput/ocw/dataset_processor.py
Author: goodale
Date: Wed Jul 24 22:45:01 2013
New Revision: 1506769
URL: http://svn.apache.org/r1506769
Log:
Progress CLIMATE-213:
- temporal_rebin(): Updated docstring and linted the function
- spatial_regrid(): Final code updates and linting
- General clean up of code using pylint as a guide
Modified:
incubator/climate/branches/RefactorInput/ocw/dataset_processor.py
Modified: incubator/climate/branches/RefactorInput/ocw/dataset_processor.py
URL: http://svn.apache.org/viewvc/incubator/climate/branches/RefactorInput/ocw/dataset_processor.py?rev=1506769&r1=1506768&r2=1506769&view=diff
==============================================================================
--- incubator/climate/branches/RefactorInput/ocw/dataset_processor.py (original)
+++ incubator/climate/branches/RefactorInput/ocw/dataset_processor.py Wed Jul 24 22:45:01 2013
@@ -23,11 +23,11 @@ import scipy.interpolate
import scipy.ndimage
from scipy.ndimage import map_coordinates
-def temporal_rebin(Dataset, temporal_resolution):
+def temporal_rebin(target_dataset, temporal_resolution):
""" Rebin a Dataset to a new temporal resolution
- :param Dataset: Dataset object that need temporal regridding applied
- :type Dataset: Open Climate Workbench Dataset Object
+ :param target_dataset: Dataset object that needs temporal regridding
+ :type target_dataset: Open Climate Workbench Dataset Object
:param temporal_resolution: The new temporal bin size
:type temporal_resolution: Python datetime.timedelta object
@@ -53,26 +53,49 @@ def spatial_regrid(target_dataset, new_l
:returns: A new spatially regridded Dataset
:rtype: Open Climate Workbench Dataset Object
"""
- # Create an empty Numpy Array of shape (len(Dataset.times), len(new_latitudes), len(new_longitudes))
-
- # Extract the Dataset.values, existing lats and lons
+ # Make masked array of shape (times, new_latitudes,new_longitudes)
+ new_values = ma.zeros([len(target_dataset.times),
+ len(new_latitudes),
+ len(new_longitudes)])
+
+ # Create grids of the given lats and lons for the underlying API
+ # NOTE: np.meshgrid() requires inputs (x, y) and returns data
+ # of shape(y|lat|rows, x|lon|columns). So we pass in lons, lats
+ # and get back data.shape(lats, lons)
+ lons, lats = np.meshgrid(target_dataset.lons, target_dataset.lats)
+ new_lons, new_lats = np.meshgrid(new_longitudes, new_latitudes)
+ # Convert all lats and lons into Numpy Masked Arrays
+ lats = ma.array(lats)
+ lons = ma.array(lons)
+ new_lats = ma.array(new_lats)
+ new_lons = ma.array(new_lons)
+ target_values = ma.array(target_dataset.values)
+
+ # Call _rcmes_spatial_regrid on each time slice
+ for i in range(len(target_dataset.times)):
+ new_values[i] = _rcmes_spatial_regrid(target_values[i],
+ lats,
+ lons,
+ new_lats,
+ new_lons)
- # Create Meshgrids of the existing lats and lons
-
- # Call _rcmes_spatial_regrid on each time slice
+ # TODO:
# This will call down to the _congrid() function and the lat and lon
# axis will be adjusted with the time axis being held constant
# Create a new Dataset Object to return using new data
- regridded_dataset = ds.Dataset(new_latitudes, new_longitudes, target_dataset.times, target_dataset.values)
+ regridded_dataset = ds.Dataset(new_latitudes,
+ new_longitudes,
+ target_dataset.times,
+ new_values)
return regridded_dataset
def ensemble(datasets):
pass
-def _rcmes_spatial_regrid(spatial_values, lat, lon, lat2, lon2, order=1, mdi=-999999999):
+def _rcmes_spatial_regrid(spatial_values, lat, lon, lat2, lon2, order=1):
'''
- Perform regridding from one set of lat,lon values onto a new set (lat2,lon2)
+ Spatial regrid from one set of lat,lon values onto a new set (lat2,lon2)
:param spatial_values: Values in a spatial grid that need to be regridded
:type spatial_values: 2d masked numpy array. shape (latitude, longitude)
@@ -86,18 +109,18 @@ def _rcmes_spatial_regrid(spatial_values
:type lon2: 2d numpy array. shape(latitudes, longitudes)
:param order: Interpolation order flag. 1=bi-linear, 3=cubic spline
:type order: [optional] Integer
- :param mdi: Fill value for missing data (used in created of masked array output)
- :type mdi: [optional] Integer *** NOT USED IN THIS FUNCTION ***
- :returns: Regridded 2d masked numpy array with a shape of (len(lat2), len(lon2))
+ :returns: 2d masked numpy array with shape(len(lat2), len(lon2))
:rtype: (float, float)
'''
nlat = spatial_values.shape[0]
nlon = spatial_values.shape[1]
+ #print nlat, nlon, "lats, lons - incoming dataset"
nlat2 = lat2.shape[0]
nlon2 = lon2.shape[1]
+ #print nlat2, nlon2, "NEW lats, lons - for the new grid output"
# To make our lives easier down the road, let's
# turn these into arrays of x & y coords
@@ -124,17 +147,22 @@ def _rcmes_spatial_regrid(spatial_values
# Deal with the "hard" break in the y-direction
lati = (nlat - 1) * (lati - lat.min()) / (lat.max() - lat.min())
- # Notes on dealing with MDI when regridding data.
- # Method adopted here:
- # Use bilinear interpolation of data by default (but user can specify other order using order=... in call)
- # Perform bilinear interpolation of data, and of mask.
- # To be conservative, new grid point which contained some missing data on the old grid is set to missing data.
- # -this is achieved by looking for any non-zero interpolated mask values.
- # To avoid issues with bilinear interpolation producing strong gradients leading into the MDI,
- # set values at MDI points to mean data value so little gradient visible = not ideal, but acceptable for now.
+ """
+ TODO: Review this docstring and see if it still holds true.
+ NOTE: This function doesn't use MDI currently. These are legacy comments
- # Set values in MDI so that similar to surroundings so don't produce large gradients when interpolating
- # Preserve MDI mask, by only changing data part of masked array object.
+ Notes on dealing with MDI when regridding data.
+ Method adopted here:
+ Use bilinear interpolation of data by default (but user can specify other order using order=... in call)
+ Perform bilinear interpolation of data, and of mask.
+ To be conservative, new grid point which contained some missing data on the old grid is set to missing data.
+ -this is achieved by looking for any non-zero interpolated mask values.
+ To avoid issues with bilinear interpolation producing strong gradients leading into the MDI,
+ set values at MDI points to mean data value so little gradient visible = not ideal, but acceptable for now.
+
+ Set values in MDI so that similar to surroundings so don't produce large gradients when interpolating
+ Preserve MDI mask, by only changing data part of masked array object.
+ """
for shift in (-1, 1):
for axis in (0, 1):
q_shifted = np.roll(spatial_values, shift=shift, axis=axis)
@@ -144,6 +172,7 @@ def _rcmes_spatial_regrid(spatial_values
# Now we actually interpolate
# map_coordinates does cubic interpolation by default,
# use "order=1" to preform bilinear interpolation instead...
+
regridded_values = map_coordinates(spatial_values, [lati, loni], order=order)
regridded_values = regridded_values.reshape([nlat2, nlon2])
@@ -244,9 +273,9 @@ def _congrid(a, newdims, method='linear'
return newa
elif method in ['spline']:
- oslices = [ slice(0,j) for j in old ]
+ oslices = [ slice(0, j) for j in old ]
oldcoords = np.ogrid[oslices]
- nslices = [ slice(0,j) for j in list(newdims) ]
+ nslices = [ slice(0, j) for j in list(newdims) ]
newcoords = np.mgrid[nslices]
newcoords_dims = range(np.rank(newcoords))