climate-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From huiky...@apache.org
Subject [1/2] climate git commit: CLIMATE-722 - Extrapolation issues in spatial_regrid
Date Mon, 25 Jan 2016 01:50:02 GMT
Repository: climate
Updated Branches:
  refs/heads/master 9c28fe6a9 -> c943988d2


CLIMATE-722 - Extrapolation issues in spatial_regrid

- ocw.dataset_processor.spatial_regrid has been updated.
- check if the new grid points are outside the original domain.
- calculate latitudes and longitudes (float) indices suitable for curvilinear grids


Project: http://git-wip-us.apache.org/repos/asf/climate/repo
Commit: http://git-wip-us.apache.org/repos/asf/climate/commit/2befe90c
Tree: http://git-wip-us.apache.org/repos/asf/climate/tree/2befe90c
Diff: http://git-wip-us.apache.org/repos/asf/climate/diff/2befe90c

Branch: refs/heads/master
Commit: 2befe90c2a056a541381ef8443b7d293fb5e38ea
Parents: d9e3c7e
Author: huikyole <huikyole@argo.jpl.nasa.gov>
Authored: Fri Jan 22 20:12:03 2016 -0800
Committer: huikyole <huikyole@argo.jpl.nasa.gov>
Committed: Fri Jan 22 20:12:03 2016 -0800

----------------------------------------------------------------------
 ocw/dataset_processor.py | 99 +++++++++++++++++++++++++++++++++----------
 1 file changed, 77 insertions(+), 22 deletions(-)
----------------------------------------------------------------------


http://git-wip-us.apache.org/repos/asf/climate/blob/2befe90c/ocw/dataset_processor.py
----------------------------------------------------------------------
diff --git a/ocw/dataset_processor.py b/ocw/dataset_processor.py
index 09edbaa..aca2a2f 100755
--- a/ocw/dataset_processor.py
+++ b/ocw/dataset_processor.py
@@ -22,8 +22,10 @@ import numpy as np
 import numpy.ma as ma
 import scipy.interpolate
 import scipy.ndimage
+from scipy.stats import rankdata
 from scipy.ndimage import map_coordinates
 import netCDF4
+from matplotlib.path import Path
 
 import logging
 
@@ -194,8 +196,10 @@ def spatial_regrid(target_dataset, new_latitudes, new_longitudes):
     #       of shape(y|lat|rows, x|lon|columns).  So we pass in lons, lats
     #       and get back data.shape(lats, lons)
     if target_dataset.lons.ndim ==1 and target_dataset.lats.ndim ==1:
+        regular_grid = True
         lons, lats = np.meshgrid(target_dataset.lons, target_dataset.lats)
     else:
+        regular_grid = False
         lons = target_dataset.lons
         lats = target_dataset.lats
     if new_longitudes.ndim ==1 and new_latitudes.ndim ==1:
@@ -204,32 +208,83 @@ def spatial_regrid(target_dataset, new_latitudes, new_longitudes):
         new_lons = new_longitudes
         new_lats = new_latitudes
 
