| # |
| # Licensed to the Apache Software Foundation (ASF) under one or more |
| # contributor license agreements. See the NOTICE file distributed with |
| # this work for additional information regarding copyright ownership. |
| # The ASF licenses this file to You under the Apache License, Version 2.0 |
| # (the "License"); you may not use this file except in compliance with |
| # the License. You may obtain a copy of the License at |
| # |
| # http://www.apache.org/licenses/LICENSE-2.0 |
| # |
| # Unless required by applicable law or agreed to in writing, software |
| # distributed under the License is distributed on an "AS IS" BASIS, |
| # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| # See the License for the specific language governing permissions and |
| # limitations under the License. |
| # |
| |
| import datetime |
| import logging |
| |
| import netCDF4 |
| import numpy as np |
| import numpy.ma as ma |
| import scipy.ndimage |
| from matplotlib.path import Path |
| from scipy.interpolate import griddata |
| from scipy.ndimage import map_coordinates |
| |
| import ocw.utils as utils |
| from ocw import dataset as ds |
| |
| logger = logging.getLogger(__name__) |
| |
| |
| def temporal_subset(target_dataset, month_start, month_end, |
| average_each_year=False): |
| """ Temporally subset data given month_index. |
| |
| :param month_start: An integer for beginning month (Jan=1) |
| :type month_start: :class:`int` |
| |
| :param month_end: An integer for ending month (Jan=1) |
| :type month_end: :class:`int` |
| |
| :param target_dataset: Dataset object that needs temporal subsetting |
| :type target_dataset: Open Climate Workbench Dataset Object |
| |
| :param average_each_year: If True, output dataset is averaged for each year |
| :type average_each_year: :class:'boolean' |
| |
| :returns: A temporal subset OCW Dataset |
| :rtype: Open Climate Workbench Dataset Object |
| """ |
| |
| if month_start > month_end: |
| month_index = list(range(month_start, 13)) |
| month_index.extend(list(range(1, month_end + 1))) |
| else: |
| month_index = range(month_start, month_end + 1) |
| |
| dates = target_dataset.times |
| months = np.array([d.month for d in dates]) |
| time_index = [] |
| for m_value in month_index: |
| time_index.extend(list(np.where(months == m_value)[0])) |
| time_index.sort() |
| |
| new_dataset = ds.Dataset(target_dataset.lats, |
| target_dataset.lons, |
| target_dataset.times[time_index], |
| target_dataset.values[time_index, :], |
| variable=target_dataset.variable, |
| units=target_dataset.units, |
| name=target_dataset.name) |
| |
| if average_each_year: |
| new_times = new_dataset.times |
| nmonth = len(month_index) |
| ntime = new_times.size |
| nyear = ntime // nmonth |
| if ntime % nmonth != 0: |
| logger.warning("Number of times in dataset ({}) does not " |
| "divide evenly into {} year(s). Trimming data..." |
| .format(ntime, nyear)) |
| slc = utils.trim_dataset(new_dataset) |
| new_dataset.values = new_dataset.values[slc] |
| new_times = new_times[slc] |
| nyear = new_times.size // nmonth |
| |
| averaged_time = [] |
| ny, nx = target_dataset.values.shape[1:] |
| averaged_values = ma.zeros([nyear, ny, nx]) |
| for iyear in np.arange(nyear): |
| # centered time index of the season between month_start and |
| # month_end in each year |
| center_index = int(nmonth / 2 + iyear * nmonth) |
| if nmonth == 1: |
| center_index = iyear |
| averaged_time.append(new_times[center_index]) |
| averaged_values[iyear, :] = ma.average(new_dataset.values[ |
| nmonth * iyear: nmonth * iyear + nmonth, :], axis=0) |
| new_dataset = ds.Dataset(target_dataset.lats, |
| target_dataset.lons, |
| np.array(averaged_time), |
| averaged_values, |
| variable=target_dataset.variable, |
| units=target_dataset.units, |
| name=target_dataset.name) |
| return new_dataset |
| |
| |
| def temporal_rebin(target_dataset, temporal_resolution): |
| """ Rebin a Dataset to a new temporal resolution |
| |
| :param target_dataset: Dataset object that needs temporal rebinned |
| :type target_dataset: :class:`dataset.Dataset` |
| |
| :param temporal_resolution: The new temporal resolution |
| :type temporal_resolution: :mod:`string` |
| |
| :returns: A new temporally rebinned Dataset |
| :rtype: :class:`dataset.Dataset` |
| """ |
| # Decode the temporal resolution into a string format that |
| # _rcmes_calc_average_on_new_time_unit_K() can understand |
| |
| binned_values, binned_dates = _rcmes_calc_average_on_new_time_unit( |
| target_dataset.values, target_dataset.times, temporal_resolution) |
| binned_dates = np.array(binned_dates) |
| new_dataset = ds.Dataset(target_dataset.lats, |
| target_dataset.lons, |
| binned_dates, |
| binned_values, |
| variable=target_dataset.variable, |
| units=target_dataset.units, |
| name=target_dataset.name, |
| origin=target_dataset.origin) |
| return new_dataset |
| |
| |
| def temporal_rebin_with_time_index(target_dataset, nt_average): |
| """ Rebin a Dataset to a new temporal resolution |
| |
| :param target_dataset: Dataset object that needs temporal rebinned |
| :type target_dataset: :class:`dataset.Dataset` |
| |
| :param nt_average: Time resolution for the output datasets. |
| It is the same as the number of time indicies to be averaged. |
| length of time dimension in the rebinned dataset) = |
| (original time dimension length/nt_average) |
| :type nt_average: integer |
| |
| :returns: A new temporally rebinned Dataset |
| :rtype: :class:`dataset.Dataset` |
| """ |
| nt = target_dataset.times.size |
| if nt % nt_average != 0: |
| logger.warning('Length of time dimension must ' |
| 'be a multiple of nt_average') |
| # nt2 is the length of time dimension in the rebinned dataset |
| nt2 = nt // nt_average |
| binned_dates = target_dataset.times[np.arange(nt2) * nt_average] |
| binned_values = ma.zeros( |
| np.insert(target_dataset.values.shape[1:], 0, nt2)) |
| for it in np.arange(nt2): |
| binned_values[it, :] = ma.average( |
| target_dataset.values[nt_average * it: |
| nt_average * it + nt_average, |
| :], |
| axis=0) |
| new_dataset = ds.Dataset(target_dataset.lats, |
| target_dataset.lons, |
| binned_dates, |
| binned_values, |
| variable=target_dataset.variable, |
| units=target_dataset.units, |
| name=target_dataset.name, |
| origin=target_dataset.origin) |
| return new_dataset |
| |
| |
| def spatial_regrid(target_dataset, new_latitudes, new_longitudes, |
| boundary_check=True): |
| """ Regrid a Dataset using the new latitudes and longitudes |
| |
| :param target_dataset: Dataset object that needs spatially regridded |
| :type target_dataset: :class:`dataset.Dataset` |
| |
| :param new_latitudes: Array of latitudes |
| :type new_latitudes: :class:`numpy.ndarray` |
| |
| :param new_longitudes: Array of longitudes |
| :type new_longitudes: :class:`numpy.ndarray` |
| |
| :param boundary_check: Check if the regriding domain's boundaries |
| are outside target_dataset's domain |
| :type boundary_check: :class:'bool' |
| |
| :returns: A new spatially regridded Dataset |
| :rtype: :class:`dataset.Dataset` |
| """ |
| |
| # Create grids of the given lats and lons for the underlying API |
| # NOTE: np.meshgrid() requires inputs (x, y) and returns data |
| # 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: |
| new_lons, new_lats = np.meshgrid(new_longitudes, new_latitudes) |
| else: |
| new_lons = new_longitudes |
| new_lats = new_latitudes |
| |
| ny_old, nx_old = lats.shape |
| ny_new, nx_new = new_lats.shape |
| |
| for iy in np.arange(ny_old): |
| if not all(x < y for x, y in zip(lons[iy, :], lons[iy, 1:])): |
| lons[iy, :][lons[iy, :] < 0] = lons[iy, :][lons[iy, :] < 0] + 360. |
| |
| # Make masked array of shape (times, new_latitudes,new_longitudes) |
| new_values = ma.zeros([len(target_dataset.times), |
| ny_new, nx_new])+1.e+20 |
| |
| # 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: |
| # from south to north along the west boundary |
| for iy in np.arange(ny_old): |
| vertices.append([lons[iy, 0], lats[iy, 0]]) |
| # from west to east along the north boundary |
| for ix in np.arange(nx_old): |
| vertices.append([lons[-1, ix], lats[-1, ix]]) |
| # from north to south along the east boundary |
| for iy in np.arange(ny_old)[::-1]: |
| vertices.append([lons[iy, -1], lats[iy, -1]]) |
| # from east to west along the south boundary |
| for ix in np.arange(nx_old)[::-1]: |
| vertices.append([lons[0, ix], lats[0, ix]]) |
| path = Path(vertices) |
| new_xy_mask = np.ones(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]]) or not boundary_check: |
| new_xy_mask[iy, ix] = 0. |
| |
| new_index = np.where(new_xy_mask == 0.) |
| # Regrid the data on each time slice |
| for i in range(len(target_dataset.times)): |
| if len(target_dataset.times) == 1 and target_dataset.values.ndim == 2: |
| values_original = ma.array(target_dataset.values) |
| else: |
| values_original = ma.array(target_dataset.values[i]) |
| new_mask = np.copy(values_original.mask) |
| for shift in (-1, 1): |
| for axis in (0, 1): |
| q_shifted = np.roll(values_original, shift=shift, axis=axis) |
| if (np.where((values_original.mask == True) & (q_shifted.mask == False)))[0].size !=0: |
| index1 =np.where((values_original.mask == True) & (q_shifted.mask == False)) |
| n_indices = len(index1[0]) |
| values_original.data[index1] = q_shifted[index1] |
| new_mask[index1] = np.repeat(False, n_indices) |
| mask_index = np.where(~new_mask) |
| if new_mask.size != 1: |
| mask_index = np.where(~new_mask) |
| else: |
| mask_index = np.where(~np.isnan(values_original)) |
| new_values_temp = griddata((lons[mask_index], lats[mask_index]), values_original[mask_index], |
| (new_lons[new_index], |
| new_lats[new_index]), |
| method='linear') |
| # Make a masking map using nearest neighbour interpolation -use this to |
| # determine locations with MDI and mask these |
| qmdi = np.zeros_like(values_original) |
| values_true_indices = np.where(values_original.mask == True) |
| values_false_indices = np.where(values_original.mask == False) |
| qmdi[values_true_indices] = 1. |
| qmdi[values_false_indices] = 0. |
| qmdi_r = griddata((lons.flatten(), lats.flatten()), qmdi.flatten(), |
| (new_lons[new_index], |
| new_lats[new_index]), |
| method='nearest') |
| new_values_temp = ma.masked_where(qmdi_r != 0.0, new_values_temp) |
| # Combine missing data mask, with outside domain mask define above. |
| new_values[i, new_index[0], new_index[1]] = new_values_temp[:] |
| new_values[i,:] = ma.masked_equal(new_values[i,:], 1.e+20) |
| |
| # TODO: |
| # This will call down to the _congrid() function and the lat and lon |
| # axis will be adjusted with the time axis being held constant |
| |
| # Create a new Dataset Object to return using new data |
| regridded_dataset = ds.Dataset(new_latitudes, |
| new_longitudes, |
| target_dataset.times, |
| new_values, |
| variable=target_dataset.variable, |
| units=target_dataset.units, |
| name=target_dataset.name, |
| origin=target_dataset.origin) |
| return regridded_dataset |
| |
| |
| def ensemble(datasets): |
| """ |
| Generate a single dataset which is the mean of the input datasets |
| |
| An ensemble datasets combines input datasets assuming the all have |
| similar shape, dimensions, and units. |
| |
| :param datasets: Datasets to be used to compose the ensemble dataset from. |
| All Datasets must be the same shape. |
| :type datasets: :class:`list` of :class:`dataset.Dataset` |
| |
| :returns: New Dataset with a name of 'Dataset Ensemble' |
| :rtype: :class:`dataset.Dataset` |
| """ |
| _check_dataset_shapes(datasets) |
| dataset_values = [dataset.values for dataset in datasets] |
| ensemble_values = ma.mean(dataset_values, axis=0) |
| |
| # Build new dataset object from the input datasets and |
| # the ensemble values and return it |
| ensemble_dataset = ds.Dataset(datasets[0].lats, |
| datasets[0].lons, |
| datasets[0].times, |
| ensemble_values, |
| units=datasets[0].units, |
| name="Dataset Ensemble") |
| |
| return ensemble_dataset |
| |
| |
| 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. |
| :type subregion: :class:`dataset.Bounds` |
| |
| :param target_dataset: The Dataset object to subset. |
| :type target_dataset: :class:`dataset.Dataset` |
| |
| :param subregion_name: The subset-ed Dataset name |
| :type subregion_name: :mod:`string` |
| |
| :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` |
| |
| :raises: ValueError |
| ''' |
| |
| if not subregion.start: |
| time_start = target_dataset.times[0] |
| time_end = target_dataset.times[-1] |
| else: |
| time_start = subregion.start |
| time_end = subregion.end |
| |
| if not subregion_name: |
| subregion_name = target_dataset.name |
| |
| if subregion.lat_min is not None: |
| _are_bounds_contained_by_dataset(target_dataset, subregion) |
| |
| if target_dataset.lats.ndim == 2 and target_dataset.lons.ndim == 2: |
| temporal_subset = temporal_slice( |
| target_dataset, time_start, time_end) |
| nt = temporal_subset.values.shape[0] |
| y_index, x_index = np.where( |
| (target_dataset.lats >= subregion.lat_max) | ( |
| target_dataset.lats <= subregion.lat_min) | |
| (target_dataset.lons >= subregion.lon_max) | ( |
| target_dataset.lons <= subregion.lon_min)) |
| for it in np.arange(nt): |
| temporal_subset.values[it, y_index, x_index] = 1.e+20 |
| new_values = ma.masked_equal( |
| temporal_subset.values, 1.e+20) |
| return ds.Dataset( |
| target_dataset.lats, |
| target_dataset.lons, |
| temporal_subset.times, |
| new_values, |
| variable=target_dataset.variable, |
| units=target_dataset.units, |
| name=subregion_name, |
| origin=target_dataset.origin) |
| |
| elif target_dataset.lats.ndim == 1 and target_dataset.lons.ndim == 1: |
| # Get subregion indices into subregion data |
| dataset_slices = _get_subregion_slice_indices(target_dataset, |
| subregion) |
| # Slice the values array with our calculated slice indices |
| if target_dataset.values.ndim == 2: |
| subset_values = ma.zeros([len(target_dataset.values[ |
| dataset_slices["lat_start"]:dataset_slices["lat_end"]]), |
| len(target_dataset.values[ |
| dataset_slices["lon_start"]:dataset_slices["lon_end"]])]) |
| |
| subset_values = target_dataset.values[ |
| dataset_slices["lat_start"]:dataset_slices["lat_end"] + 1, |
| dataset_slices["lon_start"]:dataset_slices["lon_end"] + 1] |
| |
| elif target_dataset.values.ndim == 3: |
| subset_values = ma.zeros([len(target_dataset.values[ |
| dataset_slices["time_start"]:dataset_slices["time_end"]]), |
| len(target_dataset.values[ |
| dataset_slices["lat_start"]:dataset_slices["lat_end"]]), |
| len(target_dataset.values[ |
| dataset_slices["lon_start"]:dataset_slices["lon_end"]])]) |
| |
| subset_values = target_dataset.values[ |
| dataset_slices["time_start"]:dataset_slices["time_end"] + 1, |
| dataset_slices["lat_start"]:dataset_slices["lat_end"] + 1, |
| dataset_slices["lon_start"]:dataset_slices["lon_end"] + 1] |
| # Build new dataset with subset information |
| return ds.Dataset( |
| # Slice the lats array with our calculated slice indices |
| target_dataset.lats[dataset_slices["lat_start"]: |
| dataset_slices["lat_end"] + 1], |
| # Slice the lons array with our calculated slice indices |
| target_dataset.lons[dataset_slices["lon_start"]: |
| dataset_slices["lon_end"] + 1], |
| # Slice the times array with our calculated slice indices |
| target_dataset.times[dataset_slices["time_start"]: |
| dataset_slices["time_end"] + 1], |
| # Slice the values array with our calculated slice indices |
| subset_values, |
| variable=target_dataset.variable, |
| units=target_dataset.units, |
| name=subregion_name, |
| origin=target_dataset.origin |
| ) |
| |
| if subregion.boundary_type == 'us_states' or subregion.boundary_type == 'countries': |
| temporal_subset = temporal_slice( |
| target_dataset, time_start, time_end) |
| spatial_mask = utils.mask_using_shapefile_info(target_dataset.lons, target_dataset.lats, |
| subregion.masked_regions, extract=extract) |
| subset_values = utils.propagate_spatial_mask_over_time( |
| temporal_subset.values, mask=spatial_mask) |
| return ds.Dataset( |
| target_dataset.lats, |
| target_dataset.lons, |
| temporal_subset.times, |
| subset_values, |
| variable=target_dataset.variable, |
| units=target_dataset.units, |
| name=subregion_name, |
| origin=target_dataset.origin) |
| |
| if subregion.boundary_type == 'user': |
| temporal_subset = temporal_slice( |
| target_dataset, time_start, time_end) |
| 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) |
| subset_values = utils.propagate_spatial_mask_over_time( |
| temporal_subset.values, mask=spatial_mask) |
| return ds.Dataset( |
| target_dataset.lats, |
| target_dataset.lons, |
| temporal_subset.times, |
| subset_values, |
| variable=target_dataset.variable, |
| units=target_dataset.units, |
| name=subregion_name, |
| origin=target_dataset.origin) |
| |
| def temporal_slice(target_dataset, start_time, end_time): |
| '''Temporally slice given dataset(s) with subregion information. This does not |
| spatially subset the target_Dataset |
| |
| :param start_time: start time |
| :type start_time: :class:'int' |
| |
| :param end_time: end time |
| :type end_time: :class:'datetime.datetime' |
| |
| :param target_dataset: The Dataset object to subset. |
| :type target_dataset: :class:`dataset.Dataset` |
| |
| :returns: The subset-ed Dataset object |
| :rtype: :class:`dataset.Dataset` |
| |
| :raises: ValueError |
| ''' |
| |
| # https://issues.apache.org/jira/browse/CLIMATE-938 |
| # netCDF datetimes allow for a variety of calendars while Python has |
| # only one. This would throw an error about a calendar mismatch when |
| # comparing a Python datetime object to a netcdf datetime object. |
| # Cast the date as best we can so the comparison will compare like |
| # data types This will still throw an excdeption if the start / end date are |
| # not valid in given calendar. February 29th in a DatetimeNoLeap calendar for example. |
| slice_start_time = start_time |
| slice_end_time = end_time |
| |
| if type(target_dataset.times[0]) != datetime.datetime: |
| slice_start_time =\ |
| type(target_dataset.times.item(0))(start_time.year, start_time.month, start_time.day, |
| start_time.hour, start_time.minute, start_time.second) |
| |
| slice_end_time =\ |
| type(target_dataset.times.item(0))(end_time.year, end_time.month, end_time.day, |
| end_time.hour, end_time.minute, end_time.second) |
| |
| start_time_index = np.where( |
| target_dataset.times >= slice_start_time)[0][0] |
| |
| end_time_index = np.where( |
| target_dataset.times <= slice_end_time)[0][-1] |
| |
| new_times = target_dataset.times[start_time_index:end_time_index + 1] |
| new_values = target_dataset.values[start_time_index:end_time_index + 1, :] |
| |
| return ds.Dataset( |
| target_dataset.lats, |
| target_dataset.lons, |
| new_times, |
| new_values, |
| variable=target_dataset.variable, |
| units=target_dataset.units, |
| origin=target_dataset.origin) |
| |
| def safe_subset(target_dataset, subregion, subregion_name=None): |
| '''Safely subset given dataset with subregion information |
| A standard subset requires that the provided subregion be entirely |
| contained within the datasets bounds. `safe_subset` returns the |
| overlap of the subregion and dataset without returning an error. |
| |
| :param subregion: The Bounds with which to subset the target Dataset. |
| :type subregion: :class:`dataset.Bounds` |
| |
| :param target_dataset: The Dataset object to subset. |
| :type target_dataset: :class:`dataset.Dataset` |
| |
| :param subregion_name: The subset-ed Dataset name |
| :type subregion_name: :mod:`string` |
| |
| :returns: The subset-ed Dataset object |
| :rtype: :class:`dataset.Dataset` |
| ''' |
| |
| lat_min, lat_max, lon_min, lon_max = target_dataset.spatial_boundaries() |
| start, end = target_dataset.temporal_boundaries() |
| if subregion.lat_min: |
| if subregion.lat_min < lat_min: |
| subregion.lat_min = lat_min |
| |
| if subregion.lat_max > lat_max: |
| subregion.lat_max = lat_max |
| |
| if subregion.lon_min < lon_min: |
| subregion.lon_min = lon_min |
| |
| if subregion.lon_max > lon_max: |
| subregion.lon_max = lon_max |
| |
| if subregion.start: |
| if subregion.start < start: |
| subregion.start = start |
| |
| if subregion.end: |
| if subregion.end > end: |
| subregion.end = end |
| |
| return subset(target_dataset, subregion, subregion_name) |
| |
| |
| def normalize_dataset_datetimes(dataset, timestep): |
| ''' Normalize Dataset datetime values. |
| |
| Force daily to an hour time value of 00:00:00. |
| Force monthly data to the first of the month at midnight. |
| |
| :param dataset: The Dataset which will have its time value normalized. |
| :type dataset: :class:`dataset.Dataset` |
| |
| :param timestep: The timestep of the Dataset's values. Either 'daily' or |
| 'monthly'. |
| :type timestep: :mod:`string` |
| |
| :returns: A new Dataset with normalized datetime values. |
| :rtype: :class:`dataset.Dataset` |
| ''' |
| new_times = _rcmes_normalize_datetimes(dataset.times, timestep) |
| return ds.Dataset( |
| dataset.lats, |
| dataset.lons, |
| np.array(new_times), |
| dataset.values, |
| variable=dataset.variable, |
| units=dataset.units, |
| name=dataset.name, |
| origin=dataset.origin |
| ) |
| |
| |
| def write_netcdf(dataset, path, compress=True): |
| ''' Write a dataset to a NetCDF file. |
| |
| :param dataset: The dataset to write. |
| :type dataset: :class:`dataset.Dataset` |
| |
| :param path: The output file path. |
| :type path: :mod:`string` |
| ''' |
| out_file = netCDF4.Dataset(path, 'w', format='NETCDF4') |
| |
| # Set attribute lengths |
| if dataset.lats.ndim == 2: |
| lat_len = dataset.lats.shape[0] |
| lon_len = dataset.lons.shape[1] |
| lat_dim_info = ('y', 'x') |
| lon_dim_info = ('y', 'x') |
| |
| else: |
| lat_len = len(dataset.lats) |
| lon_len = len(dataset.lons) |
| lat_dim_info = ('y',) |
| lon_dim_info = ('x',) |
| |
| time_len = len(dataset.times) |
| |
| # Create attribute dimensions |
| out_file.createDimension('y', lat_len) |
| out_file.createDimension('x', lon_len) |
| out_file.createDimension('time', None) |
| |
| # Create variables |
| lats = out_file.createVariable('lat', 'f8', lat_dim_info, zlib=compress) |
| lons = out_file.createVariable('lon', 'f8', lon_dim_info, zlib=compress) |
| times = out_file.createVariable('time', 'f8', ('time'), zlib=compress) |
| |
| var_name = dataset.variable if dataset.variable else 'var' |
| values = out_file.createVariable(var_name, |
| 'f8', |
| ('time', 'y', 'x'), |
| zlib=compress) |
| |
| # Set the time variable units |
| # We don't deal with hourly/minutely/anything-less-than-a-day data so |
| # we can safely stick with a 'days since' offset here. Note that the |
| # NetCDF4 helper date2num doesn't support 'months' or 'years' instead |
| # of days. |
| times.units = "days since %s" % dataset.times[0] |
| |
| # Store the dataset's values |
| lats[:] = dataset.lats |
| lons[:] = dataset.lons |
| times[:] = netCDF4.date2num(dataset.times, times.units) |
| values[:] = dataset.values |
| values.units = dataset.units |
| |
| out_file.close() |
| |
| |
| def write_netcdf_multiple_datasets_with_subregions(ref_dataset, ref_name, |
| model_dataset_array, |
| model_names, |
| path, |
| subregions=None, |
| subregion_array=None, |
| ref_subregion_mean=None, |
| ref_subregion_std=None, |
| model_subregion_mean=None, |
| model_subregion_std=None): |
| # Write multiple reference and model datasets and their subregional means |
| # and standard deivations in a NetCDF file. |
| |
| # :To be updated |
| # |
| out_file = netCDF4.Dataset(path, 'w', format='NETCDF4') |
| |
| dataset = ref_dataset |
| # Set attribute lenghts |
| nobs = 1 |
| nmodel = len(model_dataset_array) |
| lat_len, lon_len = dataset.values.shape[1:] |
| lat_ndim = dataset.lats.ndim |
| lon_ndim = dataset.lons.ndim |
| time_len = len(dataset.times) |
| |
| if subregions is not None: |
| nsubregion = len(subregions) |
| |
| # Create attribute dimensions |
| out_file.createDimension('y', lat_len) |
| out_file.createDimension('x', lon_len) |
| out_file.createDimension('time', time_len) |
| |
| # Create variables and store the values |
| if lat_ndim == 2: |
| lats = out_file.createVariable('lat', 'f8', ('y', 'x')) |
| else: |
| lats = out_file.createVariable('lat', 'f8', ('y')) |
| lats[:] = dataset.lats |
| if lon_ndim == 2: |
| lons = out_file.createVariable('lon', 'f8', ('y', 'x')) |
| else: |
| lons = out_file.createVariable('lon', 'f8', ('x')) |
| lons[:] = dataset.lons |
| times = out_file.createVariable('time', 'f8', ('time',)) |
| times.units = "days since %s" % dataset.times[0] |
| times[:] = netCDF4.date2num(dataset.times, times.units) |
| |
| # mask_array = np.zeros([time_len, lat_len, lon_len]) |
| # for iobs in np.arange(nobs): |
| # index = np.where(ref_dataset_array[iobs].values.mask[:] == True) |
| # mask_array[index] = 1 |
| out_file.createVariable(ref_name, 'f8', ('time', 'y', 'x')) |
| out_file.variables[ref_name][:] = ref_dataset.values |
| out_file.variables[ref_name].units = ref_dataset.units |
| for imodel in np.arange(nmodel): |
| out_file.createVariable(model_names[imodel], 'f8', ('time', 'y', 'x')) |
| # out_file.variables[model_names[imodel]][:] = ma.array( |
| # model_dataset_array[imodel].values, mask = mask_array) |
| out_file.variables[model_names[imodel]][:] = model_dataset_array[ |
| imodel].values |
| out_file.variables[model_names[imodel]].units = model_dataset_array[ |
| imodel].units |
| |
| if subregions is not None: |
| out_file.createVariable('subregion_array', 'i4', ('y', 'x')) |
| out_file.variables['subregion_array'][:] = subregion_array[:] |
| nsubregion = len(subregions) |
| out_file.createDimension('nsubregion', nsubregion) |
| out_file.createDimension('nobs', nobs) |
| out_file.createDimension('nmodel', nmodel) |
| out_file.createVariable('obs_subregion_mean', 'f8', ('nobs', |
| 'time', |
| 'nsubregion')) |
| out_file.variables['obs_subregion_mean'][:] = ref_subregion_mean[:] |
| out_file.createVariable('obs_subregion_std', 'f8', ('nobs', |
| 'time', |
| 'nsubregion')) |
| out_file.variables['obs_subregion_std'][:] = ref_subregion_std[:] |
| out_file.createVariable('model_subregion_mean', 'f8', ('nmodel', |
| 'time', |
| 'nsubregion')) |
| out_file.variables['model_subregion_mean'][:] = model_subregion_mean[:] |
| out_file.createVariable('model_subregion_std', 'f8', ('nmodel', |
| 'time', |
| 'nsubregion')) |
| out_file.variables['model_subregion_std'][:] = model_subregion_std[:] |
| |
| out_file.close() |
| |
| |
| def water_flux_unit_conversion(dataset): |
| ''' Convert water flux variables units as necessary |
| |
| Convert full SI units water flux units to more common units. |
| |
| :param dataset: The dataset to convert. |
| :type dataset: :class:`dataset.Dataset` |
| |
| :returns: A Dataset with values converted to new units. |
| :rtype: :class:`dataset.Dataset` |
| ''' |
| water_flux_variables = ['pr', 'prec', 'evspsbl', 'mrro', 'swe'] |
| variable = dataset.variable.lower() |
| |
| if any(sub_string in variable for sub_string in water_flux_variables): |
| dataset_units = dataset.units.lower() |
| if variable in 'swe': |
| if any(unit in dataset_units for unit in ['m', 'meter']): |
| dataset.values = 1.e3 * dataset.values |
| dataset.units = 'km' |
| else: |
| if any(unit in dataset_units |
| for unit in ['kg m-2 s-1', 'mm s-1', 'mm/sec']): |
| dataset.values = 86400. * dataset.values |
| dataset.units = 'mm/day' |
| |
| return dataset |
| |
| |
| def temperature_unit_conversion(dataset): |
| ''' Convert temperature units as necessary \ |
| Automatically convert Celcius to Kelvin in the given dataset. |
| |
| :param dataset: The dataset for which units should be updated. \ |
| :type dataset; :class:`dataset.Dataset` |
| |
| :returns: The dataset with (potentially) updated units. \ |
| :rtype: :class:`dataset.Dataset` |
| |
| ''' |
| |
| temperature_variables = ['temp', 'tas', 'tasmax', 'taxmin', 'T', 'tg'] |
| variable = dataset.variable.lower() |
| |
| if any(sub_string in variable for sub_string in temperature_variables): |
| dataset_units = dataset.units.lower() |
| if dataset_units == 'c' or dataset_units == 'celsius': |
| dataset.values = 273.15 + dataset.values |
| dataset.units = 'K' |
| |
| return dataset |
| |
| |
| def variable_unit_conversion(dataset): |
| ''' Convert water flux or temperature variables units as necessary |
| |
| For water flux variables, convert full SI units water flux units |
| to more common units. |
| For temperature, convert Celcius to Kelvin. |
| |
| :param dataset: The dataset to convert. |
| :type dataset: :class:`dataset.Dataset` |
| |
| :returns: A Dataset with values converted to new units. |
| :rtype: :class:`dataset.Dataset` |
| ''' |
| |
| dataset = water_flux_unit_conversion(dataset) |
| dataset = temperature_unit_conversion(dataset) |
| |
| return dataset |
| |
| |
| def _rcmes_normalize_datetimes(datetimes, timestep): |
| """ Normalize Dataset datetime values. |
| |
| Force daily to an hour time value of 00:00:00. |
| Force monthly data to the first of the month at midnight. |
| |
| :param datetimes: The datetimes to normalize. |
| :type datetimes: List of `datetime` values. |
| |
| :param timestep: The flag for how to normalize the datetimes. |
| :type timestep: String |
| """ |
| |
| normalDatetimes = [] |
| if timestep.lower() == 'monthly': |
| for inputDatetime in datetimes: |
| if inputDatetime.day != 1: |
| # Clean the inputDatetime |
| inputDatetimeString = inputDatetime.strftime('%Y%m%d') |
| normalInputDatetimeString = inputDatetimeString[:6] + '01' |
| inputDatetime = datetime.datetime.strptime( |
| normalInputDatetimeString, '%Y%m%d') |
| |
| normalDatetimes.append(inputDatetime) |
| |
| elif timestep.lower() == 'daily': |
| for inputDatetime in datetimes: |
| if (inputDatetime.hour != 0 or inputDatetime.minute != 0 or |
| inputDatetime.second != 0): |
| datetimeString = inputDatetime.strftime('%Y%m%d%H%M%S') |
| normalDatetimeString = datetimeString[:8] + '000000' |
| inputDatetime = datetime.datetime.strptime( |
| normalDatetimeString, '%Y%m%d%H%M%S') |
| |
| normalDatetimes.append(inputDatetime) |
| |
| return normalDatetimes |
| |
| |
| def mask_missing_data(dataset_array): |
| ''' Check missing values in observation and model datasets. |
| If any of dataset in dataset_array has missing values at a grid point, |
| the values at the grid point in all other datasets are masked. |
| :param dataset_array: an array of OCW datasets |
| ''' |
| |
| mask_array = np.zeros(dataset_array[0].values.shape) |
| for dataset in dataset_array: |
| # CLIMATE-797 - Not every array passed in will be a masked array. |
| # For those that are, action based on the mask passed in. |
| # For those that are not, take no action (else AttributeError). |
| if hasattr(dataset.values, 'mask'): |
| index = np.where(dataset.values.mask == True) |
| if index[0].size > 0: |
| mask_array[index] = 1 |
| masked_array = [] |
| for dataset in dataset_array: |
| dataset.values = ma.array(dataset.values, mask=mask_array) |
| masked_array.append(dataset) |
| return [masked_dataset for masked_dataset in masked_array] |
| |
| |
| def deseasonalize_dataset(dataset): |
| '''Calculate daily climatology and subtract the climatology from |
| the input dataset |
| |
| :param dataset: The dataset to convert. |
| :type dataset: :class:`dataset.Dataset` |
| |
| :returns: A Dataset with values deseasonalized. |
| :rtype: :class:`dataset.Dataset` |
| ''' |
| |
| days = [d.month * 100. + d.day for d in dataset.times] |
| days_sorted = np.unique(days) |
| ndays = days_sorted.size |
| nt, ny, nx = dataset.values.shape |
| values_clim = ma.zeros([ndays, ny, nx]) |
| for iday, day in enumerate(days_sorted): |
| t_index = np.where(days == day)[0] |
| values_clim[iday, :] = ma.mean(dataset.values[t_index, :], axis=0) |
| for iday, day in enumerate(days_sorted): |
| t_index = np.where(days == day)[0] |
| dataset.values[t_index, :] = dataset.values[ |
| t_index, :] - values_clim[iday, :] |
| return dataset |
| |
| |
| def _rcmes_spatial_regrid(spatial_values, lat, lon, lat2, lon2, order=1): |
| ''' |
| Spatial regrid from one set of lat,lon values onto a new set (lat2,lon2) |
| |
| :param spatial_values: Values in a spatial grid that need to be regridded |
| :type spatial_values: 2d masked numpy array. shape (latitude, longitude) |
| :param lat: Grid of latitude values which map to the spatial values |
| :type lat: 2d numpy array. shape(latitudes, longitudes) |
| :param lon: Grid of longitude values which map to the spatial values |
| :type lon: 2d numpy array. shape(latitudes, longitudes) |
| :param lat2: Grid of NEW latitude values to regrid the spatial_values onto |
| :type lat2: 2d numpy array. shape(latitudes, longitudes) |
| :param lon2: Grid of NEW longitude values to regrid the spatial_values onto |
| :type lon2: 2d numpy array. shape(latitudes, longitudes) |
| :param order: Interpolation order flag. 1=bi-linear, 3=cubic spline |
| :type order: [optional] Integer |
| |
| :returns: 2d masked numpy array with shape(len(lat2), len(lon2)) |
| :rtype: (float, float) |
| ''' |
| |
| nlat = spatial_values.shape[0] |
| nlon = spatial_values.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_max_indices = np.where(loni > lon.max()) |
| loni_min_indices = np.where(loni < lon.min()) |
| loni[loni_max_indices] = lon.max() |
| loni[loni_min_indices] = lon.min() |
| |
| # To deal with the "hard" break, we'll have to treat y differently, |
| # so we're just setting the min here... |
| lati_max_indices = np.where(lati > lat.max()) |
| lati_min_indices = np.where(lati < lat.min()) |
| lati[lati_max_indices] = lat.max() |
| lati[lati_min_indices] = 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()) |
| |
| """ |
| TODO: Review this docstring and see if it still holds true. |
| NOTE: This function doesn't use MDI currently. These are legacy comments |
| |
| 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(spatial_values, shift=shift, axis=axis) |
| idx = ~q_shifted.mask * spatial_values.mask |
| indices = np.where(idx) |
| spatial_values.data[indices] = q_shifted[indices] |
| |
| # Now we actually interpolate |
| # map_coordinates does cubic interpolation by default, |
| # use "order=1" to preform bilinear interpolation instead... |
| |
| regridded_values = map_coordinates(spatial_values, |
| [lati, loni], |
| order=order) |
| regridded_values = regridded_values.reshape([nlat2, nlon2]) |
| |
| # Set values to missing data outside of original domain |
| regridded_values = ma.masked_array(regridded_values, |
| 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(spatial_values) |
| spatial_indices_true = np.where(spatial_values.mask) |
| spatial_indices_false = np.where(not spatial_values.mask) |
| qmdi[spatial_indices_true] = 1. |
| qmdi[spatial_indices_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. |
| regridded_values.mask = np.logical_or(mdimask, regridded_values.mask) |
| |
| return regridded_values |
| |
| |
| def _rcmes_create_mask_using_threshold(masked_array, threshold=0.5): |
| '''Mask an array if percent of values missing data is above a threshold. |
| |
| For each value along an axis, if the proportion of steps that are missing |
| data is above ``threshold`` then the value is marked as missing data. |
| |
| ..note:: The 0th axis is currently always used. |
| |
| :param masked_array: Masked array of data |
| :type masked_array: Numpy Masked Array |
| :param threshold: (optional) Threshold proportion above which a value is |
| marked as missing data. |
| :type threshold: Float |
| |
| :returns: A Numpy array describing the mask for masked_array. |
| ''' |
| |
| # try, except used as some model files don't have a full mask, but a single |
| # bool the except catches this situation and deals with it appropriately. |
| try: |
| nT = masked_array.mask.shape[0] |
| |
| # For each pixel, count how many times are masked. |
| nMasked = masked_array.mask[:, :, :].sum(axis=0) |
| |
| # Define new mask as when a pixel has over a defined threshold ratio of |
| # masked data |
| # e.g. if the threshold is 75%, and there are 10 times, |
| # then a pixel will be masked if more than 5 times are masked. |
| mymask = nMasked > (nT * threshold) |
| |
| except: |
| mymask = np.zeros_like(masked_array.data[0, :, :]) |
| |
| return mymask |
| |
| |
| def _rcmes_calc_average_on_new_time_unit(data, dates, unit): |
| """ Rebin 3d array and list of dates using the provided unit parameter |
| |
| :param data: Input data that needs to be averaged |
| :type data: 3D masked numpy array of shape (times, lats, lons) |
| :param dates: List of dates that correspond to the given data values |
| :type dates: Python datetime objects |
| :param unit: Time unit to average the data into |
| :type unit: String matching one of these values : |
| full | annual | monthly | daily |
| |
| :returns: meanstorem, newTimesList |
| :rtype: 3D numpy masked array the same shape as the input array, |
| list of python datetime objects |
| :raises: ValueError |
| """ |
| |
| # Check if the user-selected temporal grid is valid. If not, EXIT |
| acceptable = ((unit == 'full') | (unit == 'annual') | |
| (unit == 'monthly') | (unit == 'daily')) |
| if not acceptable: |
| raise ValueError('Error: unknown unit type selected ' |
| 'for time averaging: EXIT') |
| |
| nt, ny, nx = data.shape |
| if unit == 'full': |
| new_data = ma.mean(data, axis=0) |
| new_date = [dates[dates.size // 2]] |
| if unit == 'annual': |
| years = [d.year for d in dates] |
| years_sorted = np.unique(years) |
| new_data = ma.zeros([years_sorted.size, ny, nx]) |
| it = 0 |
| new_date = [] |
| for year in years_sorted: |
| index = np.where(years == year)[0] |
| new_data[it, :] = ma.mean(data[index, :], axis=0) |
| new_date.append(datetime.datetime(year=year, month=7, day=2)) |
| it = it + 1 |
| if unit == 'monthly': |
| years = [d.year for d in dates] |
| years_sorted = np.unique(years) |
| months = [d.month for d in dates] |
| months_sorted = np.unique(months) |
| |
| new_data = ma.zeros([years_sorted.size * months_sorted.size, ny, nx]) |
| it = 0 |
| new_date = [] |
| for year in years_sorted: |
| for month in months_sorted: |
| index = np.where((years == year) & (months == month))[0] |
| new_data[it, :] = ma.mean(data[index, :], axis=0) |
| new_date.append(datetime.datetime(year=year, |
| month=month, |
| day=15)) |
| it = it + 1 |
| if unit == 'daily': |
| days = [d.year * 10000. + d.month * 100. + d.day for d in dates] |
| days_sorted = np.unique(days) |
| |
| new_data = ma.zeros([days_sorted.size, ny, nx]) |
| it = 0 |
| new_date = [] |
| for day in days_sorted: |
| index = np.where(days == day)[0] |
| new_data[it, :] = ma.mean(data[index, :], axis=0) |
| y = int(day // 10000) |
| m = int(day % 10000) // 100 |
| d = int(day % 100) |
| new_date.append(datetime.datetime(year=y, month=m, day=d)) |
| it = it + 1 |
| |
| return new_data, np.array(new_date) |
| |
| |
| def _rcmes_calc_average_on_new_time_unit_K(data, dates, unit): |
| """ Rebin 3d array and list of dates using the provided unit parameter |
| |
| :param data: Input data that needs to be averaged |
| :type data: 3D masked numpy array of shape (times, lats, lons) |
| :param dates: List of dates that correspond to the given data values |
| :type dates: Python datetime objects |
| :param unit: Time unit to average the data into |
| :type unit: String matching one of these values : |
| full | annual | monthly | daily |
| |
| :returns: meanstorem, newTimesList |
| :rtype: 3D numpy masked array the same shape as the input array, |
| list of python datetime objects |
| :raises: ValueError |
| """ |
| |
| # Check if the user-selected temporal grid is valid. If not, EXIT |
| acceptable = ((unit == 'full') | (unit == 'annual') | |
| (unit == 'monthly') | (unit == 'daily')) |
| if not acceptable: |
| raise ValueError('Error: unknown unit type selected ' |
| 'for time averaging: EXIT') |
| |
| # Calculate arrays of: annual timeseries: year (2007,2007), |
| # monthly time series: year-month (200701,200702), |
| # daily timeseries: year-month-day (20070101,20070102) |
| # depending on user-selected averaging period. |
| |
| # Year list |
| if unit == 'annual': |
| timeunits = np.array([int(d.strftime("%Y")) for d in dates]) |
| unique_times = np.unique(timeunits) |
| |
| # YearMonth format list |
| if unit == 'monthly': |
| timeunits = np.array([int(d.strftime("%Y%m")) for d in dates]) |
| unique_times = np.unique(timeunits) |
| |
| # YearMonthDay format list |
| if unit == 'daily': |
| timeunits = np.array([int(d.strftime("%Y%m%d")) for d in dates]) |
| unique_times = np.unique(timeunits) |
| |
| # TODO: add pentad setting using Julian days? |
| |
| # Full list: a special case |
| if unit == 'full': |
| # Calculating means data over the entire time range: |
| # i.e., annual-mean climatology |
| timeunits = [] |
| for i in np.arange(len(dates)): |
| # i.e. we just want the same value for all times. |
| timeunits.append(999) |
| timeunits = np.array(timeunits, dtype=int) |
| unique_times = np.unique(timeunits) |
| |
| # empty list to store new times |
| newTimesList = [] |
| |
| # Decide whether or not you need to do any time averaging. |
| # i.e. if data are already on required time unit then just |
| # pass data through and calculate and return representative |
| # datetimes. |
| processing_required = True |
| if len(timeunits) == (len(unique_times)): |
| processing_required = False |
| |
| # 1D data arrays, i.e. time series |
| if data.ndim == 1: |
| # Create array to store the resulting data |
| meanstore = np.zeros(len(unique_times)) |
| |
| # Calculate the means across each unique time unit |
| i = 0 |
| for myunit in unique_times: |
| if processing_required: |
| datam = ma.masked_array(data, timeunits != myunit) |
| meanstore[i] = ma.average(datam) |
| |
| # construct new times list |
| yyyy, mm, dd = _create_new_year_month_day(myunit, dates) |
| newTimesList.append(datetime.datetime(yyyy, mm, dd, 0, 0, 0, 0)) |
| i = i + 1 |
| |
| # 3D data arrays |
| if data.ndim == 3: |
| # Create array to store the resulting data |
| meanstore = np.zeros([len(unique_times), data.shape[1], data.shape[2]]) |
| |
| # Calculate the means across each unique time unit |
| i = 0 |
| datamask_store = [] |
| for myunit in unique_times: |
| if processing_required: |
| mask = np.zeros_like(data) |
| mask[timeunits != myunit, :, :] = 1.0 |
| # Calculate missing data mask within each time unit... |
| datamask_at_this_timeunit = np.zeros_like(data) |
| datamask_at_this_timeunit[:] = _rcmes_create_mask_using_threshold( |
| data[timeunits == myunit, :, :], threshold=0.75) |
| # Store results for masking later |
| datamask_store.append(datamask_at_this_timeunit[0]) |
| # Calculate means for each pixel in this time unit, ignoring |
| # missing data (using masked array). |
| datam = ma.masked_array(data, np.logical_or(mask, |
| datamask_at_this_timeunit)) |
| meanstore[i, :, :] = ma.average(datam, axis=0) |
| # construct new times list |
| yyyy, mm, dd = _create_new_year_month_day(myunit, dates) |
| newTimesList.append(datetime.datetime(yyyy, mm, dd)) |
| i += 1 |
| |
| if not processing_required: |
| meanstorem = data |
| |
| if processing_required: |
| # Create masked array (using missing data mask defined above) |
| datamask_store = np.array(datamask_store) |
| meanstorem = ma.masked_array(meanstore, datamask_store) |
| |
| return meanstorem, newTimesList |
| |
| |
| def _create_new_year_month_day(time_unit, dates): |
| smyunit = str(time_unit) |
| if len(smyunit) == 4: # YYYY |
| yyyy = int(smyunit[0:4]) |
| mm = 1 |
| dd = 1 |
| elif len(smyunit) == 6: # YYYYMM |
| yyyy = int(smyunit[0:4]) |
| mm = int(smyunit[4:6]) |
| dd = 1 |
| elif len(smyunit) == 8: # YYYYMMDD |
| yyyy = int(smyunit[0:4]) |
| mm = int(smyunit[4:6]) |
| dd = int(smyunit[6:8]) |
| elif len(smyunit) == 3: # Full time range |
| # Need to set an appropriate time representing the mid-point of the |
| # entire time span |
| dt = dates[-1] - dates[0] |
| halfway = dates[0] + (dt / 2) |
| yyyy = int(halfway.year) |
| mm = int(halfway.month) |
| dd = int(halfway.day) |
| |
| return (yyyy, mm, dd) |
| |
| |
| def _congrid(a, newdims, method='linear', centre=False, minusone=False): |
| ''' |
| This function is from http://wiki.scipy.org/Cookbook/Rebinning - Example 3 |
| It has been refactored and changed a bit, but the original functionality |
| has been preserved. |
| Arbitrary resampling of source array to new dimension sizes. |
| Currently only supports maintaining the same number of dimensions. |
| To use 1-D arrays, first promote them to shape (x,1). |
| |
| Uses the same parameters and creates the same co-ordinate lookup points |
| as IDL''s congrid routine, which apparently originally came from a VAX/VMS |
| routine of the same name. |
| |
| method: |
| neighbour - closest value from original data |
| nearest and linear - uses n x 1-D interpolations using |
| scipy.interpolate.interp1d |
| (see Numerical Recipes for validity of use of n 1-D interpolations) |
| spline - uses ndimage.map_coordinates |
| |
| centre: |
| True - interpolation points are at the centres of the bins |
| False - points are at the front edge of the bin |
| |
| minusone: |
| For example- inarray.shape = (i,j) & new dimensions = (x,y) |
| False - inarray is resampled by factors of (i/x) * (j/y) |
| True - inarray is resampled by(i-1)/(x-1) * (j-1)/(y-1) |
| This prevents extrapolation one element beyond bounds of input array. |
| ''' |
| if a.dtype not in [np.float64, np.float32]: |
| a = np.cast[float](a) |
| |
| # this will merely take the True/False input and convert it to an |
| # array(1) or array(0) |
| m1 = np.cast[int](minusone) |
| # this also casts the True False input into a floating point number |
| # of 0.5 or 0.0 |
| ofs = np.cast[int](centre) * 0.5 |
| |
| old = np.array(a.shape) |
| ndims = len(a.shape) |
| if len(newdims) != ndims: |
| print("[congrid] dimensions error. " |
| "This routine currently only supports " |
| "rebinning to the same number of dimensions.") |
| return None |
| newdims = np.asarray(newdims, dtype=float) |
| dimlist = [] |
| |
| if method == 'neighbour': |
| newa = _congrid_neighbor(a, newdims, m1, ofs) |
| |
| elif method in ['nearest', 'linear']: |
| # calculate new dims |
| for i in range(ndims): |
| base = np.arange(newdims[i]) |
| dimlist.append((old[i] - m1) / (newdims[i] - m1) * |
| (base + ofs) - ofs) |
| # specify old dims |
| olddims = [np.arange(i, dtype=np.float) for i in list(a.shape)] |
| |
| # first interpolation - for ndims = any |
| mint = scipy.interpolate.interp1d(olddims[-1], a, kind=method) |
| newa = mint(dimlist[-1]) |
| |
| trorder = [ndims - 1] + range(ndims - 1) |
| for i in range(ndims - 2, -1, -1): |
| newa = newa.transpose(trorder) |
| |
| mint = scipy.interpolate.interp1d(olddims[i], newa, kind=method) |
| newa = mint(dimlist[i]) |
| |
| if ndims > 1: |
| # need one more transpose to return to original dimensions |
| newa = newa.transpose(trorder) |
| |
| return newa |
| elif method in ['spline']: |
| nslices = [slice(0, j) for j in list(newdims)] |
| newcoords = np.mgrid[nslices] |
| |
| newcoords_dims = range(np.rank(newcoords)) |
| # make first index last |
| newcoords_dims.append(newcoords_dims.pop(0)) |
| newcoords_tr = newcoords.transpose(newcoords_dims) |
| # makes a view that affects newcoords |
| |
| newcoords_tr += ofs |
| |
| deltas = (np.asarray(old) - m1) / (newdims - m1) |
| newcoords_tr *= deltas |
| |
| newcoords_tr -= ofs |
| |
| newa = scipy.ndimage.map_coordinates(a, newcoords) |
| return newa |
| else: |
| print("Congrid error: Unrecognized interpolation type.\n", |
| "Currently only \'neighbour\', \'nearest\',\'linear\',", |
| "and \'spline\' are supported.") |
| return None |
| |
| |
| def _check_dataset_shapes(datasets): |
| """ If the datasets are not the same shape throw a ValueError Exception |
| |
| :param datasets: OCW Datasets to check for a consistent shape |
| :type datasets: List of OCW Dataset Objects |
| |
| :raises: ValueError |
| """ |
| dataset_shape = None |
| for dataset in datasets: |
| if dataset_shape is None: |
| dataset_shape = dataset.values.shape |
| else: |
| if dataset.values.shape != dataset_shape: |
| msg = "%s != %s" % (dataset.values.shape, dataset_shape) |
| raise ValueError("Input datasets must be the same shape " |
| "for an ensemble :: ", msg) |
| else: |
| pass |
| |
| |
| def _congrid_neighbor(values, new_dims, minus_one, offset): |
| """ Use the nearest neighbor to create a new array |
| |
| :param values: Array of values that need to be interpolated |
| :type values: Numpy ndarray |
| :param new_dims: Longitude resolution in degrees |
| :type new_dims: float |
| :param lat_resolution: Latitude resolution in degrees |
| :type lat_resolution: float |
| |
| :returns: A new spatially regridded Dataset |
| :rtype: Open Climate Workbench Dataset Object |
| """ |
| ndims = len(values.shape) |
| dimlist = [] |
| old_dims = np.array(values.shape) |
| for i in range(ndims): |
| base = np.indices(new_dims)[i] |
| dimlist.append((old_dims[i] - minus_one) / (new_dims[i] - minus_one) * |
| (base + offset) - offset) |
| cd = np.array(dimlist).round().astype(int) |
| new_values = values[list(cd)] |
| return new_values |
| |
| |
| def _are_bounds_contained_by_dataset(dataset, bounds): |
| '''Check if a Dataset fully contains a bounds. |
| |
| :param bounds: The Bounds object to check. |
| :type bounds: Bounds |
| :param dataset: The Dataset that should be fully contain the |
| Bounds |
| :type dataset: Dataset |
| |
| :returns: True if the Bounds are contained by the Dataset, Raises |
| a ValueError otherwise |
| ''' |
| lat_min, lat_max, lon_min, lon_max = dataset.spatial_boundaries() |
| start, end = dataset.temporal_boundaries() |
| |
| errors = [] |
| |
| # TODO: THIS IS TERRIBLY inefficent and we need to use a geometry |
| # lib instead in the future |
| if (lat_min > bounds.lat_max): |
| error = ("bounds.lat_max: %s is not between lat_min: %s and" |
| " lat_max: %s" % (bounds.lat_max, lat_min, lat_max)) |
| errors.append(error) |
| |
| if (lon_min > bounds.lon_max): |
| error = ("bounds.lon_max: %s is not between lon_min: %s and" |
| " lon_max: %s" % (bounds.lon_max, lon_min, lon_max)) |
| errors.append(error) |
| if not (bounds.start is None): |
| if not start <= bounds.start <= end: |
| error = ("bounds.start: %s is not between start: %s and end: %s" % |
| (bounds.start, start, end)) |
| errors.append(error) |
| |
| if not (bounds.end is None): |
| if not start <= bounds.end <= end: |
| error = ("bounds.end: %s is not between start: %s and end: %s" % |
| (bounds.end, start, end)) |
| errors.append(error) |
| |
| if len(errors) == 0: |
| return True |
| else: |
| error_message = '\n'.join(errors) |
| raise ValueError(error_message) |
| |
| |
| def _get_subregion_slice_indices(target_dataset, subregion): |
| '''Get the indices for slicing Dataset values to generate the subregion. |
| |
| :param subregion: The Bounds that specify the subset of the Dataset |
| that should be extracted. |
| :type subregion: Bounds |
| :param target_dataset: The Dataset to subset. |
| :type target_dataset: Dataset |
| |
| :returns: The indices to slice the Datasets arrays as a Dictionary. |
| ''' |
| latStart = min(np.nonzero(target_dataset.lats >= subregion.lat_min)[0]) |
| latEnd = max(np.nonzero(target_dataset.lats <= subregion.lat_max)[0]) |
| |
| lonStart = min(np.nonzero(target_dataset.lons >= subregion.lon_min)[0]) |
| lonEnd = max(np.nonzero(target_dataset.lons <= subregion.lon_max)[0]) |
| |
| if not (subregion.start is None): |
| timeStart = min(np.nonzero(target_dataset.times >= subregion.start)[0]) |
| timeEnd = max(np.nonzero(target_dataset.times <= subregion.end)[0]) |
| else: |
| timeStart = 0 |
| timeEnd = len(target_dataset.times) -1 |
| |
| return { |
| "lat_start": latStart, |
| "lat_end": latEnd, |
| "lon_start": lonStart, |
| "lon_end": lonEnd, |
| "time_start": timeStart, |
| "time_end": timeEnd |
| } |
