Merge remote-tracking branch 'origin/develop-3.X' into main

podpac/core/coordinates/array_coordinates1d.py

@@ -377,7 +377,46 @@ def _select(self, bounds, return_index, outer):
             gt = self.coordinates >= bounds[0]
             lt = self.coordinates <= bounds[1]
             b = gt & lt
-
+            b2 = gt | lt
+            if b2.sum() == b2.size and b.sum() == 0 and self.is_monotonic:
+                # bounds between data points
+                indlt = np.argwhere(lt).squeeze()
+                indgt = np.argwhere(gt).squeeze()
+                if self._is_descending:
+                    if indlt.size > 0:
+                        indlt = indlt[0]
+                    else:
+                        indlt = b.size - 1
+                    if indgt.size > 0:
+                        indgt = indgt[-1]
+                    else:
+                        indgt = 0
+                else:
+                    if indlt.size > 0:
+                        indlt = indlt[-1]
+                    else:
+                        indlt = 0
+                    if indgt.size > 0:
+                        indgt = indgt[0]
+                    else:
+                        indgt = b.size - 1
+
+                ind0 = min(indlt, indgt)
+                ind1 = max(indlt, indgt) + 1
+                b[ind0:ind1] = True
+                if b.sum() > 1:
+                    # These two coordinates are candidates, we need
+                    # to make sure that the bounds cross the edge between
+                    # the two points (selects both) or not (only selects)
+                    crds = self.coordinates[b]
+                    step = np.diff(self.coordinates[b])[0]
+                    edge = crds[0] + step / 2
+                    bounds_lt = bounds <= edge
+                    bounds_gt = bounds > edge
+                    keep_point = [np.any(bounds_lt), np.any(bounds_gt)]
+                    if self._is_descending:
+                        keep_point = keep_point[::-1]
+                    b[ind0:ind1] = keep_point
         elif self.is_monotonic:
             gt = np.where(self.coordinates >= bounds[0])[0]
             lt = np.where(self.coordinates <= bounds[1])[0]

podpac/core/coordinates/base_coordinates.py

@@ -69,6 +69,11 @@ def full_definition(self):
         """Coordinates definition, containing all properties. For internal use."""
         raise NotImplementedError
 
+    @property
+    def is_stacked(self):
+        """stacked or unstacked property"""
+        raise NotImplementedError
+
     @classmethod
     def from_definition(cls, d):
         """Get Coordinates from a coordinates definition."""
@@ -106,6 +111,10 @@ def issubset(self, other):
         """Report if these coordinates are a subset of other coordinates."""
         raise NotImplementedError
 
+    def horizontal_resolution(self, latitude, ellipsoid_tuple, coordinate_name, restype="nominal", units="meter"):
+        """Get horizontal resolution of coordiantes."""
+        raise NotImplementedError
+
     def __getitem__(self, index):
         raise NotImplementedError
 

podpac/core/coordinates/coordinates.py

@@ -19,6 +19,7 @@
 from six import string_types
 import pyproj
 import logging
+from scipy import spatial
 
 import podpac
 from podpac.core.settings import settings
@@ -478,9 +479,16 @@ def from_url(cls, url):
             r = 1
 
         # Extract bounding box information and translate to PODPAC coordinates
+        # Not, size does not get re-ordered, Height == Lat and width = lon -- ALWAYS
+        size = np.array([_get_param(params, "HEIGHT"), _get_param(params, "WIDTH")], int)
         start = bbox[:2][::r]
