#
#  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.
#

from ocw import dataset as ds

import datetime
import numpy as np
import numpy.ma as ma
import scipy.interpolate
import scipy.ndimage
from scipy.ndimage import map_coordinates
import netCDF4

import logging

logger = logging.getLogger(__name__)

def temporal_subset(month_start, month_end, target_dataset, 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 = range(month_start,13)
        month_index.extend(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 = np.append(time_index, np.where(months == m_value)[0])
        if m_value == month_index[0]:
            time_index_first = np.min(np.where(months == m_value)[0])
        if m_value == month_index[-1]:
            time_index_last = np.max(np.where(months == m_value)[0])

    time_index = np.sort(time_index)

    time_index = time_index[np.where((time_index >= time_index_first) & (time_index <= time_index_last))]

    time_index = list(time_index)

    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:
        nmonth = len(month_index)
        ntime  = new_dataset.times.size
        nyear = ntime/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_dataset.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 bin size
    :type temporal_resolution: :class:`datetime.timedelta`
    
    :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
    day_count = temporal_resolution.days
    time_unit = None
    if day_count == 1:
        time_unit = 'daily'
    elif day_count > 1 and day_count <= 31:
        time_unit = 'monthly'
    elif day_count > 31 and day_count <= 366:
        time_unit = 'annual'
    else:
        time_unit = 'full'

    masked_values = target_dataset.values.view(ma.MaskedArray)
    binned_values, binned_dates = _rcmes_calc_average_on_new_time_unit_K(masked_values, target_dataset.times, time_unit)
    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 spatial_regrid(target_dataset, new_latitudes, new_longitudes):
    """ 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`

    :returns: A new spatially regridded Dataset
    :rtype: :class:`dataset.Dataset`
    """
    # Make masked array of shape (times, new_latitudes,new_longitudes)
    new_values = ma.zeros([len(target_dataset.times), 
                           len(new_latitudes), 
                           len(new_longitudes)])

    # 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)
    lons, lats = np.meshgrid(target_dataset.lons, target_dataset.lats)
    new_lons, new_lats = np.meshgrid(new_longitudes, new_latitudes)
    # 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)):
        new_values[i] = _rcmes_spatial_regrid(target_values[i],
                                              lats,
                                              lons,
                                              new_lats,
                                              new_lons)
    
    # 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(subregion, target_dataset, subregion_name=None):
    '''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`

    :returns: The subset-ed Dataset object
    :rtype: :class:`dataset.Dataset`

    :raises: ValueError
    '''

    if not subregion.start:
        subregion.start = target_dataset.times[0] 
        subregion.end = target_dataset.times[-1]
        
    # Ensure that the subregion information is well formed
    _are_bounds_contained_by_dataset(subregion, target_dataset)

    # Get subregion indices into subregion data
    dataset_slices = _get_subregion_slice_indices(subregion, target_dataset)

    if not subregion_name:
        subregion_name = target_dataset.name

   # 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
    )


def safe_subset(subregion, target_dataset, 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.time_range()

    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(subregion, target_dataset, 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 lenghts
    lat_len = len(dataset.lats)
    lon_len = len(dataset.lons)
    time_len = len(dataset.times)

    # Create attribute dimensions
    lat_dim = out_file.createDimension('lat', lat_len)
    lon_dim = out_file.createDimension('lon', lon_len)
    time_dim = out_file.createDimension('time', time_len)

    # Create variables
    lats = out_file.createVariable('lat', 'f8', ('lat',), zlib=compress)
    lons = out_file.createVariable('lon', 'f8', ('lon',), 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', 'lat', 'lon'),
                                    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 = len(dataset.lats)
    lon_len = len(dataset.lons)
    time_len = len(dataset.times)

    if not subregions == None:
        nsubregion = len(subregions)
                 
    # Create attribute dimensions
    lat_dim = out_file.createDimension('y', lat_len)
    lon_dim = out_file.createDimension('x', lon_len)
    time_dim = out_file.createDimension('time', time_len)

    # Create variables and store the values
    lats = out_file.createVariable('lat', 'f8', ('y'))
    lats[:] = dataset.lats
    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 not subregions == 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']
    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':
            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:
        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 _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]
    #print nlat, nlon, "lats, lons - incoming dataset"

    nlat2 = lat2.shape[0]
    nlon2 = lon2.shape[1]
    #print nlat2, nlon2, "NEW lats, lons - for the new grid output"

    # 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[loni > lon.max()] = lon.max()
    loni[loni < lon.min()] = lon.min()
    
    # To deal with the "hard" break, we'll have to treat y differently,
    # so we're just setting the min here...
    lati[lati > lat.max()] = lat.max()
    lati[lati < lat.min()] = 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
            spatial_values.data[idx] = q_shifted[idx]

    # 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)
    qmdi[spatial_values.mask == True] = 1.
    qmdi[spatial_values.mask == 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_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
    """

    # Check if the user-selected temporal grid is valid. If not, EXIT
    acceptable = (unit=='full')|(unit=='annual')|(unit=='monthly')|(unit=='daily')
    if not acceptable:
        print 'Error: unknown unit type selected for time averaging: EXIT'
        return -1,-1,-1,-1

    # 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)):
            timeunits.append(999)  # i.e. we just want the same value for all times.
        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 not a.dtype 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']:
        oslices = [ slice(0, j) for j in old ]
        oldcoords = np.ogrid[oslices]
        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 == 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(bounds, dataset):
    '''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.time_range()
    errors = []

    # TODO:  THIS IS TERRIBLY inefficent and we need to use a geometry lib instead in the future
    if not lat_min <= bounds.lat_min <= lat_max:
        error = "bounds.lat_min: %s is not between lat_min: %s and lat_max: %s" % (bounds.lat_min, lat_min, lat_max)
        errors.append(error)

    if not lat_min <= bounds.lat_max <= 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 not lon_min <= bounds.lon_min <= lon_max:
        error = "bounds.lon_min: %s is not between lon_min: %s and lon_max: %s" % (bounds.lon_min, lon_min, lon_max)
        errors.append(error)

    if not lon_min <= bounds.lon_max <= 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 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 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(subregion, target_dataset):
    '''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])


    timeStart = min(np.nonzero(target_dataset.times >= subregion.start)[0])
    timeEnd = max(np.nonzero(target_dataset.times <= subregion.end)[0])

    return {
        "lat_start"  : latStart,
        "lat_end"    : latEnd,
        "lon_start"  : lonStart,
        "lon_end"    : lonEnd,
        "time_start" : timeStart,
        "time_end"   : timeEnd
    }