You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@climate.apache.org by hu...@apache.org on 2016/01/25 02:50:02 UTC
[1/2] climate git commit: CLIMATE-722 - Extrapolation issues in
spatial_regrid
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 <hu...@argo.jpl.nasa.gov>
Authored: Fri Jan 22 20:12:03 2016 -0800
Committer: huikyole <hu...@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
[2/2] climate git commit: CLIMATE-723 - Extrapolation issues in
spatial_regrid
Posted by hu...@apache.org.
CLIMATE-723 - 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/c943988d
Tree: http://git-wip-us.apache.org/repos/asf/climate/tree/c943988d
Diff: http://git-wip-us.apache.org/repos/asf/climate/diff/c943988d
Branch: refs/heads/master
Commit: c943988d29220fb29b61f2822fb2f08417910d0f
Parents: 9c28fe6 2befe90
Author: huikyole <hu...@argo.jpl.nasa.gov>
Authored: Sun Jan 24 17:49:28 2016 -0800
Committer: huikyole <hu...@argo.jpl.nasa.gov>
Committed: Sun Jan 24 17:49:28 2016 -0800
----------------------------------------------------------------------
ocw/dataset_processor.py | 99 +++++++++++++++++++++++++++++++++----------
1 file changed, 77 insertions(+), 22 deletions(-)
----------------------------------------------------------------------
http://git-wip-us.apache.org/repos/asf/climate/blob/c943988d/ocw/dataset_processor.py
----------------------------------------------------------------------