-        stop = bbox[2::][::r]
-        size = np.array([_get_param(params, "WIDTH"), _get_param(params, "HEIGHT")], int)[::r]
+        stop = bbox[2:][::r]
+
+        # The bbox gives the edges of the pixels, but our coordinates use the
+        # box centers -- so we have to shrink the start/stop portions
+        dx = (stop - start) / (size)  # This should take care of signs
+        start = start + dx / 2
+        stop = stop - dx / 2
 
         coords["coords"] = [
             {"name": "lat", "start": stop[0], "stop": start[0], "size": size[0]},
@@ -1501,12 +1509,84 @@ def issubset(self, other):
 
         return all(c.issubset(other) for c in self.values())
 
-    def is_stacked(self, dim):
-        if dim not in self.udims:
+    def is_stacked(self, dim):  # re-wrote to be able to iterate through c.dims
+        value = (dim in self.dims) + (dim in self.udims)
+        if value == 0:
             raise ValueError("Dimension {} is not in self.dims={}".format(dim, self.dims))
-        elif dim not in self.dims:
+        elif value == 1:  # one true, one false
             return True
-        return False
+        elif value == 2:  # both true
+            return False
+
+    def horizontal_resolution(self, units="meter", restype="nominal"):
+        """
+        Returns horizontal resolution of coordinate system.
+
+        Parameters
+        ----------
+        units : str
+            The desired unit the returned resolution should be in. Supports any unit supported by podpac.units (i.e. pint). Default is 'meter'.
+        restype : str
+            The kind of horizontal resolution that should be returned. Supported values are:
+            - "nominal" <-- Returns a number. Gives a 'nominal' resolution over the entire domain. This is wrong but fast.
+            - "summary" <-- Returns a tuple (mean, standard deviation). Gives the exact mean and standard deviation for unstacked coordinates, some error for stacked coordinates
+            - "full" <-- Returns a 1 or 2-D array. Gives exact grid differences if unstacked coordinates or distance matrix if stacked coordinates
+
+        Returns
+        -------
+        OrderedDict
+            A dictionary with:
+            keys : str
+                dimension names
+            values
+                resolution (format determined by 'type' parameter)
+
+        Raises
+        ------
+        ValueError
+            If the 'restype' is not one of the supported resolution types
+
+
+        """
+        # This function handles mainly edge case sanitation.
+        # It calls StackedCoordinates and Coordinates1d 'horizontal_resolution' methods to get the actual values.
+
+        if "lat" not in self.udims:  # require latitude
+            raise ValueError("Latitude required for horizontal resolution.")
+
+        # ellipsoid tuple to pass to geodesic
+        ellipsoid_tuple = (
+            self.CRS.ellipsoid.semi_major_metre / 1000,
+            self.CRS.ellipsoid.semi_minor_metre / 1000,
+            1 / self.CRS.ellipsoid.inverse_flattening,
+        )
+
+        # main execution loop
+        resolutions = OrderedDict()  # To return
+        for name, dim in self.items():
+            if dim.is_stacked:
+                if "lat" in dim.dims and "lon" in dim.dims:
+                    resolutions[name] = dim.horizontal_resolution(
+                        None, ellipsoid_tuple, self.CRS.coordinate_system.name, restype, units
+                    )
+                elif "lat" in dim.dims:
+                    # Calling self['lat'] forces UniformCoordinates1d, even if stacked
+                    resolutions["lat"] = self["lat"].horizontal_resolution(
+                        self["lat"], ellipsoid_tuple, self.CRS.coordinate_system.name, restype, units
+                    )
+                elif "lon" in dim.dims:
+                    # Calling self['lon'] forces UniformCoordinates1d, even if stacked
+                    resolutions["lon"] = self["lon"].dim.horizontal_resolution(
+                        self["lat"], ellipsoid_tuple, self.CRS.coordinate_system.name, restype, units
+                    )
+            elif (
+                name == "lat" or name == "lon"
+            ):  # need to do this inside of loop in case of stacked [[alt,time]] but unstacked [lat, lon]
+                resolutions[name] = dim.horizontal_resolution(
+                    self["lat"], ellipsoid_tuple, self.CRS.coordinate_system.name, restype, units
+                )
+
+        return resolutions
 
     # ------------------------------------------------------------------------------------------------------------------
     # Operators/Magic Methods

podpac/core/coordinates/coordinates1d.py

@@ -10,10 +10,12 @@
 import numpy as np
 import traitlets as tl
 
+import podpac
 from podpac.core.utils import ArrayTrait, TupleTrait
 from podpac.core.coordinates.utils import make_coord_value, make_coord_delta, make_coord_delta_array
 from podpac.core.coordinates.utils import add_coord, divide_delta, lower_precision_time_bounds
 from podpac.core.coordinates.utils import Dimension
+from podpac.core.coordinates.utils import calculate_distance
 from podpac.core.coordinates.base_coordinates import BaseCoordinates
 
 
@@ -187,6 +189,10 @@ def _get_definition(self, full=True):
     def _full_properties(self):
         return {"name": self.name}
 
+    @property
+    def is_stacked(self):
+        return False
+
     # ------------------------------------------------------------------------------------------------------------------
     # Methods
     # ------------------------------------------------------------------------------------------------------------------
@@ -423,3 +429,118 @@ def issubset(self, other):
                 other_coordinates = other_coordinates.astype(my_coordinates.dtype)
 
         return set(my_coordinates).issubset(other_coordinates)
+
+    def horizontal_resolution(self, latitude, ellipsoid_tuple, coordinate_name, restype="nominal", units="meter"):
+        """Return the horizontal resolution of a Uniform 1D Coordinate
+
+        Parameters
+        ----------
+        latitude: Coordinates1D
+            The latitude coordiantes of the current coordinate system, required for both lat and lon resolution.
+        ellipsoid_tuple: tuple
+            a tuple containing ellipsoid information from the the original coordinates to pass to geopy
+        coordinate_name: str
+            "cartesian" or "ellipsoidal", to tell calculate_distance what kind of calculation to do
+        restype: str
+            The kind of horizontal resolution that should be returned. Supported values are:
+            - "nominal" <-- returns average resolution based upon bounds and number of grid squares
+            - "summary" <-- Gives the exact mean and standard deviation of grid square resolutions
+            - "full" <-- Gives exact grid differences.
+
+        Returns
+        -------
+        float * (podpac.unit)
+            if restype == "nominal", return average distance in desired units
+        tuple * (podpac.unit)
+            if restype == "summary", return average distance and standard deviation in desired units
+        np.ndarray * (podpac.unit)
+            if restype == "full", give full resolution. 1D array for latitude, MxN matrix for longitude.
+        """
+
+        def nominal_unstacked_resolution():
+            """Return resolution for unstacked coordiantes using the bounds
+
+            Returns
+            --------
+            The average distance between each grid square for this dimension
+            """
+            if self.name == "lat":
+                return (
+                    calculate_distance(
+                        (self.bounds[0], 0), (self.bounds[1], 0), ellipsoid_tuple, coordinate_name, units
+                    ).magnitude
+                    / (self.size - 1)
+                ) * podpac.units(units)
+            elif self.name == "lon":
+                median_lat = ((latitude.bounds[1] - latitude.bounds[0]) / 2) + latitude.bounds[0]
+                return (
+                    calculate_distance(
+                        (median_lat, self.bounds[0]),
+                        (median_lat, self.bounds[1]),
+                        ellipsoid_tuple,
+                        coordinate_name,
+                        units,
+                    ).magnitude
+                    / (self.size - 1)
+                ) * podpac.units(units)
+            else:
+                return ValueError("Unknown dim: {}".format(self.name))
+
+        def summary_unstacked_resolution():
+            """Return summary resolution for the dimension.
+
+            Returns
+            -------
+            tuple
+                the average distance between grid squares
+                the standard deviation of those distances
+            """
+            if self.name == "lat" or self.name == "lon":
+                full_res = full_unstacked_resolution().magnitude
+                return (np.average(full_res) * podpac.units(units), np.std(full_res) * podpac.units(units))
+            else:
+                return ValueError("Unknown dim: {}".format(self.name))
+
+        def full_unstacked_resolution():
+            """Calculate full resolution of unstacked dimension
+
+            Returns
+            -------
+            A matrix of distances
+            """
+            if self.name == "lat":
+                top_bounds = np.stack(
+                    [latitude.coordinates[1:], np.full((latitude.coordinates[1:]).shape[0], 0)], axis=1
+                )  # use exact lat values
+                bottom_bounds = np.stack(
+                    [latitude.coordinates[:-1], np.full((latitude.coordinates[:-1]).shape[0], 0)], axis=1
+                )  # use exact lat values
+                return calculate_distance(top_bounds, bottom_bounds, ellipsoid_tuple, coordinate_name, units)
+            elif self.name == "lon":
+                M = latitude.coordinates.size
+                N = self.size
+                diff = np.zeros((M, N - 1))
+                for i in range(M):
+                    lat_value = latitude.coordinates[i]
+                    lon_values = self.coordinates
+                    top_bounds = np.stack(
+                        [np.full((lon_values[1:]).shape[0], lat_value), lon_values[1:]], axis=1
+                    )  # use exact lat values
+                    bottom_bounds = np.stack(
+                        [np.full((lon_values[:-1]).shape[0], lat_value), lon_values[:-1]], axis=1
+                    )  # use exact lat values
+                    diff[i] = calculate_distance(
+                        top_bounds, bottom_bounds, ellipsoid_tuple, coordinate_name, units
+                    ).magnitude
+                return diff * podpac.units(units)
+            else:
+                raise ValueError("Unknown dim: {}".format(self.name))
+
+        if restype == "nominal":
+            return nominal_unstacked_resolution()
+        elif restype == "summary":
+            return summary_unstacked_resolution()
+        elif restype == "full":
+            return full_unstacked_resolution()
+        else:
+            raise ValueError("Invalid value for type: {}".format(restype))