diff --git a/podpac/core/coordinates/array_coordinates1d.py b/podpac/core/coordinates/array_coordinates1d.py
index 2d4c2f4a..2bc03897 100644
--- a/podpac/core/coordinates/array_coordinates1d.py
+++ b/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]
diff --git a/podpac/core/coordinates/base_coordinates.py b/podpac/core/coordinates/base_coordinates.py
index 47976403..2b14e4b9 100644
--- a/podpac/core/coordinates/base_coordinates.py
+++ b/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
 
diff --git a/podpac/core/coordinates/coordinates.py b/podpac/core/coordinates/coordinates.py
index 75653fd9..b5885099 100644
--- a/podpac/core/coordinates/coordinates.py
+++ b/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
diff --git a/podpac/core/coordinates/coordinates1d.py b/podpac/core/coordinates/coordinates1d.py
index e855255e..38972497 100644
--- a/podpac/core/coordinates/coordinates1d.py
+++ b/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))
diff --git a/podpac/core/coordinates/stacked_coordinates.py b/podpac/core/coordinates/stacked_coordinates.py
index 9ffd45b4..3d7a1861 100644
--- a/podpac/core/coordinates/stacked_coordinates.py
+++ b/podpac/core/coordinates/stacked_coordinates.py
@@ -9,12 +9,15 @@
 import traitlets as tl
 from six import string_types
 import lazy_import
+from scipy import spatial
 
+import podpac
 from podpac.core.coordinates.base_coordinates import BaseCoordinates
 from podpac.core.coordinates.coordinates1d import Coordinates1d
 from podpac.core.coordinates.array_coordinates1d import ArrayCoordinates1d
 from podpac.core.coordinates.uniform_coordinates1d import UniformCoordinates1d
 from podpac.core.coordinates.utils import make_coord_value
+from podpac.core.coordinates.utils import calculate_distance
 
 
 class StackedCoordinates(BaseCoordinates):
@@ -358,6 +361,10 @@ def full_definition(self):
 
         return [c.full_definition for c in self._coords]
 
+    @property
+    def is_stacked(self):
+        return True
+
     # -----------------------------------------------------------------------------------------------------------------
     # Methods
     # -----------------------------------------------------------------------------------------------------------------
@@ -700,3 +707,94 @@ def is_affine(self):
             return False
 
         return True
+
+    def horizontal_resolution(self, latitude, ellipsoid_tuple, coordinate_name, restype="nominal", units="meter"):
+        """Return the horizontal resolution of a Uniform 1D Coordinate
+
+        Parameters
+        ----------
+        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" <-- Gives average nearest distance of all points, with some error
+            - "summary" <-- Gives the mean and standard deviation of nearest distance betweem points, with some error
+            - "full" <-- Gives exact distance matrix
+        units: str
+            desired unit to return
+
+        Returns
+        -------
+        float * (podpac.unit)
+            If restype == "nominal", return the average nearest distance with some error
+        tuple * (podpac.unit)
+            If restype == "summary", return average and std.dev of nearest distances, with some error
+        np.ndarray * (podpac.unit)
+            if restype == "full", return exact distance matrix
+        ValueError
+            if unknown restype
+
+        """
+        order = tuple([self.dims.index(d) for d in ["lat", "lon"]])
+
+        def nominal_stacked_resolution():
+            """Use a KDTree to return approximate stacked resolution with some loss of accuracy.
+
+            Returns
+            -------
+            The average min distance of every point
+
+            """
+            tree = spatial.KDTree(self.coordinates[:, order] + [0, 180.0], boxsize=[0.0, 360.0000000000001])
+            return np.average(
+                calculate_distance(
+                    tree.data - [0, 180.0],
+                    tree.data[tree.query(tree.data, k=2)[1][:, 1]] - [0, 180.0],
+                    ellipsoid_tuple,
+                    coordinate_name,
+                    units,
+                )
+            )
+
+        def summary_stacked_resolution():
+            """Return the approximate mean resolution and std.deviation using a KDTree
+
+            Returns
+            -------
+            tuple
+                Average min distance of every point and standard deviation of those min distances
+            """
+            tree = spatial.KDTree(self.coordinates[:, order] + [0, 180.0], boxsize=[0.0, 360.0000000000001])
+            distances = calculate_distance(
+                tree.data - [0, 180.0],
+                tree.data[tree.query(tree.data, k=2)[1][:, 1]] - [0, 180.0],
+                ellipsoid_tuple,
+                coordinate_name,
+                units,
+            )
+            return (np.average(distances), np.std(distances))
+
+        def full_stacked_resolution():
+            """Returns the exact distance between every point using brute force
+
+            Returns
+            -------
+            distance matrix of size (NxN), where N is the number of points in the dimension
+            """
+            distance_matrix = np.zeros((len(self.coordinates), len(self.coordinates)))
+            for i in range(len(self.coordinates)):
+                distance_matrix[i, :] = calculate_distance(
+                    self.coordinates[i, order], self.coordinates[:, order], ellipsoid_tuple, coordinate_name, units
+                ).magnitude
+            return distance_matrix * podpac.units(units)
+
+        if restype == "nominal":
+            return nominal_stacked_resolution()
+        elif restype == "summary":
+            return summary_stacked_resolution()
+        elif restype == "full":
+            return full_stacked_resolution()
+        else:
+            raise ValueError("Invalid value for type: {}".format(restype))
diff --git a/podpac/core/coordinates/test/test_array_coordinates1d.py b/podpac/core/coordinates/test/test_array_coordinates1d.py
index 211e62fa..804393bf 100644
--- a/podpac/core/coordinates/test/test_array_coordinates1d.py
+++ b/podpac/core/coordinates/test/test_array_coordinates1d.py
@@ -836,6 +836,47 @@ def test_select(self):
         assert_equal(s.coordinates, [])
         assert_equal(c.coordinates[I], [])
 
