Return the indices to crop dataset

        Assuming that the dataset have the dimensions given by
          dims, this function return the indices to conform with
          the given coordinates (lat, lon, ...)
187618
def cropIndices(dims, lat, lon, depth=None, doy=None):
    
    dims_out = {}
    idx = {}
    
    yn = slice(
            np.nonzero(dims['lat'] <= lat.min())[0].max(),
            np.nonzero(dims['lat'] >= lat.max())[0].min() + 1)
    dims_out['lat'] = np.atleast_1d(dims['lat'][yn])
    idx['yn'] = yn

    lon_ext = np.array(
            (dims['lon'] - 2*360).tolist() +
            (dims['lon'] - 360).tolist() +
            dims['lon'].tolist() +
            (dims['lon'] + 360).tolist())
    xn_ext = list(4 * list(range(dims['lon'].shape[0])))
    xn_start = np.nonzero(lon_ext <= lon.min())[0].max()
    xn_end = np.nonzero(lon_ext >= lon.max())[0].min()
    xn = xn_ext[xn_start:xn_end+1]
    dims_out['lon'] = np.atleast_1d(lon_ext[xn_start:xn_end+1])
    idx['xn'] = xn

    if depth is not None:
        zn = slice(
                np.nonzero(dims['depth'] <= depth.min())[0].max(),
                np.nonzero(dims['depth'] >= min(dims['depth'].max(), depth.max())
                    )[0].min() + 1
                                                        )
        # If a higher degree interpolation system uses more than one data
        #   point in the edge, I should extend this selection one point on
        #   each side, without go beyond 0
        # if zn.start < 0:
        #    zn = slice(0, zn.stop, zn.step)
        dims_out['depth'] = np.atleast_1d(dims['depth'][zn])
        idx['zn'] = zn

    if doy is not None:
        # Source has only one time, like total mean field, or annual mean.
        if dims['time'].shape == (1,):
            dims_out['time'] = dims['time']
            idx['tn'] = [0]
        else:
            time_ext = np.array(
                    [dims['time'][-1] - 365.25] +
                            dims['time'].tolist() +
                            [dims['time'][0] + 365.25])
            tn_ext = list(range(dims['time'].size))
            tn_ext = [tn_ext[-1]] + tn_ext + [tn_ext[0]]
            tn_start = np.nonzero(time_ext <= doy.min())[0].max()
            tn_end = np.nonzero(time_ext >= doy.max())[0].min()
            dims_out['time'] = np.atleast_1d(time_ext[tn_start:tn_end+1])
            idx['tn'] = tn_ext[tn_start:tn_end+1]

    return dims_out, idx
553310
def indexes(self):
        
        if self._indexes is None:
            self._indexes = default_indexes(self._coords, self.dims)
        return Indexes(self._indexes)
1307156
def _fetch_dimensions(self, dataset):
        
        for dimension in super(SCB, self)._fetch_dimensions(dataset):
            if dimension.id == "Region":
                yield Dimension(dimension.id,
                                datatype="region",
                                dialect="skatteverket",
                                label=dimension.label)
            else:
                yield dimension
1139393
def nn_dims(self, x, y, dims_x, dims_y, k=1, radius=np.inf, eps=0.0, p=2):
        
        if self.size > 0:
            dists, idxes = Dataset.nn_dims(self, x, y, dims_x, dims_y, k, radius, eps, p)
        else:
            return self.buffer.nn_dims(x, y, dims_x, dims_y, k, radius, eps, p)
        if self.buffer.size > 0:
            buffer_dists, buffer_idxes = self.buffer.nn_dims(x, y, dims_x, dims_y, k, radius, eps, p)
            buffer_idxes = [i + self.size for i in buffer_idxes]
            ziped = zip(dists, idxes)
            buffer_ziped = zip(buffer_dists, buffer_idxes)
            sorted_dists_idxes = sorted(ziped + buffer_ziped, key=lambda di:di[0])
            knns = sorted_dists_idxes[:k]
            return [knn[0] for knn in knns], [knn[1] for knn in knns] 
        else:
            return dists, idxes
245225
def area2lonlat(dataarray):
    """"""
    area = dataarray.attrs['area']
    lons, lats = area.get_lonlats_dask()
    lons = xr.DataArray(lons, dims=['y', 'x'],
                        attrs={'name': "longitude",
                               'standard_name': "longitude",
                               'units': 'degrees_east'},
                        name='longitude')
    lats = xr.DataArray(lats, dims=['y', 'x'],
                        attrs={'name': "latitude",
                               'standard_name': "latitude",
                               'units': 'degrees_north'},
                        name='latitude')
    dataarray.attrs['coordinates'] = 'longitude latitude'

    return [dataarray, lons, lats]
540525
def load_coordinates(self):
        
        coord_file = Dataset(self.coordinate_file)
        if "lon" in coord_file.variables.keys():
            self.coordinates["lon"] = coord_file.variables["lon"][:]
            self.coordinates["lat"] = coord_file.variables["lat"][:]
        else:
            self.coordinates["lon"] = coord_file.variables["XLONG"][0]
            self.coordinates["lat"] = coord_file.variables["XLAT"][0]
        coord_file.close()
1303691
def get_dimension_indices(self, query):
        
        ids = self['id'] if self.get('id') else self['dimension']['id']
        indices = []

        for idx, id in enumerate(ids):
            indices.append(self.get_dimension_index(id,
                                                    [d.get(id) for d in query
                                                     if id in d][0]))

        return indices
1076794
def olindices(order, dim):
    
    indices = [olindex(o, dim) for o in range(order+1)]
    indices = numpy.vstack(indices)
    return indices
