climate-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From good...@apache.org
Subject svn commit: r1506604 - in /incubator/climate/branches/RefactorInput: ocw/dataset_processor.py rcmet/src/main/python/rcmes/toolkit/do_data_prep.py rcmet/src/main/python/rcmes/toolkit/process.py
Date Wed, 24 Jul 2013 16:01:08 GMT
Author: goodale
Date: Wed Jul 24 16:01:08 2013
New Revision: 1506604

URL: http://svn.apache.org/r1506604
Log:
Progress toward CLIMATE-213: Porting rcmes spatial regrid to ocw/dataset_processor

toolkit/do_data_prep.py - Moved a comment block closer to the functionality it references

toolkit/process.py - Moved the logic within do_regrid() over to the ocw/dataset_processor
and wired it up to call the new function in the new location.

ocw/dataset_processor.py - Added in the do_regrid() function as _rcmes_do_regrid()

Modified:
    incubator/climate/branches/RefactorInput/ocw/dataset_processor.py
    incubator/climate/branches/RefactorInput/rcmet/src/main/python/rcmes/toolkit/do_data_prep.py
    incubator/climate/branches/RefactorInput/rcmet/src/main/python/rcmes/toolkit/process.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=1506604&r1=1506603&r2=1506604&view=diff
==============================================================================
--- incubator/climate/branches/RefactorInput/ocw/dataset_processor.py (original)
+++ incubator/climate/branches/RefactorInput/ocw/dataset_processor.py Wed Jul 24 16:01:08
2013
@@ -16,8 +16,10 @@
 #
 
 import numpy as np
+import numpy.ma as ma
 import scipy.interpolate
 import scipy.ndimage
