| # 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 sys |
| import os |
| import datetime as dt |
| import numpy as np |
| import numpy.ma as ma |
| |
| from mpl_toolkits.basemap import shiftgrid, Basemap |
| from matplotlib.path import Path |
| from dateutil.relativedelta import relativedelta |
| from netCDF4 import num2date, date2num |
| import scipy.interpolate as interpolate |
| import scipy.stats as stats |
| from scipy.ndimage import map_coordinates |
| |
| |
| def decode_time_values(dataset, time_var_name): |
| ''' Decode NetCDF time values into Python datetime objects. |
| |
| :param dataset: The dataset from which time values should be extracted. |
| :type dataset: netCDF4.Dataset |
| :param time_var_name: The name of the time variable in dataset. |
| :type time_var_name: :mod:`string` |
| |
| :returns: The list of converted datetime values. |
| |
| :raises ValueError: If the time units value couldn't be parsed, if the |
| base time value couldn't be parsed, or if the time_var_name could not |
| be found in the dataset. |
| ''' |
| time_data = dataset.variables[time_var_name] |
| time_format = time_data.units |
| if time_format[-1].lower() == 'z': |
| time_format = time_format[:-1] |
| if time_format[-3:].lower() == 'utc': |
| time_format = time_format[:-3] |
| |
| time_units = parse_time_units(time_format) |
| |
| times = [] |
| if time_units == 'months': |
| # datetime.timedelta doesn't support a 'months' option. To remedy |
| # this, a month == 30 days for our purposes. |
| time_base = parse_time_base(time_format) |
| for time_val in time_data: |
| times.append(time_base + relativedelta(months=int(time_val))) |
| else: |
| try: |
| times_calendar = time_data.calendar |
| except: |
| times_calendar = 'standard' |
| |
| times = num2date( |
| time_data[:], units=time_format, calendar=times_calendar) |
| if times_calendar == '360_day': |
| for it, t in enumerate(times): |
| times[it] = dt.datetime.strptime(str(t), '%Y-%m-%d 00:00:00') |
| return times |
| |
| |
| def parse_time_units(time_format): |
| ''' Parse units value from time units string. |
| |
| The only units that are supported are: seconds, minutes, hours, days, |
| months, or years. |
| |
| :param time_format: The time data units string from the dataset |
| being processed. The string should be of the format |
| '<units> since <base time date>' |
| :type time_format: :mod:`string` |
| |
| :returns: The unit substring from the time units string |
| |
| :raises ValueError: If the units present in the time units string doesn't |
| match one of the supported unit value. |
| ''' |
| for unit in ['seconds', 'minutes', 'hours', 'days', 'months', 'years']: |
| if unit in time_format: |
| return unit |
| else: |
| cur_frame = sys._getframe().f_code |
| err = "{}.{}: Unable to parse valid units from {}".format( |
| cur_frame.co_filename, |
| cur_frame.co_name, |
| time_format |
| ) |
| raise ValueError(err) |
| |
| |
| def parse_time_base(time_format): |
| ''' Parse time base object from the time units string. |
| |
| :param time_format: The time data units string from the dataset |
| being processed. The string should be of the format |
| '<units> since <base time date>' |
| :type time_format: :mod:`string` |
| |
| :returns: The base time as a datetime object. |
| |
| :raises ValueError: When the base time string couldn't be parsed from the |
| units time_format string or if the date string didn't match any of the |
| expected formats. |
| ''' |
| base_time_string = parse_base_time_string(time_format) |
| |
| time_format = time_format.strip() |
| |
| possible_time_formats = [ |
| '%Y:%m:%d %H:%M:%S', '%Y-%m-%d %H-%M-%S', '%Y/%m/%d %H/%M/%S', |
| '%Y-%m-%d %H:%M:%S', '%Y/%m/%d %H:%M:%S', '%Y%m%d %H:%M:%S', |
| '%Y%m%d%H%M%S', '%Y-%m-%d-%H-%M-%S', '%Y/%m/%d/%H/%M/%S', |
| '%Y:%m:%d:%H:%M:%S', '%Y-%m-%d-%H:%M:%S', '%Y-%m-%d %H:%M:%S', |
| '%Y/%m/%d%H:%M:%S', '%Y-%m-%d %H:%M', '%Y/%m/%d %H:%M', |
| '%Y:%m:%d %H:%M', '%Y%m%d %H:%M', '%Y-%m-%d', '%Y/%m/%d', |
| '%Y:%m:%d', '%Y%m%d', '%Y-%m-%d %H:%M:%S.%f', '%Y-%m-%d %H', |
| '%Y-%m-%dT%H:%M:%S', '%Y-%m-%dT%H:%M:%SZ' |
| ] |
| |
| # Attempt to match the base time string with a possible format parsing |
| # string. |
| for time_format in possible_time_formats: |
| try: |
| stripped_time = dt.datetime.strptime(base_time_string, time_format) |
| break |
| except ValueError: |
| # This exception means that the time format attempted was |
| # incorrect. No need to report or raise this, |
| # simply try the next one! |
| pass |
| # If we got through the entire loop without a break, we couldn't parse the |
| # date string with our known formats. |
| else: |
| cur_frame = sys._getframe().f_code |
| err = "{}.{}: Unable to parse valid date from {}".format( |
| cur_frame.co_filename, |
| cur_frame.co_name, |
| base_time_string |
| ) |
| |
| raise ValueError(err) |
| |
| return stripped_time |
| |
| |
| def parse_base_time_string(time_format): |
| ''' Retrieve base time string from time data units information. |
| |
| :param time_format: The time data units string from the dataset |
| being processed. The string should be of the format |
| '<units> since <base time date>' |
| :type time_format: :mod:`string` |
| |
| :returns: The base time string split out of the time units information. |
| |
| :raises ValueError: When the time_format parameter is malformed. |
| ''' |
| if 'since' not in time_format: |
| cur_frame = sys._getframe().f_code |
| err = "{}.{}: Invalid time_format value {} given".format( |
| cur_frame.co_filename, |
| cur_frame.co_name, |
| time_format |
| ) |
| |
| raise ValueError(err) |
| |
| return time_format.split('since')[1].strip() |
| |
| |
| def normalize_lat_lon_values(lats, lons, values): |
| ''' Normalize lat/lon values |
| |
| Ensure that lat/lon values are within [-180, 180)/[-90, 90) as well |
| as sorted. If the values are off the grid they are shifted into the |
| expected range. |
| |
| :param lats: A 1D numpy array of sorted lat values. |
| :type lats: :class:`numpy.ndarray` |
| :param lons: A 1D numpy array of sorted lon values. |
| :type lons: :class:`numpy.ndarray` |
| :param values: A 3D array of data values. |
| |
| :returns: A :func:`tuple` of the form (adjust_lats, adjusted_lons, |
| adjusted_values) |
| |
| :raises ValueError: If the lat/lon values are not sorted. |
| ''' |
| if lats.ndim == 1 and lons.ndim == 1: |
| # Avoid unnecessary shifting if all lons are higher than 180 |
| if lons.min() > 180: |
| lons -= 360 |
| |
| # Make sure lats and lons are monotonically increasing |
| lats_decreasing = np.diff(lats) < 0 |
| lons_decreasing = np.diff(lons) < 0 |
| |
| # If all values are decreasing then they just need to be reversed |
| lats_reversed, lons_reversed = (lats_decreasing.all(), |
| lons_decreasing.all()) |
| |
| # If the lat values are unsorted then raise an exception |
| if not lats_reversed and lats_decreasing.any(): |
| raise ValueError('Latitudes must be sorted.') |
| |
| # Perform same checks now for lons |
| if not lons_reversed and lons_decreasing.any(): |
| raise ValueError('Longitudes must be sorted.') |
| |
| # Also check if lons go from [0, 360), and convert to [-180, 180) |
| # if necessary |
| lons_shifted = lons.max() > 180 |
| lats_out, lons_out, data_out = lats[:], lons[:], values[:] |
| # Now correct data if latlon grid needs to be shifted |
| if lats_reversed: |
| lats_out = lats_out[::-1] |
| data_out = data_out[..., ::-1, :] |
| |
| if lons_reversed: |
| lons_out = lons_out[::-1] |
| data_out = data_out[..., :, ::-1] |
| |
| if lons_shifted: |
| data_out, lons_out = shiftgrid( |
| 180, data_out, lons_out, start=False) |
| |
| return lats_out, lons_out, data_out |
| else: |
| lons[lons > 180] = lons[lons > 180] - 360. |
| |
| return lats, lons, values |
| |
| |
| def reshape_monthly_to_annually(dataset): |
| ''' Reshape monthly binned dataset to annual bins. |
| |
| Reshape a monthly binned dataset's 3D value array with shape |
| (num_months, num_lats, num_lons) to a 4D array with shape |
| (num_years, 12, num_lats, num_lons). This causes the data to be binned |
| annually while retaining its original shape. |
| |
| It is assumed that the number of months in the dataset is evenly |
| divisible by 12. If it is not you will receive error due to |
| an invalid shape. |
| |
| Example change of a dataset's shape: |
| (24, 90, 180) -> (2, 12, 90, 180) |
| |
| :param dataset: Dataset object with full-year format |
| :type dataset: :class:`dataset.Dataset` |
| |
| :returns: Dataset values array with shape (num_year, 12, num_lat, num_lon) |
| ''' |
| |
| values = dataset.values[:] |
| data_shape = values.shape |
| num_total_month = data_shape[0] |
| num_year = num_total_month // 12 |
| if num_total_month % 12 != 0: |
| raise ValueError("Number of months in dataset ({}) does not " |
| "divide evenly into {} year(s)." |
| .format(num_total_month, num_year)) |
| |
| num_month = 12 |
| year_month_shape = num_year, num_month |
| lat_lon_shape = data_shape[1:] |
| # Make new shape (num_year, 12, num_lats, num_lons) |
| new_shape = tuple(year_month_shape + lat_lon_shape) |
| # Reshape data with new shape |
| values.shape = new_shape |
| |
| return values |
| |
| |
| def calc_temporal_mean(dataset): |
| ''' Calculate temporal mean of dataset's values |
| |
| :param dataset: OCW Dataset whose first dimension is time |
| :type dataset: :class:`dataset.Dataset` |
| |
| :returns: Mean values averaged for the first dimension (time) |
| ''' |
| return ma.mean(dataset.values, axis=0) |
| |
| |
| def calc_climatology_year(dataset): |
| ''' Calculate climatology of dataset's values for each year |
| |
| :param dataset: Monthly binned Dataset object with an evenly divisible |
| number of months. |
| :type dataset: :class:`dataset.Dataset` |
| |
| :returns: Mean values for each year (annual_mean) and mean values for all |
| years (total_mean) |
| |
| :raise ValueError: If the number of monthly bins is not evenly divisible |
| by 12. |
| ''' |
| |
| values_shape = dataset.values.shape |
| time_shape = values_shape[0] |
| if time_shape % 12: |
| raise ValueError('The dataset should be in full-time format.') |
| else: |
| # Get values reshaped to (num_year, 12, num_lats, num_lons) |
| values = reshape_monthly_to_annually(dataset) |
| # Calculate mean values over year (num_year, num_lats, num_lons) |
| annually_mean = values.mean(axis=1) |
| # Calculate mean values over all years (num_lats, num_lons) |
| total_mean = annually_mean.mean(axis=0) |
| |
| return annually_mean, total_mean |
| |
| |
| def calc_climatology_monthly(dataset): |
| ''' Calculate monthly mean values for a dataset. |
| Follow COARDS climo stats calculation, the year can be given as 0 |
| but the min year allowed in Python is 1 |
| http://www.cgd.ucar.edu/cms/eaton/netcdf/CF-20010629.htm#climatology |
| |
| :param dataset: Monthly binned Dataset object with the number of months |
| divisible by 12 |
| :type dataset: :class:`dataset.Dataset` |
| |
| :returns: Mean values for each month of the year of shape |
| (12, num_lats, num_lons) and times array of datetime objects |
| of length 12 |
| |
| :raise ValueError: If the number of monthly bins is not divisible by 12 |
| ''' |
| |
| if dataset.values.shape[0] % 12: |
| error = ( |
| "The length of the time axis in the values array should be " |
| "divisible by 12." |
| ) |
| raise ValueError(error) |
| else: |
| values = reshape_monthly_to_annually(dataset).mean(axis=0) |
| |
| # A year can commence from any month |
| first_month = dataset.times[0].month |
| times = np.array([dt.datetime(1, first_month, 1) + |
| relativedelta(months=x) |
| for x in range(12)]) |
| return values, times |
| |
| |
| def calc_time_series(dataset): |
| ''' Calculate time series mean values for a dataset |
| |
| :param dataset: Dataset object |
| :type dataset: :class:`dataset.Dataset` |
| |
| :returns: time series for the dataset of shape (nT) |
| ''' |
| |
| t_series = [] |
| for t in range(dataset.values.shape[0]): |
| t_series.append(dataset.values[t, :, :].mean()) |
| |
| return t_series |
| |
| |
| def get_temporal_overlap(dataset_array): |
| ''' Find the maximum temporal overlap across the observation and model datasets |
| |
| :param dataset_array: an array of OCW datasets |
| ''' |
| start_times = [] |
| end_times = [] |
| for dataset in dataset_array: |
| slc = trim_dataset(dataset) |
| times = dataset.times[slc] |
| start_time, end_time = times[0], times[-1] |
| start_times.append(start_time) |
| end_times.append(end_time) |
| |
| return np.max(start_times), np.min(end_times) |
| |
| def adjust_model_years_for_climatology_calculation(dataset_array): |
| ''' Using the time length of the first element in the input dataset_array, |
| adjust years in the rest ofi the dataset_array so that every dataset ends in the same year. |
| :param dataset_array: an array of OCW datasets |
| ''' |
| slc = trim_dataset(dataset_array[0]) |
| obs_times = dataset_array[0].times[slc] |
| for idata, dataset in enumerate(dataset_array[1:]): |
| year_diff = obs_times[-1].year - dataset.times[-1].year |
| nt = dataset.times.size |
| for it in np.arange(nt): |
| dataset.times[it] = dataset.times[it].replace(year = dataset.times[it].year + year_diff) |
| return dataset_array |
| |
| def trim_dataset(dataset): |
| ''' Trim datasets such that first and last year of data have all 12 months |
| |
| :param dataset: Dataset object |
| :type dataset: :class:`dataset.Dataset` |
| |
| :returns: Slice index for trimmed dataset |
| ''' |
| start_time, end_time = dataset.temporal_boundaries() |
| start_month = 13 if start_time.month == 1 else start_time.month |
| end_month = 0 if end_time.month == 12 else end_time.month |
| slc = slice(13 - start_month, len(dataset.times) - end_month) |
| return slc |
| |
| |
| def calc_subregion_area_mean_and_std(dataset_array, subregions): |
| ''' Calculate area mean and standard deviation values for a given \ |
| subregions using datasets on common grid points |
| |
| :param dataset_array: An array of OCW Dataset Objects \ |
| :type list: :mod:'list' |
| |
| :param subregions: list of subregions \ |
| :type subregions: :class:`numpy.ma.array` |
| |
| :returns: area averaged time series for the dataset of shape \ |
| (ntime, nsubregion) |
| ''' |
| |
| ndata = len(dataset_array) |
| dataset0 = dataset_array[0] |
| if dataset0.lons.ndim == 1: |
| lons, lats = np.meshgrid(dataset0.lons, dataset0.lats) |
| else: |
| lons = dataset0.lons |
| lats = dataset0.lats |
| subregion_array = np.zeros(lons.shape) |
| mask_array = dataset_array[0].values[0, :].mask |
| # dataset0.values.shsape[0]: length of the time dimension |
| # spatial average |
| t_series = ma.zeros([ndata, dataset0.values.shape[0], len(subregions)]) |
| # spatial standard deviation |
| spatial_std = ma.zeros([ndata, dataset0.values.shape[0], len(subregions)]) |
| |
| for iregion, subregion in enumerate(subregions): |
| lat_min, lat_max, lon_min, lon_max = subregion[1] |
| y_index, x_index = np.where((lats >= lat_min) & ( |
| lats <= lat_max) & (lons >= lon_min) & (lons <= lon_max)) |
| subregion_array[y_index, x_index] = iregion + 1 |
| for idata in np.arange(ndata): |
| t_series[idata, :, iregion] = ma.mean(dataset_array[idata].values[ |
| :, y_index, x_index], axis=1) |
| spatial_std[idata, :, iregion] = ma.std( |
| dataset_array[idata].values[:, y_index, x_index], axis=1) |
| subregion_array = ma.array(subregion_array, mask=mask_array) |
| return t_series, spatial_std, subregion_array |
| |
| |
| def calc_area_weighted_spatial_average(dataset, area_weight=False): |
| '''Calculate area weighted average of the values in OCW dataset |
| |
| :param dataset: Dataset object |
| :type dataset: :class:`dataset.Dataset` |
| |
| :returns: time series for the dataset of shape (nT) |
| ''' |
| |
| if dataset.lats.ndim == 1: |
| lons, lats = np.meshgrid(dataset.lons, dataset.lats) |
| else: |
| lats = dataset.lats |
| weights = np.cos(lats * np.pi / 180.) |
| |
| nt, ny, nx = dataset.values.shape |
| spatial_average = ma.zeros(nt) |
| for it in np.arange(nt): |
| if area_weight: |
| spatial_average[it] = ma.average( |
| dataset.values[it, :], weights=weights) |
| else: |
| spatial_average[it] = ma.average(dataset.values[it, :]) |
| |
| return spatial_average |
| |
| |
| def shapefile_boundary(boundary_type, region_names): |
| ''' |
| :param boundary_type: The type of spatial subset boundary |
| :type boundary_type: :mod:'string' |
| |
| :param region_names: An array of regions for spatial subset |
| :type region_names: :mod:'list' |
| ''' |
| # Read the shapefile |
| map_read = Basemap() |
| regions = [] |
| shapefile_dir = os.sep.join([os.path.dirname(__file__), 'shape']) |
| map_read.readshapefile(os.path.join(shapefile_dir, boundary_type), |
| boundary_type, drawbounds=False) |
| # Note: The shapefile may contain countries or states with latin letters. |
| # Hence, the output of readshapefile can be a mix of ascii and unicode |
| # strings. Because python 2 and 3 treat them differently, we must |
| # explicitly check. |
| if boundary_type == 'us_states': |
| for region_name in region_names: |
| region_name = _force_unicode(region_name) |
| for iregion, region_info in enumerate(map_read.us_states_info): |
| state1 = _force_unicode(region_info['st'], 'latin-1') |
| state2 = _force_unicode(region_info['state'], 'latin-1') |
| if state1 == region_name or state2 == region_name: |
| regions.append(np.array(map_read.us_states[iregion])) |
| elif boundary_type == 'countries': |
| for region_name in region_names: |
| region_name = _force_unicode(region_name) |
| for iregion, region_info in enumerate(map_read.countries_info): |
| country = _force_unicode(region_info['COUNTRY'], 'latin-1') |
| if (country.replace(" ", "").lower() == |
| region_name.replace(" ", "").lower()): |
| regions.append(np.array(map_read.countries[iregion])) |
| return regions |
| |
| |
| def CORDEX_boundary(domain_name): |
| ''' |
| :param domain_name: CORDEX domain name (http://www.cordex.org/) |
| :type domain_name: :mod:'string' |
| ''' |
| if domain_name == 'southamerica' or domain_name == 'sam': |
| return -57.61, 18.50, 254.28 - 360., 343.02 - 360. |
| elif domain_name == 'centralamerica' or domain_name == 'cam': |
| return -19.46, 34.83, 235.74 - 360., 337.78 - 360. |
| elif domain_name == 'northamerica' or domain_name == 'nam': |
| return 12.55, 75.88, 189.26 - 360., 336.74 - 360. |
| elif domain_name == 'europe' or domain_name == 'eur': |
| return 22.20, 71.84, 338.23 - 360., 64.4 |
| elif domain_name == 'africa' or domain_name == 'afr': |
| return -45.76, 42.24, 335.36 - 360., 60.28 |
| elif domain_name == 'southasia' or domain_name == 'was': |
| return -15.23, 45.07, 19.88, 115.55 |
| elif domain_name == 'eastasia' or domain_name == 'eas': |
| return -0.10, 61.90, 51.59, 179.99 |
| elif domain_name == 'centralasia' or domain_name == 'cas': |
| return 18.34, 69.37, 11.05, 139.13 |
| elif domain_name == 'australasia' or domain_name == 'aus': |
| return -52.36, 12.21, 89.25, 179.99 |
| elif domain_name == 'antartica' or domain_name == 'ant': |
| return -89.48, -56.00, -179.00, 179.00 |
| elif domain_name == 'artic' or domain_name == 'arc': |
| return 46.06, 89.50, -179.00, 179.00 |
| elif domain_name == 'mediterranean' or domain_name == 'med': |
| return 25.63, 56.66, 339.79 - 360.00, 50.85 |
| elif domain_name == 'middleeastnorthafrica' or domain_name == 'mna': |
| return -7.00, 45.00, 333.00 - 360.00, 76.00 |
| elif domain_name == 'southeastasia' or domain_name == 'sea': |
| return -15.14, 27.26, 89.26, 146.96 |
| else: |
| err = "Invalid CORDEX domain name" |
| raise ValueError(err) |
| |
| |
| def mask_using_shapefile_info(lons, lats, masked_regions, extract=True): |
| if lons.ndim == 2 and lats.ndim == 2: |
| lons_2d = lons |
| lats_2d = lats |
| else: |
| lons_2d, lats_2d = np.meshgrid(lons, lats) |
| |
| points = np.array((lons_2d.flatten(), lats_2d.flatten())).T |
| for iregion, region in enumerate(masked_regions): |
| mpath = Path(region) |
| mask0 = mpath.contains_points(points).reshape(lons_2d.shape) |
| if iregion == 0: |
| mask = mask0 |
| else: |
| 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) |
| |
| 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) |
| |
| 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] |
| new_data_array = ma.zeros(data_array.shape) |
| if isinstance(data_array.mask, bool): |
| new_mask = mask |
| for it in np.arange(nt): |
| new_data_array[it,:] = ma.array(data_array[it,:], |
| mask=new_mask) |
| else: |
| for it in np.arange(nt): |
| new_mask = data_array[it, :].mask | mask |
| new_data_array[it,:] = ma.array(data_array[it,:], |
| mask=new_mask) |
| |
| return new_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 |
| |
| def _force_unicode(s, encoding='utf-8'): |
| ''' |
| If the input is bytes, convert to unicode, otherwise return the input |
| ''' |
| if hasattr(s, 'decode'): |
| s = s.decode(encoding=encoding) |
| |
| return s |
| |
| def calculate_temporal_trends(dataset): |
| ''' Calculate temporal trends in dataset.values |
| :param dataset: The dataset from which time values should be extracted. |
| :type dataset: :class:`dataset.Dataset` |
| |
| :returns: Arrays of the temporal trend and standard error |
| :rtype: :class:`numpy.ma.core.MaskedArray` |
| ''' |
| |
| nt, ny, nx = dataset.values.shape |
| x = np.arange(nt) |
| |
| trend = np.zeros([ny, nx])-999. |
| slope_err = np.zeros([ny, nx])-999. |
| for iy in np.arange(ny): |
| for ix in np.arange(nx): |
| if dataset.values[:,iy,ix].count() == nt: |
| trend[iy,ix], slope_err[iy,ix] = calculate_temporal_trend_of_time_series( |
| x, dataset.values[:,iy,ix]) |
| |
| return ma.masked_equal(trend, -999.), ma.masked_equal(slope_err, -999.) |
| |
| def calculate_ensemble_temporal_trends(timeseries_array, number_of_samples=1000): |
| ''' Calculate temporal trends in an ensemble of time series |
| :param timeseries_array: Two dimensional array. 1st index: model, 2nd index: time. |
| :type timeseries_array: :class:`numpy.ndarray` |
| |
| :param sampling: A list whose elements are one-dimensional numpy arrays |
| :type timeseries_array: :class:`list` |
| |
| :returns: temporal trend and estimated error from bootstrapping |
| :rtype: :class:`float`, :class:`float` |
| ''' |
| |
| nmodels, nt = timeseries_array.shape |
| x = np.arange(nt) |
| sampled_trend = np.zeros(number_of_samples) |
| ensemble_trend, _ = calculate_temporal_trend_of_time_series( |
| x, np.mean(timeseries_array, axis=0)) |
| |
| for isample in np.arange(number_of_samples): |
| index = np.random.choice(nmodels, size=nmodels, replace=True) |
| random_ensemble = np.mean(timeseries_array[index, :], axis=0) |
| sampled_trend[isample], _ = calculate_temporal_trend_of_time_series( |
| x, random_ensemble) |
| return ensemble_trend, np.std(sampled_trend, ddof=1) |
| |
| |
| def calculate_temporal_trend_of_time_series(x,y): |
| ''' Calculate least-square trends (a) in y = ax+b and a's standard error |
| :param x: time series |
| :type x: :class:`numpy.ndarray` |
| |
| :param x: time series |
| :type x: :class:`numpy.ndarray` |
| |
| :returns: temporal trend and standard error |
| :rtype: :class:`float`, :class:`float` |
| ''' |
| slope, intercept, r_value, p_value, std_err = stats.linregress(x,y) |
| return slope, std_err |
| |
| def calculate_daily_climatology(dataset): |
| '''Calculate daily climatology from the input dataset |
| :param dataset: The dataset to convert. |
| :type dataset: :class:`dataset.Dataset` |
| :returns: values_clim |
| :rtype: 3d masked numpy array.shape (number of unique days, y, x) |
| ''' |
| |
| 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) |
| return values_clim |