1309108
def assert_coordinate_consistent(obj, coords):
    
    for k in obj.dims:
        # make sure there are no conflict in dimension coordinates
        if k in coords and k in obj.coords:
            if not coords[k].equals(obj[k].variable):
                raise IndexError(
                    'dimension coordinate {!r} conflicts between '
                    'indexed and indexing objects:\n{}\nvs.\n{}'
                    .format(k, obj[k], coords[k]))
374692
def mapping(self, kdims=None, vdims=None, **kwargs):
        ""
        if config.future_deprecations:
            self.param.warning(
                "The mapping method is deprecated and should no "
                "longer be used. Use another one of the common "
                "formats instead, e.g. .dframe, .array or .columns.")

        length = len(self)
        if not kdims: kdims = self.kdims
        if kdims:
            keys = zip(*[self.dimension_values(dim.name)
                         for dim in self.kdims])
        else:
            keys = [()]*length

        if not vdims: vdims = self.vdims
        if vdims:
            values = zip(*[self.dimension_values(dim.name)
                           for dim in vdims])
        else:
            values = [()]*length
        return OrderedDict(zip(keys, values))
44985
def get_value_index(self, indices):
        
        size = self['size'] if self.get('size') else self['dimension']['size']
        ndims = len(size)
        mult = 1
        num = 0
        for idx, dim in enumerate(size):
            mult *= size[ndims - idx] if (idx > 0) else 1
            num += mult * indices[ndims - idx - 1]
        return num
1158965
def _get_coordinates_for_dataset_keys(self, dsids):
        """"""
        coordinates = {}
        for dsid in dsids:
            cids = self._get_coordinates_for_dataset_key(dsid)
            coordinates.setdefault(dsid, []).extend(cids)
        return coordinates
19461
def get_axes(x, y):
    
    n=int(np.sqrt(len(x)))
    if n**2 != len(x):
        raise ValueError("lenght of input array given as %s is not square of integer number!" % (len(x)))
    x_image = x.reshape(n,n)
    y_image = y.reshape(n,n)
    x_axes = x_image[0,:]
    y_axes = y_image[:,0]
    return x_axes, y_axes
1353433
def dimension_values(self, dim, expanded=True, flat=True):
        
        dim_idx = self.get_dimension_index(dim)
        if not expanded and dim_idx == 0:
            return np.array(range(self.data.shape[1]))
        elif not expanded and dim_idx == 1:
            return np.array(range(self.data.shape[0]))
        elif dim_idx in [0, 1]:
            values = np.mgrid[0:self.data.shape[1], 0:self.data.shape[0]][dim_idx]
            return values.flatten() if flat else values
        elif dim_idx == 2:
            arr = self.data.T
            return arr.flatten() if flat else arr
        else:
            return super(Raster, self).dimension_values(dim)
821968
def check_dimensions(self, dataset):
        '''
        
        '''
        required_ctx = TestCtx(BaseCheck.HIGH, 'All geophysical variables are point feature types')
        t = util.get_time_variable(dataset)

        # Exit prematurely
        if not t:
            required_ctx.assert_true(False, 'A dimension representing time is required for point feature types')
            return required_ctx.to_result()
        t_dims = dataset.variables[t].dimensions
        o = None or (t_dims and t_dims[0])

        message = '{} must be a valid timeseries feature type. It must have dimensions of ({}), and all coordinates must have dimensions of ({})'
        for variable in util.get_geophysical_variables(dataset):
            is_valid = util.is_point(dataset, variable)
            required_ctx.assert_true(
                is_valid,
                message.format(variable, o, o)
            )
        return required_ctx.to_result()
92056
def get_grid_spatial_dimensions(self, variable):
        """"""

        data = self.open_dataset(self.service).variables[variable.variable]
        dimensions = list(data.dimensions)
        return data.shape[dimensions.index(variable.x_dimension)], data.shape[dimensions.index(variable.y_dimension)]
593353
def interpolate(self, lat, lon, var):
        

        subset, dims = self.crop(lat, lon, var)

        if np.all([y in dims['lat'] for y in lat]) & \
                np.all([x in dims['lon'] for x in lon]):
                    yn = np.nonzero([y in lat for y in dims['lat']])[0]
                    xn = np.nonzero([x in lon for x in dims['lon']])[0]
                    output = {}
                    for v in subset:
                        # output[v] = subset[v][dn, zn, yn, xn]
                        # Seriously that this is the way to do it?!!??
                        output[v] = subset[v][:, xn][yn]
                    return output

        # The output coordinates shall be created only once.
        points_out = []
        for latn in lat:
            for lonn in lon:
                points_out.append([latn, lonn])
        points_out = np.array(points_out)

        output = {}
        for v in var:
            output[v] = ma.masked_all(
                    (lat.size, lon.size),
                    dtype=subset[v].dtype)

            # The valid data
            idx = np.nonzero(~ma.getmaskarray(subset[v]))

            if idx[0].size > 0:
                points = np.array([
                    dims['lat'][idx[0]], dims['lon'][idx[1]]]).T
                values = subset[v][idx]

                # Interpolate along the dimensions that have more than one
                #   position, otherwise it means that the output is exactly
                #   on that coordinate.
                ind = np.array(
                        [np.unique(points[:, i]).size > 1 for i in
                            range(points.shape[1])])
                assert ind.any()

                values_out = griddata(
                        np.atleast_1d(np.squeeze(points[:, ind])),
                        values,
                        np.atleast_1d(np.squeeze(points_out[:, ind]))
                        )

                # Remap the interpolated value back into a 4D array
                idx = np.isfinite(values_out)
                for [y, x], out in zip(points_out[idx], values_out[idx]):
                    output[v][y==lat, x==lon] = out

        return output