+    def test_select_one_between_coords(self):
+        # Ascending
+        c = ArrayCoordinates1d([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
+        c1, inds1 = c.select([5.6, 6.1], return_index=True, outer=False)
+        assert np.argwhere(inds1).squeeze() == 6
+
+        c2, inds2 = c.select([5.6, 5.6], return_index=True, outer=False)
+        assert np.argwhere(inds2).squeeze() == 6
+
+        c3, inds3 = c.select([5.4, 5.4], return_index=True, outer=False)
+        assert np.argwhere(inds3).squeeze() == 5
+
+        c3b, inds3b = c.select([5.4, 5.6], return_index=True, outer=False)
+        assert np.all(np.argwhere(inds3b).squeeze() == [5, 6])
+
+        c4, inds4 = c.select([9.1, 9.1], return_index=True, outer=False)
+        assert inds4 == []
+
+        c5, inds5 = c.select([-0.1, -0.1], return_index=True, outer=False)
+        assert inds5 == []
+
+        # Decending
+        c = ArrayCoordinates1d([0, 1, 2, 3, 4, 5, 6, 7, 8, 9][::-1])
+        c1, inds1 = c.select([5.6, 6.1], return_index=True, outer=False)
+        assert np.argwhere(inds1).squeeze() == 3
+
+        c2, inds2 = c.select([5.6, 5.6], return_index=True, outer=False)
+        assert np.argwhere(inds2).squeeze() == 3
+
+        c3, inds3 = c.select([5.4, 5.4], return_index=True, outer=False)
+        assert np.argwhere(inds3).squeeze() == 4
+
+        c3b, inds3b = c.select([5.4, 5.6], return_index=True, outer=False)
+        assert np.all(np.argwhere(inds3b).squeeze() == [3, 4])
+
+        c4, inds4 = c.select([9.1, 9.1], return_index=True, outer=False)
+        assert inds4 == []
+
+        c5, inds5 = c.select([-0.1, -0.1], return_index=True, outer=False)
+        assert inds5 == []
+
     def test_select_outer_ascending(self):
         c = ArrayCoordinates1d([10.0, 20.0, 40.0, 50.0, 60.0, 90.0])
 
diff --git a/podpac/core/coordinates/test/test_coordinates.py b/podpac/core/coordinates/test/test_coordinates.py
index 97db6e36..b4dd1d7a 100644
--- a/podpac/core/coordinates/test/test_coordinates.py
+++ b/podpac/core/coordinates/test/test_coordinates.py
@@ -8,6 +8,8 @@
 from numpy.testing import assert_equal, assert_array_equal
 import xarray as xr
 import pyproj
+from collections import OrderedDict
+import pint
 
 import podpac
 from podpac.core.coordinates.coordinates1d import Coordinates1d
@@ -378,57 +380,63 @@ def test_from_url(self):
 
         url = (
             r"http://testwms/?map=map&&service=WMS&request=GetMap&layers=layer&styles=&format=image%2Fpng"
-            r"&transparent=true&version={version}&transparency=true&width=256&height=256&srs=EPSG%3A{epsg}"
+            r"&transparent=true&version={version}&transparency=true&width=128&height=256&srs=EPSG%3A{epsg}"
             r"&bbox={},{},{},{}&time={}"
         )
 
         # version 1.1.1
         c = Coordinates.from_url(
             url.format(
-                min(crds2.bounds["lon"]),
-                min(crds2.bounds["lat"]),
-                max(crds2.bounds["lon"]),
-                max(crds2.bounds["lat"]),
+                min(crds2.bounds["lon"]) - np.abs(crds2["lon"].step / 127 / 2),
+                min(crds2.bounds["lat"]) - np.abs(crds2["lat"].step / 255 / 2),
+                max(crds2.bounds["lon"]) + np.abs(crds2["lon"].step / 127 / 2),
+                max(crds2.bounds["lat"]) + np.abs(crds2["lat"].step / 255 / 2),
                 crds2.bounds["time"][0],
                 version="1.1.1",
                 epsg="3857",
             )
         )
         assert c.bounds == crds2.bounds
+        assert crds2.bounds["lon"][0] - c["lon"].step / 2 == c.geotransform[0]
+        assert crds2.bounds["lat"][1] - c["lat"].step / 2 == c.geotransform[3]
 
         c = Coordinates.from_url(
             url.format(
-                min(crds.bounds["lon"]),
-                min(crds.bounds["lat"]),
-                max(crds.bounds["lon"]),
-                max(crds.bounds["lat"]),
+                min(crds.bounds["lon"]) - np.abs(crds["lon"].step / 127 / 2),
+                min(crds.bounds["lat"]) - np.abs(crds["lat"].step / 255 / 2),
+                max(crds.bounds["lon"]) + np.abs(crds["lon"].step / 127 / 2),
+                max(crds.bounds["lat"]) + np.abs(crds["lat"].step / 255 / 2),
                 crds.bounds["time"][0],
                 version="1.1.1",
                 epsg="4326",
             )
         )
         assert c.bounds == crds.bounds
+        assert crds.bounds["lon"][0] - c["lon"].step / 2 == c.geotransform[0]
+        assert crds.bounds["lat"][1] - c["lat"].step / 2 == c.geotransform[3]
 
         # version 1.3
         c = Coordinates.from_url(
             url.format(
-                min(crds2.bounds["lon"]),
-                min(crds2.bounds["lat"]),
-                max(crds2.bounds["lon"]),
-                max(crds2.bounds["lat"]),
+                min(crds2.bounds["lon"]) - np.abs(crds2["lon"].step / 127 / 2),
+                min(crds2.bounds["lat"]) - np.abs(crds2["lat"].step / 255 / 2),
+                max(crds2.bounds["lon"]) + np.abs(crds2["lon"].step / 127 / 2),
+                max(crds2.bounds["lat"]) + np.abs(crds2["lat"].step / 255 / 2),
                 crds2.bounds["time"][0],
                 version="1.3",
                 epsg="3857",
             )
         )
         assert c.bounds == crds2.bounds
+        assert crds2.bounds["lon"][0] - c["lon"].step / 2 == c.geotransform[0]
+        assert crds2.bounds["lat"][1] - c["lat"].step / 2 == c.geotransform[3]
 
         c = Coordinates.from_url(
             url.format(
-                min(crds.bounds["lat"]),
-                min(crds.bounds["lon"]),
-                max(crds.bounds["lat"]),
-                max(crds.bounds["lon"]),
+                min(crds.bounds["lat"]) - np.abs(crds["lat"].step / 255 / 2),
+                min(crds.bounds["lon"]) - np.abs(crds["lon"].step / 127 / 2),
+                max(crds.bounds["lat"]) + np.abs(crds["lat"].step / 255 / 2),
+                max(crds.bounds["lon"]) + np.abs(crds["lon"].step / 127 / 2),
                 crds.bounds["time"][0],
                 version="1.3",
                 epsg="4326",
@@ -436,6 +444,8 @@ def test_from_url(self):
         )
 
         assert c.bounds == crds.bounds
+        assert crds.bounds["lon"][0] - c["lon"].step / 2 == c.geotransform[0]
+        assert crds.bounds["lat"][1] - c["lat"].step / 2 == c.geotransform[3]
 
         # WCS version
         crds = Coordinates([[41, 40], [-71, -70], "2018-05-19"], dims=["lat", "lon", "time"])
@@ -443,23 +453,24 @@ def test_from_url(self):
 
         url = (
             r"http://testwms/?map=map&&service=WCS&request=GetMap&layers=layer&styles=&format=image%2Fpng"
-            r"&transparent=true&version={version}&transparency=true&width=256&height=256&srs=EPSG%3A{epsg}"
+            r"&transparent=true&version={version}&transparency=true&width=128&height=256&srs=EPSG%3A{epsg}"
             r"&bbox={},{},{},{}&time={}"
         )
 
         c = Coordinates.from_url(
             url.format(
-                min(crds2.bounds["lon"]),
-                min(crds2.bounds["lat"]),
-                max(crds2.bounds["lon"]),
-                max(crds2.bounds["lat"]),
+                min(crds2.bounds["lon"]) - np.abs(crds2["lon"].step / 127 / 2),
+                min(crds2.bounds["lat"]) - np.abs(crds2["lat"].step / 255 / 2),
+                max(crds2.bounds["lon"]) + np.abs(crds2["lon"].step / 127 / 2),
+                max(crds2.bounds["lat"]) + np.abs(crds2["lat"].step / 255 / 2),
                 crds2.bounds["time"][0],
                 version="1.1",
                 epsg="3857",
             )
         )
         assert c.bounds == crds2.bounds
-
+        assert crds2.bounds["lon"][0] - c["lon"].step / 2 == c.geotransform[0]
+        assert crds2.bounds["lat"][1] - c["lat"].step / 2 == c.geotransform[3]
         # Based on all the documentation I've read, this should be correct, but
         # based on the server's I've checked, this does not seem correct
         # c = Coordinates.from_url(
@@ -476,10 +487,10 @@ def test_from_url(self):
 
         c = Coordinates.from_url(
             url.format(
-                min(crds.bounds["lon"]),
-                min(crds.bounds["lat"]),
-                max(crds.bounds["lon"]),
-                max(crds.bounds["lat"]),
+                min(crds.bounds["lon"]) - np.abs(crds["lon"].step / 127 / 2),
+                min(crds.bounds["lat"]) - np.abs(crds["lat"].step / 255 / 2),
+                max(crds.bounds["lon"]) + np.abs(crds["lon"].step / 127 / 2),
+                max(crds.bounds["lat"]) + np.abs(crds["lat"].step / 255 / 2),
                 crds.bounds["time"][0],
                 version="1.1",
                 epsg="4326",
@@ -487,6 +498,8 @@ def test_from_url(self):
         )
 
         assert c.bounds == crds.bounds
+        assert crds.bounds["lon"][0] - c["lon"].step / 2 == c.geotransform[0]
+        assert crds.bounds["lat"][1] - c["lat"].step / 2 == c.geotransform[3]
 
     def test_from_xarray(self):
         lat = [0, 1, 2]
@@ -1445,7 +1458,7 @@ def test_intersect(self):
         assert ct2.size == 3
 
         ct2 = ct.intersect(cti, outer=False)
-        assert ct2.size == 0  # Is this behavior desired?
+        assert ct2.size == 1  #
 
     def test_intersect_dims(self):
         lat = ArrayCoordinates1d([0, 1, 2, 3, 4, 5], name="lat")
@@ -2202,3 +2215,53 @@ def test_affine_to_uniform(self):
 
         # unstacked uniform -> unstacked uniform
         assert c3.simplify() == c3
+
+
+class TestResolutions(object):
+    def test_horizontal_resolution(self):
+        """Test edge cases of resolutions, and that Resolutions are returned correctly in an OrderedDict.
+        StackedCoordinates and Coordiantes1d handle the resolution calculations, so correctness of the resolutions are tested there.
+        """
+
+        # Dimensions
+        lat = podpac.clinspace(-80, 80, 5)
+        lon = podpac.clinspace(-180, 180, 5)
+        time = ["2018-01-01", "2018-01-02", "2018-01-03", "2018-01-04", "2018-01-05"]
+
+        # Invalid stacked dims check:
+        c = Coordinates([[lon, time]], dims=["lon_time"])
+        with pytest.raises(ValueError):
+            c.horizontal_resolution()
+
+        # Require latitude
+        c = Coordinates([lon], dims=["lon"])
+        with pytest.raises(ValueError):
+            c.horizontal_resolution()
+
+        # Valid dims check:
+        c = Coordinates([[lat, lon]], dims=["lat_lon"])
+        c.horizontal_resolution()
+
+        c = Coordinates([[lon, lat]], dims=["lon_lat"])
+        c.horizontal_resolution()
+
+        # Corect name for restype:
+        with pytest.raises(ValueError):
+            c.horizontal_resolution(restype="whateverIwant")
+
+        # Unstacked
+        c = Coordinates([lat, lon], dims=["lat", "lon"])
+        c.horizontal_resolution()
+
+        c = Coordinates([lat, lon], dims=["lat", "lon"])
+        c.horizontal_resolution(restype="summary")
+
+        # Mixed stacked, unstacked, but still valid
+        # Stacked and Unstacked valid Check:
+        c = Coordinates([[lat, time], lon], dims=["lat_time", "lon"])
+        c2 = Coordinates([lat, lon], dims=["lat", "lon"])
+        np.testing.assert_array_equal(c.horizontal_resolution(), c2.horizontal_resolution())
+        # Lat only
+        c = Coordinates([[lat, time]], dims=["lat_time"])
+        c2 = Coordinates([lat], dims=["lat"])
+        np.testing.assert_array_equal(c.horizontal_resolution(), c2.horizontal_resolution())
diff --git a/podpac/core/coordinates/test/test_coordinates1d.py b/podpac/core/coordinates/test/test_coordinates1d.py
index 35d1422f..b7f7d3d5 100644
--- a/podpac/core/coordinates/test/test_coordinates1d.py
+++ b/podpac/core/coordinates/test/test_coordinates1d.py
@@ -1,6 +1,9 @@
 import pytest
+import numpy as np
 
+import podpac
 from podpac.core.coordinates.coordinates1d import Coordinates1d
+from podpac import clinspace
 
 
 class TestCoordinates1d(object):
@@ -82,3 +85,59 @@ def test_common_api(self):
             c.issubset(c)
         except NotImplementedError:
             pass
+
+    def test_horizontal_resolution(self):
+        """Test horizontal resolution implentation for Coordinates1d. Edge cases are handled in Coordinates.py"""
+        # Latitude
+        lat = clinspace(-80, 80, 5)
+        lat.name = "lat"  # normally assigned when creating Coords object
+        assert type(lat) == podpac.core.coordinates.uniform_coordinates1d.UniformCoordinates1d
+
+        # Longitude
+        lon = podpac.clinspace(-180, 180, 5)
+        lon.name = "lon"
+        assert type(lon) == podpac.core.coordinates.uniform_coordinates1d.UniformCoordinates1d
+
+        # Sample Ellipsoid Tuple
+        ell_tuple = (6378.137, 6356.752314245179, 0.0033528106647474805)
+
+        # Sample Coordinate name:
+        coord_name = "ellipsoidal"
+
+        # Resolution: nominal
+        assert lat.horizontal_resolution(lat, ell_tuple, coord_name) == 4442569.935968436 * podpac.units("meter")
+        assert lon.horizontal_resolution(lat, ell_tuple, coord_name) == 0.0 * podpac.units("meter")
+
+        # Resolution: summary
+        assert lat.horizontal_resolution(lat, ell_tuple, coord_name, restype="summary") == (
+            4442569.935968436 * podpac.units("meter"),
+            13040.905617921147 * podpac.units("meter"),
+        )
+        assert lon.horizontal_resolution(lat, ell_tuple, coord_name, restype="summary") == (
+            5558704.3695234 * podpac.units("meter"),
+            3399219.0171971265 * podpac.units("meter"),
+        )
+
+        # Resolution: full
+        lat_answer = [4455610.84158636, 4429529.03035052, 4429529.03035052, 4455610.84158636]
+
+        lon_answer = [
+            [1575399.99090356, 1575399.99090356, 1575399.99090356, 1575399.99090356],
+            [7311983.84720763, 7311983.84720763, 7311983.84720763, 7311983.84720763],
+            [10018754.17139462, 10018754.17139462, 10018754.17139462, 10018754.17139462],
+            [7311983.84720763, 7311983.84720763, 7311983.84720763, 7311983.84720763],
+            [1575399.99090356, 1575399.99090356, 1575399.99090356, 1575399.99090356],
+        ]
+
+        np.testing.assert_array_almost_equal(
+            lat.horizontal_resolution(lat, ell_tuple, coord_name, restype="full").magnitude, lat_answer
+        )
+        np.testing.assert_array_almost_equal(
+            lon.horizontal_resolution(lat, ell_tuple, coord_name, restype="full").magnitude, lon_answer
+        )
+
+        # Different Units
+        np.testing.assert_almost_equal(
+            lat.horizontal_resolution(lat, ell_tuple, coord_name).to(podpac.units("feet")).magnitude,
+            lat.horizontal_resolution(lat, ell_tuple, coord_name, units="feet").magnitude,
+        )
diff --git a/podpac/core/coordinates/test/test_stacked_coordinates.py b/podpac/core/coordinates/test/test_stacked_coordinates.py
index d9460a05..c2b78440 100644
--- a/podpac/core/coordinates/test/test_stacked_coordinates.py
+++ b/podpac/core/coordinates/test/test_stacked_coordinates.py
@@ -981,3 +981,69 @@ def test_reshape(self):
 
         assert c.reshape((4, 3)) == StackedCoordinates([lat.reshape((4, 3)), lon.reshape((4, 3))])
         assert c.flatten().reshape((3, 4)) == c
+
+    def test_horizontal_resolution(self):
+        """Test Horizontal Resolution of Stacked Coordinates. Edge cases are handled in Coordinates.py"""
+        lat = podpac.clinspace(-80, 80, 5)
+        lat.name = "lat"  # normally assigned when creating Coords object
+        lon = podpac.clinspace(-180, 180, 5)
+        lon.name = "lon"
+        c = StackedCoordinates([lat, lon])
+
+        # Sample Ellipsoid Tuple
+        ell_tuple = (6378.137, 6356.752314245179, 0.0033528106647474805)
+
+        # Sample Coordinate name:
+        coord_name = "ellipsoidal"
+
+        # Nominal resolution:
+        np.testing.assert_almost_equal(
+            c.horizontal_resolution(None, ell_tuple, coord_name, restype="nominal").magnitude, 7397047.845631437
+        )
+
+        # Summary resolution
+        np.testing.assert_almost_equal(
+            c.horizontal_resolution(None, ell_tuple, coord_name, restype="summary")[0].magnitude, 7397047.845631437
+        )
+        np.testing.assert_almost_equal(
+            c.horizontal_resolution(None, ell_tuple, coord_name, restype="summary")[1].magnitude, 2134971.4571846593
+        )
+
+        # Full resolution
+        distance_matrix = [
+            [0.0, 5653850.95046188, 11118791.58668857, 14351078.11393555, 17770279.74387375],
+            [5653850.95046188, 0.0, 10011843.18838578, 20003931.45862544, 14351078.11393555],
+            [11118791.58668857, 10011843.18838578, 0.0, 10011843.18838578, 11118791.58668857],
+            [14351078.11393555, 20003931.45862544, 10011843.18838578, 0.0, 5653850.95046188],
+            [17770279.74387375, 14351078.11393555, 11118791.58668857, 5653850.95046188, 0.0],
+        ]
+
+        np.testing.assert_array_almost_equal(
+            c.horizontal_resolution(None, ell_tuple, coord_name, restype="full"), distance_matrix
+        )
+
+        # Test different order of lat/lon still works
+        c2 = StackedCoordinates([lon, lat])
+        np.testing.assert_equal(
+            c.horizontal_resolution(None, ell_tuple, "lat").magnitude,
+            c2.horizontal_resolution(None, ell_tuple, "lat").magnitude,
+        )
+
+        # Test multiple stacked coordinates
+        alt = podpac.clinspace(0, 1, 5, "alt")
+
+        c2 = StackedCoordinates([lon, alt, lat])
+        np.testing.assert_equal(
+            c.horizontal_resolution(None, ell_tuple, "lat").magnitude,
+            c2.horizontal_resolution(None, ell_tuple, "lat").magnitude,
+        )
+        c2 = StackedCoordinates([alt, lon, lat])
+        np.testing.assert_equal(
+            c.horizontal_resolution(None, ell_tuple, "lat").magnitude,
+            c2.horizontal_resolution(None, ell_tuple, "lat").magnitude,
+        )
+        c2 = StackedCoordinates([lat, alt, lon])
+        np.testing.assert_equal(
+            c.horizontal_resolution(None, ell_tuple, "lat").magnitude,
+            c2.horizontal_resolution(None, ell_tuple, "lat").magnitude,
+        )
diff --git a/podpac/core/coordinates/test/test_uniform_coordinates1d.py b/podpac/core/coordinates/test/test_uniform_coordinates1d.py
index 51f1924a..a25b40a0 100644
--- a/podpac/core/coordinates/test/test_uniform_coordinates1d.py
+++ b/podpac/core/coordinates/test/test_uniform_coordinates1d.py
@@ -944,13 +944,13 @@ def test_select_ascending(self):
 
         # between coordinates
         s = c.select([52, 55])
-        assert isinstance(s, ArrayCoordinates1d)
-        assert_equal(s.coordinates, [])
+        assert isinstance(s, UniformCoordinates1d)
+        assert_equal(s.coordinates, c[3:5].coordinates)
 
         s, I = c.select([52, 55], return_index=True)
-        assert isinstance(s, ArrayCoordinates1d)
-        assert_equal(s.coordinates, [])
-        assert_equal(c.coordinates[I], [])
+        assert isinstance(s, UniformCoordinates1d)
+        assert_equal(s.coordinates, c[3:5].coordinates)
+        assert_equal(c.coordinates[I], c[3:5].coordinates)
 
         # backwards bounds
         s = c.select([70, 30])
@@ -1015,13 +1015,13 @@ def test_select_descending(self):
 
         # between coordinates
         s = c.select([52, 55])
-        assert isinstance(s, ArrayCoordinates1d)
-        assert_equal(s.coordinates, [])
+        assert isinstance(s, UniformCoordinates1d)
+        assert_equal(s.coordinates, [60, 50])
 
         s, I = c.select([52, 55], return_index=True)
-        assert isinstance(s, ArrayCoordinates1d)
-        assert_equal(s.coordinates, [])
-        assert_equal(c.coordinates[I], [])
+        assert isinstance(s, UniformCoordinates1d)
+        assert_equal(s.coordinates, [60, 50])
+        assert_equal(c.coordinates[I], [60, 50])
 
         # backwards bounds
         s = c.select([70, 30])
@@ -1116,6 +1116,26 @@ def test_select_time_variable_precision(self):
         assert s1.size == 0
         assert s2.size == 1
 
+    def test_select_one_floating_point_error(self):
+        c = UniformCoordinates1d(0, 9, 1)  # FIX THE PROBLEM HERE!!
+        c1, inds1 = c.select([5.6, 6.1], return_index=True, outer=False)
+        assert inds1 == slice(6, 7)
+
+        c2, inds2 = c.select([5.6, 5.6], return_index=True, outer=False)
+        assert inds2 == slice(6, 7)
+
+        c3, inds3 = c.select([5.4, 5.4], return_index=True, outer=False)
+        assert inds3 == slice(5, 6)
+
+        c3, inds3b = c.select([5.4, 5.6], return_index=True, outer=False)
+        assert inds3b == slice(5, 7)
+
+        c4, inds4 = c.select([9.1, 9.1], return_index=True, outer=False)
+        assert inds4 == []
+
+        c5, inds5 = c.select([-0.1, -0.1], return_index=True, outer=False)
+        assert inds5 == []
+
 
 class TestUniformCoordinatesMethods(object):
     def test_unique(self):
diff --git a/podpac/core/coordinates/uniform_coordinates1d.py b/podpac/core/coordinates/uniform_coordinates1d.py
index f972c2d1..3f7ae039 100644
--- a/podpac/core/coordinates/uniform_coordinates1d.py
+++ b/podpac/core/coordinates/uniform_coordinates1d.py
@@ -530,8 +530,14 @@ def _select(self, bounds, return_index, outer):
         imin = np.clip(imin, 0, self.size)
 
         # empty case
-        if imin >= imax:
+        if imin > imax:
             return self._select_empty(return_index)
+        if imax == imin:
+            # could have been selected between two existing coordinates
+            imin = int(np.round(fmin))
+            imax = int(np.round(fmax)) + 1
+            if imin >= (self.size - 1) | imin < 0:
+                return self._select_empty(return_index)
 
         if self.is_descending:
             imax, imin = self.size - imin, self.size - imax
diff --git a/podpac/core/coordinates/utils.py b/podpac/core/coordinates/utils.py
index ef4f0b8d..b0639262 100644
--- a/podpac/core/coordinates/utils.py
+++ b/podpac/core/coordinates/utils.py
@@ -20,6 +20,11 @@
 from six import string_types
 import pyproj
 
+from lazy_import import lazy_function
+geodesic = lazy_function("geopy.distance.geodesic")
+
+import podpac
+
 
 def get_timedelta(s):
     """
@@ -602,3 +607,35 @@ def has_alt_units(crs):
     with warnings.catch_warnings():
         warnings.simplefilter("ignore")
         return crs.is_vertical or "vunits" in crs.to_dict() or any(axis.direction == "up" for axis in crs.axis_info)
+
+
+def calculate_distance(point1, point2, ellipsoid_tuple, coordinate_name, units="meter"):
+    """Return distance of 2 points in desired unit measurement
+
+    Parameters
+    ----------
+    point1 : tuple
+    point2 : tuple
+
+    Returns
+    -------
+    float
+        The distance between point1 and point2, according to the current coordinate system's distance metric, using the desired units
+    """
+    if coordinate_name == "cartesian":
+        return np.linalg.norm(point1 - point2, axis=-1, units="meter") * podpac.units(units)
+    else:
+        if not isinstance(point1, tuple) and point1.size > 2:
+            distances = np.empty(len(point1))
+            for i in range(len(point1)):
+                distances[i] = geodesic(point1[i], point2[i], ellipsoid=ellipsoid_tuple).m
+            return distances * podpac.units("metre").to(podpac.units(units))
+        if not isinstance(point2, tuple) and point2.size > 2:
+            distances = np.empty(len(point2))
+            for i in range(len(point2)):
+                distances[i] = geodesic(point1, point2[i], ellipsoid=ellipsoid_tuple).m
+            return distances * podpac.units("metre").to(podpac.units(units))
+        else:
+            return (geodesic(point1, point2, ellipsoid=ellipsoid_tuple).m) * podpac.units("metre").to(
+                podpac.units(units)
+            )
diff --git a/podpac/core/interpolation/interpolation_manager.py b/podpac/core/interpolation/interpolation_manager.py
index daaa3f23..9ed2e92d 100644
--- a/podpac/core/interpolation/interpolation_manager.py
+++ b/podpac/core/interpolation/interpolation_manager.py
@@ -536,7 +536,14 @@ def select_coordinates(self, source_coordinates, eval_coordinates, index_type="n
                 if isinstance(selected_coords_idx[i], tuple):
                     selected_coords_idx2.extend(selected_coords_idx[i])
                 else:
-                    selected_coords_idx2.append(selected_coords_idx[i])
+                    if isinstance(selected_coords_idx[i], np.ndarray):
+                        # This happens when the interpolator_queue is empty, so we have to turn the
+                        # initialized coordinates into slices instead of numpy arrays
+                        selected_coords_idx2.append(
+                            slice(selected_coords_idx[i].min(), selected_coords_idx[i].max() + 1)
+                        )
+                    else:
+                        selected_coords_idx2.append(selected_coords_idx[i])
 
             selected_coords_idx2 = tuple(selected_coords_idx2)
         else:
diff --git a/podpac/core/node.py b/podpac/core/node.py
index d04fa4f9..739ab0ea 100644
--- a/podpac/core/node.py
+++ b/podpac/core/node.py
@@ -206,8 +206,8 @@ def __init__(self, **kwargs):
                         self.set_trait(name, tkwargs.pop(name))
                     trait.read_only = True
 
-            # Call traitlets constructor
-            super(Node, self).__init__(**tkwargs)
+        # Call traitlets constructor
+        super(Node, self).__init__(**tkwargs)
 
         self._traits_initialized_guard = True
 
@@ -843,71 +843,29 @@ def from_definition(cls, definition):
 
         # parse node definitions in order
         nodes = OrderedDict()
+        output_node = None
         for name, d in definition.items():
+            if name == "podpac_output_node":
+                output_node = d
+                continue
             if name == "podpac_version":
                 continue
 
             if "node" not in d:
                 raise ValueError("Invalid definition for node '%s': 'node' property required" % name)
 
-            # get node class
-            module_root = d.get("plugin", "podpac")
-            node_string = "%s.%s" % (module_root, d["node"])
-            module_name, node_name = node_string.rsplit(".", 1)
-            try:
-                module = importlib.import_module(module_name)
-            except ImportError:
-                raise ValueError("Invalid definition for node '%s': no module found '%s'" % (name, module_name))
-            try:
-                node_class = getattr(module, node_name)
-            except AttributeError:
-                raise ValueError(
-                    "Invalid definition for node '%s': class '%s' not found in module '%s'"
-                    % (name, node_name, module_name)
-                )
-
-            # parse and configure kwargs
-            kwargs = {}
-            for k, v in d.get("attrs", {}).items():
-                kwargs[k] = v
-
-            for k, v in d.get("inputs", {}).items():
-                kwargs[k] = _lookup_input(nodes, name, v)
-
-            for k, v in d.get("lookup_attrs", {}).items():
-                kwargs[k] = _lookup_attr(nodes, name, v)
-
-            if "style" in d:
-                style_class = getattr(node_class, "style", Style)
-                if isinstance(style_class, tl.TraitType):
-                    # Now we actually have to look through the class to see
-                    # if there is a custom initializer for style
-                    for attr in dir(node_class):
-                        atr = getattr(node_class, attr)
-                        if not isinstance(atr, tl.traitlets.DefaultHandler) or atr.trait_name != "style":
-                            continue
-                        try:
-                            style_class = atr(node_class)
-                        except Exception as e:
-                            # print ("couldn't make style from class", e)
-                            try:
-                                style_class = atr(node_class())
-                            except:
-                                # print ("couldn't make style from class instance", e)
-                                style_class = style_class.klass
-                try:
-                    kwargs["style"] = style_class.from_definition(d["style"])
-                except Exception as e:
-                    kwargs["style"] = Style.from_definition(d["style"])
-                    # print ("couldn't make style from inferred style class", e)
-
-            for k in d:
-                if k not in ["node", "inputs", "attrs", "lookup_attrs", "plugin", "style"]:
-                    raise ValueError("Invalid definition for node '%s': unexpected property '%s'" % (name, k))
+            _process_kwargs(name, d, definition, nodes)
 
-            nodes[name] = node_class(**kwargs)
+        # look for podpac_output_node attribute
+        if output_node is None:
+            return list(nodes.values())[-1]
 
-        return list(nodes.values())[-1]
+        if output_node not in nodes:
+            raise ValueError(
+                "Invalid definition for value 'podpac_output_node': reference to nonexistent node '%s' in lookup_attrs"
+                % (output_node)
+            )
+        return nodes[output_node]
 
     @classmethod
     def from_json(cls, s):
@@ -1205,13 +1163,36 @@ def get_ui_spec(cls, help_as_html=False):
         return spec
 
 
-def _lookup_input(nodes, name, value):
+def _lookup_input(nodes, name, value, definition):
+    """check if inputs of a node are stored in nodes, if not add them
+
+    Parameters
+    -----------
+    nodes: OrderedDict
+        Keys: Node names (strings)
+        Values: Node objects
+
+    name: string
+        the Node whose inputs are being checked
+
+    value: string, list, dictionary:
+        the Node (or collection of Nodes) which is being looked
+
+    definition: pipeline definition
+
+    Returns
+    --------
+    node: the node searched for
+
+    Note: this function calles _process_kwargs, which alters nodes by loading a Node if it is not yet in nodes.
+
+    """
     # containers
     if isinstance(value, list):
-        return [_lookup_input(nodes, name, elem) for elem in value]
+        return [_lookup_input(nodes, name, elem, definition) for elem in value]
 
     if isinstance(value, dict):
-        return {k: _lookup_input(nodes, name, v) for k, v in value.items()}
+        return {k: _lookup_input(nodes, name, v, definition) for k, v in value.items()}
 
     # node reference
     if not isinstance(value, six.string_types):
@@ -1219,12 +1200,21 @@ def _lookup_input(nodes, name, value):
             "Invalid definition for node '%s': invalid reference '%s' of type '%s' in inputs"
             % (name, value, type(value))
         )
+    # node not yet discovered yet
+    if not value in nodes:
+        # Look for it in the definition items:
+        for found_name, d in definition.items():
+            if value != found_name:
+                continue
+            # Load the node into nodes
+            _process_kwargs(found_name, d, definition, nodes)
+
+            break
 
     if not value in nodes:
         raise ValueError(
             "Invalid definition for node '%s': reference to nonexistent node '%s' in inputs" % (name, value)
         )
-
     node = nodes[value]
 
     # copy in debug mode
@@ -1240,7 +1230,7 @@ def _lookup_attr(nodes, name, value):
         return [_lookup_attr(nodes, name, elem) for elem in value]
 
     if isinstance(value, dict):
-        return {_k: _lookup_attr(nodes, name, v) for k, v in value.items()}
+        return {k: _lookup_attr(nodes, name, v) for k, v in value.items()}
 
     if not isinstance(value, six.string_types):
         raise ValueError(
@@ -1272,6 +1262,84 @@ def _lookup_attr(nodes, name, value):
     return attr
 
 
+def _process_kwargs(name, d, definition, nodes):
+    """create a node and add it to nodes
+
+    Parameters
+    -----------
+    nodes: OrderedDict
+        Keys: Node names (strings)
+        Values: Node objects
+
+    name: string
+        the Node which will be created
+
+    d: the definition of the node to be created
+
+    definition: pipeline definition
+
+    Returns
+    --------
+    Nothing, but loads the node with name "name" and definition "d" into nodes
+
+
+    """
+    # get node class
+    module_root = d.get("plugin", "podpac")
+    node_string = "%s.%s" % (module_root, d["node"])
+    module_name, node_name = node_string.rsplit(".", 1)
+    try:
+        module = importlib.import_module(module_name)
+    except ImportError:
+        raise ValueError("Invalid definition for node '%s': no module found '%s'" % (name, module_name))
+    try:
+        node_class = getattr(module, node_name)
+    except AttributeError:
+        raise ValueError(
+            "Invalid definition for node '%s': class '%s' not found in module '%s'" % (name, node_name, module_name)
+        )
+
+    kwargs = {}
+    for k, v in d.get("attrs", {}).items():
+        kwargs[k] = v
+
+    for k, v in d.get("inputs", {}).items():
+        kwargs[k] = _lookup_input(nodes, name, v, definition)
+
+    for k, v in d.get("lookup_attrs", {}).items():
+        kwargs[k] = _lookup_attr(nodes, name, v)
+
+    if "style" in d:
+        style_class = getattr(node_class, "style", Style)
+        if isinstance(style_class, tl.TraitType):
+            # Now we actually have to look through the class to see
+            # if there is a custom initializer for style
+            for attr in dir(node_class):
+                atr = getattr(node_class, attr)
+                if not isinstance(atr, tl.traitlets.DefaultHandler) or atr.trait_name != "style":
+                    continue
+                try:
+                    style_class = atr(node_class)
+                except Exception as e:
+                    # print ("couldn't make style from class", e)
+                    try:
+                        style_class = atr(node_class())
+                    except:
+                        # print ("couldn't make style from class instance", e)
+                        style_class = style_class.klass
+        try:
+            kwargs["style"] = style_class.from_definition(d["style"])
+        except Exception as e:
+            kwargs["style"] = Style.from_definition(d["style"])
+            # print ("couldn't make style from inferred style class", e)
+
+    for k in d:
+        if k not in ["node", "inputs", "attrs", "lookup_attrs", "plugin", "style"]:
+            raise ValueError("Invalid definition for node '%s': unexpected property '%s'" % (name, k))
+
+    nodes[name] = node_class(**kwargs)
+
+
 # --------------------------------------------------------#
 #  Mixins
 # --------------------------------------------------------#
diff --git a/podpac/core/test/test_node.py b/podpac/core/test/test_node.py
index 03853927..055c5972 100644
--- a/podpac/core/test/test_node.py
+++ b/podpac/core/test/test_node.py
@@ -1228,6 +1228,236 @@ def test_from_definition_version_warning(self):
         with pytest.warns(UserWarning, match="node definition version mismatch"):
             node = Node.from_json(s)
 
+    def test_from_proper_json(self):
+        not_ordered_json = """
+        {
+            "Arithmetic": {
+                "node": "core.algorithm.generic.Arithmetic",
+                "attrs": {
+                    "eqn": "a+b",
+                    "params": {
+
+                    }
+                },
+                "inputs": {
+                    "a": "SinCoords",
+                    "b": "Arange"
+                }
+            },
+            "SinCoords": {
+                "node": "core.algorithm.utility.SinCoords",
+                "style": {
+                    "colormap": "jet",
+                    "clim": [
+                        -1.0,
+                        1.0
+                    ]
+                }
+            },
+            "Arange": {
+                "node": "core.algorithm.utility.Arange"
+            },
+            "podpac_version": "3.2.0"
+        }
+        """
+        not_ordered_json_2 = """
+        {
+            "SinCoords": {
+                "node": "core.algorithm.utility.SinCoords",
+                "style": {
+                    "colormap": "jet",
+                    "clim": [
+                        -1.0,
+                        1.0
+                    ]
+                }
+            },
+            "Arithmetic": {
+                "node": "core.algorithm.generic.Arithmetic",
+                "attrs": {
+                    "eqn": "a+b",
+                    "params": {
+
+                    }
+                },
+                "inputs": {
+                    "a": "SinCoords",
+                    "b": "Arange"
+                }
+            },
+            "Arange": {
+                "node": "core.algorithm.utility.Arange"
+            },
+            "podpac_version": "3.2.0"
+        }
+        """
+        ordered_json = """
+        {
+            "SinCoords": {
+                "node": "core.algorithm.utility.SinCoords",
+                "style": {
+                    "colormap": "jet",
+                    "clim": [
+                        -1.0,
+                        1.0
+                    ]
+                }
+            },
+            "Arange": {
+                "node": "core.algorithm.utility.Arange"
+            },
+             "Arithmetic": {
+                "node": "core.algorithm.generic.Arithmetic",
+                "attrs": {
+                    "eqn": "a+b",
+                    "params": {
+
+                    }
+                },
+                "inputs": {
+                    "a": "SinCoords",
+                    "b": "Arange"
+                }
+            },
+            "podpac_version": "3.2.0"
+        }
+        """
+        # Check that the order doesn't matter. Because .from_json returns the output node, also checks correct output_node is returned
+        not_ordered_pipe = Node.from_json(not_ordered_json)
+        not_ordered_pipe_2 = Node.from_json(not_ordered_json_2)
+        ordered_pipe = Node.from_json(ordered_json)
+        assert not_ordered_pipe.definition == ordered_pipe.definition == not_ordered_pipe_2.definition
+        assert not_ordered_pipe.hash == ordered_pipe.hash
+
+        # Check that incomplete json will throw ValueError:
+        incomplete_json = """
+        {
+            "Arange": {
+                "node": "core.algorithm.utility.Arange"
+            },
+             "Arithmetic": {
+                "node": "core.algorithm.generic.Arithmetic",
+                "attrs": {
+                    "eqn": "a+b",
+                    "params": {
+
+                    }
+                },
+                "inputs": {
+                    "a": "SinCoords",
+                    "b": "Arange"
+                }
+            },
+            "podpac_version": "3.2.0"
+        }
+        """
+        with pytest.raises(ValueError):
+            Node.from_json(incomplete_json)
+
+    def test_output_node(self):
+        included_json = """
+        {
+            "Arithmetic": {
+                "node": "core.algorithm.generic.Arithmetic",
+                "attrs": {
+                    "eqn": "a+b",
+                    "params": {
+
+                    }
+                },
+                "inputs": {
+                    "a": "SinCoords",
+                    "b": "Arange"
+                }
+            },
+            "SinCoords": {
+                "node": "core.algorithm.utility.SinCoords",
+                "style": {
+                    "colormap": "jet",
+                    "clim": [
+                        -1.0,
+                        1.0
+                    ]
+                }
+            },
+            "Arange": {
+                "node": "core.algorithm.utility.Arange"
+            },
+            "podpac_version": "3.2.0",
+            "podpac_output_node": "Arithmetic"
+        }
+        """
+        ordered_json = """
+        {
+            "SinCoords": {
+                "node": "core.algorithm.utility.SinCoords",
+                "style": {
+                    "colormap": "jet",
+                    "clim": [
+                        -1.0,
+                        1.0
+                    ]
+                }
+            },
+            "Arange": {
+                "node": "core.algorithm.utility.Arange"
+            },
+             "Arithmetic": {
+                "node": "core.algorithm.generic.Arithmetic",
+                "attrs": {
+                    "eqn": "a+b",
+                    "params": {
+
+                    }
+                },
+                "inputs": {
+                    "a": "SinCoords",
+                    "b": "Arange"
+                }
+            },
+            "podpac_version": "3.2.0"
+        }
+        """
+        included_pipe = Node.from_json(included_json)
+        ordered_pipe = Node.from_json(ordered_json)
+        assert included_pipe.definition == ordered_pipe.definition
+        assert included_pipe.hash == ordered_pipe.hash
+
+        wrong_name_json = """
+        {
+            "SinCoords": {
+                "node": "core.algorithm.utility.SinCoords",
+                "style": {
+                    "colormap": "jet",
+                    "clim": [
+                        -1.0,
+                        1.0
+                    ]
+                }
+            },
+            "Arange": {
+                "node": "core.algorithm.utility.Arange"
+            },
+             "Arithmetic": {
+                "node": "core.algorithm.generic.Arithmetic",
+                "attrs": {
+                    "eqn": "a+b",
+                    "params": {
+
+                    }
+                },
+                "inputs": {
+                    "a": "SinCoords",
+                    "b": "Arange"
+                }
+            },
+            "podpac_version": "3.2.0",
+            "podpac_output_node": "Sum"
+        }
+        """
+        with pytest.raises(ValueError):
+            Node.from_json(wrong_name_json)
+
 
 class TestNoCacheMixin(object):
     class NoCacheNode(NoCacheMixin, Node):
diff --git a/setup.py b/setup.py
index 405c7283..ce510f01 100644
--- a/setup.py
+++ b/setup.py
@@ -26,6 +26,7 @@
     "lazy-import>=0.2.2",
     "psutil",
     "affine",
+    "geopy",
 ]
 
 if sys.version_info.major == 2:
@@ -38,7 +39,7 @@
         "beautifulsoup4>=4.6",
         "h5py>=2.9",
         "lxml>=4.2",
-        "pydap>=3.2",
+        "pydap>=3.3",
         "rasterio>=1.0",
         "zarr>=2.3",
         "owslib",