climate-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From good...@apache.org
Subject svn commit: r1506769 - /incubator/climate/branches/RefactorInput/ocw/dataset_processor.py
Date Wed, 24 Jul 2013 22:45:01 GMT
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))



Mime
View raw message