+from scipy.ndimage import map_coordinates
 
 def temporal_rebin(Dataset, temporal_resolution):
     """ Rebin a Dataset to a new temporal resolution
@@ -56,6 +58,95 @@ def spatial_regrid(Dataset, lon_resoluti
 def ensemble(datasets):
     pass
 
+def _rcmes_spatial_regrid(q, lat, lon, lat2, lon2, order=1, mdi=-999999999):
+    '''
+     Perform regridding from one set of lat,lon values onto a new set (lat2,lon2)
+    
+     Input::
+         q          - the variable to be regridded
+         lat,lon    - original co-ordinates corresponding to q values
+         lat2,lon2  - new set of latitudes and longitudes that you want to regrid q onto

+         order      - (optional) interpolation order 1=bi-linear, 3=cubic spline
+         mdi          - (optional) fill value for missing data (used in creation of masked
array)
+      
+     Output::
+         q2  - q regridded onto the new set of lat2,lon2 
+    
+    '''
+
+    nlat = q.shape[0]
+    nlon = q.shape[1]
+
+    nlat2 = lat2.shape[0]
+    nlon2 = lon2.shape[1]
+
+    # To make our lives easier down the road, let's 
+    # turn these into arrays of x & y coords
+    loni = lon2.ravel()
+    lati = lat2.ravel()
+
+    loni = loni.copy() # NB. it won't run unless you do this...
+    lati = lati.copy()
+
+    # Now, we'll set points outside the boundaries to lie along an edge
+    loni[loni > lon.max()] = lon.max()
+    loni[loni < lon.min()] = lon.min()
+    
+    # To deal with the "hard" break, we'll have to treat y differently,
+    # so we're just setting the min here...
+    lati[lati > lat.max()] = lat.max()
+    lati[lati < lat.min()] = lat.min()
+    
+    
+    # We need to convert these to (float) indicies
+    #   (xi should range from 0 to (nx - 1), etc)
+    loni = (nlon - 1) * (loni - lon.min()) / (lon.max() - lon.min())
+    
+    # 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.
+    
+    # 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(q, shift=shift, axis=axis)
+            idx = ~q_shifted.mask * q.mask
+            q.data[idx] = q_shifted[idx]
+
+    # Now we actually interpolate
+    # map_coordinates does cubic interpolation by default, 
+    # use "order=1" to preform bilinear interpolation instead...
+    q2 = map_coordinates(q, [lati, loni], order=order)
+    q2 = q2.reshape([nlat2, nlon2])
+
+    # Set values to missing data outside of original domain
+    q2 = ma.masked_array(q2, mask=np.logical_or(np.logical_or(lat2 >= lat.max(), 
+                                                              lat2 <= lat.min()), 
+                                                np.logical_or(lon2 <= lon.min(), 
+                                                              lon2 >= lon.max())))
+    
+    # Make second map using nearest neighbour interpolation -use this to determine locations
with MDI and mask these
+    qmdi = np.zeros_like(q)
+    qmdi[q.mask == True] = 1.
+    qmdi[q.mask == False] = 0.
+    qmdi_r = map_coordinates(qmdi, [lati, loni], order=order)
+    qmdi_r = qmdi_r.reshape([nlat2, nlon2])
+    mdimask = (qmdi_r != 0.0)
+    
+    # Combine missing data mask, with outside domain mask define above.
+    q2.mask = np.logical_or(mdimask, q2.mask)
+
+    return q2
+
 def _congrid(a, newdims, method='linear', centre=False, minusone=False):
     '''
     This function is from http://wiki.scipy.org/Cookbook/Rebinning - Example 3

Modified: incubator/climate/branches/RefactorInput/rcmet/src/main/python/rcmes/toolkit/do_data_prep.py
URL: http://svn.apache.org/viewvc/incubator/climate/branches/RefactorInput/rcmet/src/main/python/rcmes/toolkit/do_data_prep.py?rev=1506604&r1=1506603&r2=1506604&view=diff
==============================================================================
--- incubator/climate/branches/RefactorInput/rcmet/src/main/python/rcmes/toolkit/do_data_prep.py
(original)
+++ incubator/climate/branches/RefactorInput/rcmet/src/main/python/rcmes/toolkit/do_data_prep.py
Wed Jul 24 16:01:08 2013
@@ -138,8 +138,7 @@ def prep_data(settings, obsDatasetList, 
         print 'No input observation data file. EXIT'
         sys.exit()
 
-    ## Part 1: retrieve observation data from the database and regrid them
-    ##       NB. automatically uses local cache if already retrieved.
+    ## Part 0: Setup the regrid variables based on the regridOption given
 
     # preparation for spatial re-gridding: define the size of horizontal array of the target
interpolation grid system (ngrdX and ngrdY)
     print 'regridOption in prep_data= ', regridOption
@@ -180,6 +179,11 @@ def prep_data(settings, obsDatasetList, 
 
     regObsData = []
     
+    
+    
+    ## Part 1: retrieve observation data from the database and regrid them
+    ##       NB. automatically uses local cache if already retrieved.
+    
     for n in np.arange(numOBSs):
         # spatial regridding
         oLats, oLons, _, oTimes, oData = db.extractData(obsDatasetId[n],

Modified: incubator/climate/branches/RefactorInput/rcmet/src/main/python/rcmes/toolkit/process.py
URL: http://svn.apache.org/viewvc/incubator/climate/branches/RefactorInput/rcmet/src/main/python/rcmes/toolkit/process.py?rev=1506604&r1=1506603&r2=1506604&view=diff
==============================================================================
--- incubator/climate/branches/RefactorInput/rcmet/src/main/python/rcmes/toolkit/process.py
(original)
+++ incubator/climate/branches/RefactorInput/rcmet/src/main/python/rcmes/toolkit/process.py
Wed Jul 24 16:01:08 2013
@@ -29,7 +29,7 @@ import time
 import netCDF4
 import numpy as np
 import numpy.ma as ma
-from scipy.ndimage import map_coordinates
+
 
 
 def extract_subregion_from_data_array(data, lats, lons, latmin, latmax, lonmin, lonmax):
@@ -157,92 +157,12 @@ def calc_area_in_grid_box(latitude, dlat
 
     return A
 
-def do_regrid(q, lat, lon, lat2, lon2, order=1, mdi= -999999999):
-    '''
-     Perform regridding from one set of lat,lon values onto a new set (lat2,lon2)
-    
-     Input::
-         q          - the variable to be regridded
-         lat,lon    - original co-ordinates corresponding to q values
-         lat2,lon2  - new set of latitudes and longitudes that you want to regrid q onto

-         order      - (optional) interpolation order 1=bi-linear, 3=cubic spline
-         mdi  	    - (optional) fill value for missing data (used in creation of masked array)
-      
-     Output::
-         q2  - q regridded onto the new set of lat2,lon2 
-    
-    '''
-
-    nlat = q.shape[0]
-    nlon = q.shape[1]
-
-    nlat2 = lat2.shape[0]
-    nlon2 = lon2.shape[1]
-
-    # To make our lives easier down the road, let's 
-    # turn these into arrays of x & y coords
-    loni = lon2.ravel()
-    lati = lat2.ravel()
-
-    loni = loni.copy() # NB. it won't run unless you do this...
-    lati = lati.copy()
-
-    # Now, we'll set points outside the boundaries to lie along an edge
-    loni[loni > lon.max()] = lon.max()
-    loni[loni < lon.min()] = lon.min()
-    
-    # To deal with the "hard" break, we'll have to treat y differently,
-    # so we're just setting the min here...
-    lati[lati > lat.max()] = lat.max()
-    lati[lati < lat.min()] = lat.min()
-    
-    
-    # We need to convert these to (float) indicies
-    #   (xi should range from 0 to (nx - 1), etc)
-    loni = (nlon - 1) * (loni - lon.min()) / (lon.max() - lon.min())
-    
-    # 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.
-    
-    # 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(q, shift=shift, axis=axis)
-            idx = ~q_shifted.mask * q.mask
-            q.data[idx] = q_shifted[idx]
-
-    # Now we actually interpolate
-    # map_coordinates does cubic interpolation by default, 
-    # use "order=1" to preform bilinear interpolation instead...
-    q2 = map_coordinates(q, [lati, loni], order=order)
-    q2 = q2.reshape([nlat2, nlon2])
-
-    # Set values to missing data outside of original domain
-    q2 = ma.masked_array(q2, mask=np.logical_or(np.logical_or(lat2 >= lat.max(), 
-                                                              lat2 <= lat.min()), 
-                                                np.logical_or(lon2 <= lon.min(), 
-                                                              lon2 >= lon.max())))
-    
-    # Make second map using nearest neighbour interpolation -use this to determine locations
with MDI and mask these
-    qmdi = np.zeros_like(q)
-    qmdi[q.mask == True] = 1.
-    qmdi[q.mask == False] = 0.
-    qmdi_r = map_coordinates(qmdi, [lati, loni], order=order)
-    qmdi_r = qmdi_r.reshape([nlat2, nlon2])
-    mdimask = (qmdi_r != 0.0)
-    
-    # Combine missing data mask, with outside domain mask define above.
-    q2.mask = np.logical_or(mdimask, q2.mask)
+def do_regrid(q, lat, lon, lat2, lon2, order=1, mdi=-999999999):
+    """ 
+    This function has been moved to the ocw/dataset_processor module
+    """
+    from ocw import dataset_processor
+    q2 = dataset_processor._rcmes_spatial_regrid(q, lat, lon, lat2, lon2, order=1, mdi=-999999999)
 
     return q2
 



Mime
View raw message