Source code for verde.projections

# Copyright (c) 2017 The Verde Developers.
# Distributed under the terms of the BSD 3-Clause License.
# SPDX-License-Identifier: BSD-3-Clause
#
# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
#
"""
Operations with projections for grids, regions, etc.
"""
import numpy as np

from .blockreduce import BlockReduce
from .chain import Chain
from .coordinates import check_region, get_region, grid_coordinates, shape_to_spacing
from .mask import convexhull_mask
from .neighbors import KNeighbors
from .scipygridder import Cubic, Linear
from .utils import grid_to_table


[docs] def project_region(region, projection): """ Calculate the bounding box of a region in projected coordinates. Parameters ---------- region : list = [W, E, S, N] The boundaries of a given region in Cartesian or geographic coordinates. projection : callable Should be a callable object (like a function) ``projection(easting, northing) -> (proj_easting, proj_northing)`` that takes in easting and northing coordinate arrays and returns projected northing and easting coordinate arrays. Returns ------- proj_region : list = [W, E, S, N] The bounding box of the projected region. Examples -------- >>> def projection(x, y): ... return (2*x, -1*y) >>> project_region((3, 5, -9, -4), projection) (6.0, 10.0, 4.0, 9.0) """ east, north = grid_coordinates(region, shape=(101, 101)) east, north = projection(east.ravel(), north.ravel()) return (east.min(), east.max(), north.min(), north.max())
[docs] def project_grid(grid, projection, method="linear", antialias=True, **kwargs): """ Apply the given map projection to a grid and re-sample it. Creates a new grid in the projected coordinates by interpolating the original values using the chosen *method* (linear by default). Before interpolation, apply a blocked mean operation (:class:`~verde.BlockReduce`) to avoid aliasing when the projected coordinates become oversampled in some regions (which would cause the interpolation to down-sample the original data). For example, applying a polar projection results in oversampled data close to the pole. Points that fall outside the convex hull of the original data will be masked (see :func:`~verde.convexhull_mask`) since they are not constrained by any data points. Any arguments that can be passed to the :meth:`~verde.base.BaseGridder.grid` method of Verde gridders can be passed to this function as well. Use this to set a region and spacing (or shape) for the projected grid. The region and spacing must be in *projected coordinates*. If no region is provided, the bounding box of the projected data will be used. If no spacing or shape is provided, the shape of the input *grid* will be used for the projected grid. By default, the ``data_names`` argument will be set to the name of the data variable of the input *grid* (if it has been set). .. note:: The interpolation methods are limited to what is available in Verde and there is only support for single 2D grids. For more sophisticated use cases, you might want to try `pyresample <https://github.com/pytroll/pyresample>`__ instead. Parameters ---------- grid : :class:`xarray.DataArray` A single 2D grid of values. The first dimension is assumed to be the northing/latitude dimension while the second is assumed to be the easting/longitude dimension. projection : callable Should be a callable object (like a function) ``projection(easting, northing) -> (proj_easting, proj_northing)`` that takes in easting and northing coordinate arrays and returns projected northing and easting coordinate arrays. method : string or Verde gridder If a string, will use it to create :class:`~verde.KNeighbors`, :class:`~verde.Linear`, or :class:`~verde.Cubic` (``"nearest"``, ``"linear"``, or ``"cubic"``). Otherwise, should be a gridder/estimator object, like :class:`~verde.Spline`. Default is ``"linear"``. antialias : bool If True, will run a :class:`~verde.BlockReduce` with a mean function to avoid aliasing when the projection results in oversampling of the data in some areas (for example, in polar projections). If False, will not run the blocked mean. Returns ------- projected_grid : :class:`xarray.DataArray` The projected grid, interpolated with the given parameters. """ if hasattr(grid, "data_vars"): raise ValueError( "Projecting xarray.Dataset is not currently supported. " "Please provide a DataArray instead." ) if len(grid.dims) != 2: raise ValueError( "Projecting grids with number of dimensions other than 2 is not " "currently supported (dimensions of the given DataArray: {}).".format( len(grid.dims) ) ) if isinstance(method, str): classes = dict( linear=Linear, nearest=KNeighbors, cubic=Cubic, ) if method not in classes: raise ValueError( f"Invalid interpolation method '{method}'. Must be one of {classes.keys()}." ) method = classes[method]() # Can be set to None for some data arrays depending on how they are created # so we can't just rely on the default value for getattr. name = getattr(grid, "name", None) if name is None: name = "scalars" data = grid_to_table(grid).dropna() coordinates = projection(data[grid.dims[1]].values, data[grid.dims[0]].values) data_region = get_region(coordinates) region = kwargs.pop("region", data_region) shape = kwargs.pop("shape", grid.shape) spacing = kwargs.pop("spacing", shape_to_spacing(region, shape)) check_region(region) steps = [] if antialias: steps.append( ("mean", BlockReduce(np.mean, spacing=spacing, region=data_region)) ) steps.append(("interpolator", method)) interpolator = Chain(steps) interpolator.fit(coordinates, data[name]) projected = interpolator.grid( region=region, spacing=spacing, data_names=kwargs.pop("data_names", [name]), **kwargs, ) projected = convexhull_mask(coordinates, grid=projected) return projected[name]