+    ny_old, nx_old = lats.shape
+    ny_new, nx_new = new_lats.shape
+
     # Make masked array of shape (times, new_latitudes,new_longitudes)
     new_values = ma.zeros([len(target_dataset.times), 
-                           new_lats.shape[0],  
-                           new_lons.shape[1]])
-
-    # Call _rcmes_spatial_regrid on each time slice
+                           ny_new, nx_new]) 
+    # Make masked array of shape (times, new_latitudes,new_longitudes)
+    new_values = ma.zeros([len(target_dataset.times), 
+                           ny_new, nx_new])
+
+    # Boundary vertices of target_dataset
+    vertices = []
+
+    if regular_grid: 
+        vertices.append([lons[0,0], lats[0,0]])
+        vertices.append([lons[-1,0], lats[-1,0]])
+        vertices.append([lons[-1,-1], lats[-1,-1]])
+        vertices.append([lons[0,-1], lats[0,-1]])
+    else: 
+        for iy in np.arange(ny_old):   # from south to north along the west boundary
+            vertices.append([lons[iy,0], lats[iy,0]]) 
+        for ix in np.arange(nx_old):   # from west to east along the north boundary
+            vertices.append([lons[-1, ix], lats[-1, ix]])
+        for iy in np.arange(ny_old)[::-1]:   # from north to south along the east boundary
+            vertices.append([lons[iy, -1], lats[iy, -1]])
+        for ix in np.arange(nx_old)[::-1]:   # from east to west along the south boundary
+            vertices.append([lons[0, ix], lats[0, ix]])
+    path = Path(vertices)
+
+    # Convert new_lats and new_lons to float indices
+    new_lons_indices = np.zeros(new_lons.shape)
+    new_lats_indices = np.zeros(new_lats.shape)
+  
+    for iy in np.arange(ny_new):
+        for ix in np.arange(nx_new):
+            if path.contains_point([new_lons[iy,ix], new_lats[iy,ix]]): 
+                if regular_grid:
+                    new_lats_indices[iy,ix] = (ny_old -1.)*(new_lats[iy,ix] - lats.min())/(lats.max()
- lats.min())  
+                    new_lons_indices[iy,ix] = (nx_old -1.)*(new_lons[iy,ix] - lons.min())/(lons.max()
- lons.min())  
+                else:
+                    distance_from_original_grids = ((lons - new_lons[iy,ix])**2.+(lats -
new_lats[iy,ix])**2.)**0.5  
+                    if np.min(distance_from_original_grids) == 0.:
+                        new_lats_indices[iy,ix], new_lons_indices[iy,ix] = np.where(distance_from_original_grids
==0)
+                    else:
+                        distance_rank = rankdata(distance_from_original_grids.flatten(),
method='ordinal').reshape(lats.shape)                
+                        iy1, ix1 = np.where(distance_rank == 1) # the nearest grid point's
indices 
+                        iy2, ix2 = np.where(distance_rank == 4) # point [iy2, ix] is diagonally
across from [iy1, ix1] 
+                        dist1 = distance_from_original_grids[iy1,ix1]
+                        dist2 = distance_from_original_grids[iy2,ix2]
+                        new_lats_indices[iy,ix] = (dist1*iy2 + dist2*iy1)/(dist1+dist2) 
+                        new_lons_indices[iy,ix] = (dist1*ix2 + dist2*ix1)/(dist1+dist2) 
+            else:
+                new_lats_indices[iy,ix] = -999.    
+                new_lats_indices[iy,ix] = -999.    
+    new_lats_indices = ma.masked_less(new_lats_indices, 0.)
+    new_lons_indices = ma.masked_less(new_lons_indices, 0.)
+                      
+    # Regrid the data on each time slice
     for i in range(len(target_dataset.times)):
-        print 'Regridding time: %d/%d' %(i+1,len(target_dataset.times))
         values_original = ma.array(target_dataset.values[i])
-        if ma.count_masked(values_original) >= 1:
-            # Make a masking map using nearest neighbour interpolation -use this to determine
locations with MDI and mask these
-            qmdi = np.zeros_like(values_original)
-            qmdi[values_original.mask == True] = 1.
-            qmdi[values_original.mask == False] = 0.
-            qmdi_r = scipy.interpolate.griddata((lons.flatten(), lats.flatten()), qmdi.flatten(),
-                                              (new_lons.flatten(), new_lats.flatten()), method='nearest').reshape([new_lats.shape[0],new_lons.shape[1]])
-            mdimask = (qmdi_r != 0.0)
-
-            index = np.where(values_original.mask == False)
-            new_values[i] = scipy.interpolate.griddata((lons[index], lats[index]), values_original[index],
-                                              (new_lons.flatten(), new_lats.flatten()), method='linear',
fill_value=1.e+20).reshape([new_lats.shape[0],new_lons.shape[1]])
-            new_values[i] = ma.masked_greater(new_values[i], 1.e+19) 
-            new_values[i] = ma.array(new_values[i], mask = mdimask)
-        else:
-            new_values[i] = scipy.interpolate.griddata((lons.flatten(), lats.flatten()),
values_original.flatten(),
-                                              (new_lons.flatten(), new_lats.flatten()), method='linear').reshape([new_lats.shape[0],new_lons.shape[1]])
+        for shift in (-1, 1):
+            for axis in (0, 1):
+                q_shifted = np.roll(values_original, shift=shift, axis=axis)
+                idx = ~q_shifted.mask * values_original.mask
+                values_original.data[idx] = q_shifted[idx]
+        new_values[i] = map_coordinates(values_original, [new_lats_indices.flatten(), new_lons_indices.flatten()],
order=1).reshape(new_lats.shape) 
+        new_values[i] = ma.array(new_values[i], mask=new_lats_indices.mask)
+        # Make a masking map using nearest neighbour interpolation -use this to determine
locations with MDI and mask these
+        qmdi = np.zeros_like(values_original)
+        qmdi[values_original.mask == True] = 1.
+        qmdi[values_original.mask == False] = 0.
+        qmdi_r = map_coordinates(qmdi, [new_lats_indices.flatten(), new_lons_indices.flatten()],
order=1).reshape(new_lats.shape)
+        mdimask = (qmdi_r != 0.0)
+
+        # Combine missing data mask, with outside domain mask define above.
+        new_values[i].mask = np.logical_or(mdimask, new_values[i].mask)
+
     
     # TODO: 
     # This will call down to the _congrid() function and the lat and lon 


Mime
View raw message