Repository: climate Updated Branches: refs/heads/master 126b3b9a7 -> be6a84337
CLIMATE-714 - Updating the regridding routine - ocw.dataset_processor._rcmes_spatial_regrid has been replaced by scipy.interpolate.griddata to properly regrid datasets in curvilinear grid. Project: http://git-wip-us.apache.org/repos/asf/climate/repo Commit: http://git-wip-us.apache.org/repos/asf/climate/commit/6397dd98 Tree: http://git-wip-us.apache.org/repos/asf/climate/tree/6397dd98 Diff: http://git-wip-us.apache.org/repos/asf/climate/diff/6397dd98 Branch: refs/heads/master Commit: 6397dd985b8a3b94a459af96fcae752a10954e01 Parents: 6436186 Author: huikyole <[email protected]> Authored: Thu Dec 10 13:19:04 2015 -0800 Committer: huikyole <[email protected]> Committed: Thu Dec 10 13:19:04 2015 -0800 ---------------------------------------------------------------------- ocw/dataset_processor.py | 27 ++++++++++++++++++--------- 1 file changed, 18 insertions(+), 9 deletions(-) ---------------------------------------------------------------------- http://git-wip-us.apache.org/repos/asf/climate/blob/6397dd98/ocw/dataset_processor.py ---------------------------------------------------------------------- diff --git a/ocw/dataset_processor.py b/ocw/dataset_processor.py index 739ef5d..effa64c 100755 --- a/ocw/dataset_processor.py +++ b/ocw/dataset_processor.py @@ -209,17 +209,26 @@ def spatial_regrid(target_dataset, new_latitudes, new_longitudes): new_lats.shape[0], new_lons.shape[1]]) - # 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)): - print 'Regridding time = %d/%d' %(i+1,len(target_dataset.times)) - new_values[i] = scipy.interpolate.griddata((lons.flatten(), lats.flatten()), target_values[i].flatten(), + 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]]) # TODO:
