The boundary_type = 'user' is supported with this commit.
Project: http://git-wip-us.apache.org/repos/asf/climate/repo Commit: http://git-wip-us.apache.org/repos/asf/climate/commit/0b753c7c Tree: http://git-wip-us.apache.org/repos/asf/climate/tree/0b753c7c Diff: http://git-wip-us.apache.org/repos/asf/climate/diff/0b753c7c Branch: refs/heads/master Commit: 0b753c7c09aac9322ee3296199e69362a55cc5fc Parents: 243a895 Author: huikyole <huiky...@argo.jpl.nasa.gov> Authored: Tue Aug 9 23:22:41 2016 -0700 Committer: huikyole <huiky...@argo.jpl.nasa.gov> Committed: Tue Aug 9 23:22:41 2016 -0700 ---------------------------------------------------------------------- ocw/dataset_processor.py | 15 ++++++++++-- ocw/utils.py | 55 +++++++++++++++++++++++++++++++++++++++---- 2 files changed, 64 insertions(+), 6 deletions(-) ---------------------------------------------------------------------- http://git-wip-us.apache.org/repos/asf/climate/blob/0b753c7c/ocw/dataset_processor.py ---------------------------------------------------------------------- diff --git a/ocw/dataset_processor.py b/ocw/dataset_processor.py index 56980c3..7379ef4 100755 --- a/ocw/dataset_processor.py +++ b/ocw/dataset_processor.py @@ -363,7 +363,7 @@ def ensemble(datasets): return ensemble_dataset -def subset(target_dataset, subregion, subregion_name=None, extract=True): +def subset(target_dataset, subregion, subregion_name=None, extract=True, user_mask_values=[1]): '''Subset given dataset(s) with subregion information :param subregion: The Bounds with which to subset the target Dataset. @@ -378,6 +378,9 @@ def subset(target_dataset, subregion, subregion_name=None, extract=True): :param extract: If False, the dataset inside regions will be masked. :type extract: :mod:`boolean` + :param user_mask_value: grid points where mask_variable == user_mask_value will be extracted or masked . + :type user_mask_value: :mod:`int` + :returns: The subset-ed Dataset object :rtype: :class:`dataset.Dataset` @@ -465,10 +468,18 @@ def subset(target_dataset, subregion, subregion_name=None, extract=True): target_dataset = temporal_slice( target_dataset, start_time_index, end_time_index) nt, ny, nx = target_dataset.values.shape - spatial_mask = utils.mask_using_shapefile_info(target_dataset.lons, target_dataset.lats, subregion.masked_regions, extract = extract) + spatial_mask = utils.mask_using_shapefile_info(target_dataset.lons, target_dataset.lats, + subregion.masked_regions, extract = extract) target_dataset.values = utils.propagate_spatial_mask_over_time(target_dataset.values, mask=spatial_mask) return target_dataset + if subregion.boundary_type == 'user': + spatial_mask = utils.regrid_spatial_mask(target_dataset.lons, target_dataset.lats, + subregion.mask_longitude, subregion.mask_latitude, subregion.mask_variable, + user_mask_values, extract = extract) + target_dataset.values = utils.propagate_spatial_mask_over_time(target_dataset.values, mask=spatial_mask) + return target_dataset + def temporal_slice(target_dataset, start_time_index, end_time_index): '''Temporally slice given dataset(s) with subregion information. This does not spatially subset the target_Dataset http://git-wip-us.apache.org/repos/asf/climate/blob/0b753c7c/ocw/utils.py ---------------------------------------------------------------------- diff --git a/ocw/utils.py b/ocw/utils.py index 0a9ba22..1820786 100755 --- a/ocw/utils.py +++ b/ocw/utils.py @@ -23,10 +23,12 @@ import datetime as dt import numpy as np import numpy.ma as ma -from mpl_toolkits.basemap import shiftgrid, Basemap, maskoceans +from mpl_toolkits.basemap import shiftgrid, Basemap from matplotlib.path import Path from dateutil.relativedelta import relativedelta from netCDF4 import num2date +import scipy.interpolate as interpolate +from scipy.ndimage import map_coordinates def decode_time_values(dataset, time_var_name): @@ -522,16 +524,61 @@ def mask_using_shapefile_info(lons, lats, masked_regions, extract = True): if iregion == 0: mask = mask0 else: - mask = mask|mask0 + mask = mask | mask0 if extract: mask = np.invert(mask) return mask +def regrid_spatial_mask(target_lon, target_lat, mask_lon, mask_lat, mask_var, + user_mask_values, extract=True): + target_lons, target_lats = convert_lat_lon_2d_array(target_lon, target_lat) + mask_lons, mask_lats = convert_lat_lon_2d_array(mask_lon, mask_lat) + + if target_lons != mask_lons or target_lats != mask_lats: + mask_var_regridded = interpolate.griddata((mask_lons.flatten(), mask_lats.flatten()), + mask_var.flatten(), + (target_lons.flatten(), target_lats.flatten()), + method='nearest', + fill_value=-9999.).reshape(target_lons.shape) + else: + mask_var_regridded = mask_var + + mask_outside = ma.masked_equal(mask_var_regridded, -9999.).mask + values_original = ma.array(mask_var) + 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] + # 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, [target_lats.flatten( + ), target_lons.flatten()], order=1).reshape(target_lons.shape) + mdimask = (qmdi_r != 0.0) + + for value in user_mask_values: + mask_var_regridded = ma.masked_equal(mask_var_regridded, value) + + if extract: + mask_out = np.invert(mask_var_regridded.mask | mdimask) + return mask_out | mask_outside + else: + mask_out = mask_var_regridded.mask | mdimask + return mask_out | mask_outside + def propagate_spatial_mask_over_time(data_array, mask): if data_array.ndim == 3 and mask.ndim == 2: nt = data_array.shape[0] for it in np.arange(nt): - mask_original = data_array[it,:].mask - data_array[it,:] = ma.array(data_array[it,:], mask=mask|mask_original) + mask_original = data_array[it, :].mask + data_array.mask[it,:] = mask | mask_original return data_array +def convert_lat_lon_2d_array(lon, lat): + if lon.ndim == 1 and lat.ndim == 1: + return np.meshgrid(lon, lat) + if lon.ndim == 2 and lat.ndim == 2: + return lon, lat
