Repository: creare-com/podpac
Issue: "You will likely lose important projection information" when using using alt_units

**Description**

The `set_vunits` coordinates utility function triggers a pyproj warning about losing information with some versions of pyproj. The Coordinates `alt_units` requires this function.

**Steps to Reproduce**

Using python 3.8 and pyproj 2.4.2:

```
from podpac.core.coordinates.utils import set_vunits
set_vunits("+proj=merc +vunits=m", "ft")
```

You can also run the tests in `test_coordinates_utils.py::test_set_vunits` and many others.

**Observed Behavior**

```
UserWarning: You will likely lose important projection information when converting to a PROJ string from another format.
```

This occurs in both locations that `set_vunits` calls `to_proj4()`.

**Expected Behavior**

The `+vunits` should be set without any warnings.

**Additional Notes**

I think the problem is that the way we handle the `alt_units` parameter uses proj4 strings internally, which are lossy. We should either, in order of preference:
 - fix `set_vunits` to avoid losing projection information, if possible
 - catch the UserWarnings in `set_vunits` and replace them with our own warning that using `alt_units` may be lossy. We should also point to an alternative.
 - remove the `alt_units` parameter. The user would have to use the `crs` string to specify altitude units.

Changed files:
podpac/core/coordinates/test/test_coordinates_utils.py, podpac/core/coordinates/utils.py, podpac/core/coordinates/coordinates.py, podpac/core/coordinates/test/test_coordinates.py

diff --git a/podpac/core/coordinates/coordinates.py b/podpac/core/coordinates/coordinates.py
index 60663c64..3c369bf9 100644
--- a/podpac/core/coordinates/coordinates.py
+++ b/podpac/core/coordinates/coordinates.py
@@ -24,7 +24,6 @@ import pyproj
 import podpac
 from podpac.core.settings import settings
 from podpac.core.utils import OrderedDictTrait, _get_query_params_from_url, _get_param
-from podpac.core.coordinates.utils import get_vunits, set_vunits, rem_vunits
 from podpac.core.coordinates.base_coordinates import BaseCoordinates
 from podpac.core.coordinates.coordinates1d import Coordinates1d
 from podpac.core.coordinates.array_coordinates1d import ArrayCoordinates1d
@@ -64,11 +63,10 @@ class Coordinates(tl.HasTraits):
     """
 
     crs = tl.Unicode(read_only=True, allow_none=True)
-    alt_units = tl.Unicode(read_only=True, allow_none=True, default_value=None)
 
     _coords = OrderedDictTrait(trait=tl.Instance(BaseCoordinates), default_value=OrderedDict())
 
-    def __init__(self, coords, dims=None, crs=None, alt_units=None, ctype=None):
+    def __init__(self, coords, dims=None, crs=None, ctype=None):
         """
         Create multidimensional coordinates.
 
@@ -87,10 +85,7 @@ class Coordinates(tl.HasTraits):
              * 'lat', 'lon', 'alt', or 'time' for unstacked coordinates
              * dimension names joined by an underscore for stacked coordinates
         crs : str, optional
-            Coordinate reference system. Supports any PROJ4 compliant string (https://proj4.org/index.html).
-        alt_units : str, optional
-            Altitude units. Supports any `PROJ4 Distance Units
-            <https://proj4.org/operations/conversions/unitconvert.html#distance-units>`. Must not contradict the crs.
+            Coordinate reference system. Supports any PROJ4 or PROJ6 compliant string (https://proj.org).
         ctype : str, optional
             Default coordinates type. One of 'point', 'midpoint', 'left', 'right'.
         """
@@ -144,8 +139,6 @@ class Coordinates(tl.HasTraits):
         self.set_trait("_coords", dcoords)
         if crs is not None:
             self.set_trait("crs", crs)
-        if alt_units is not None:
-            self.set_trait("alt_units", alt_units)
         super(Coordinates, self).__init__()
 
     @tl.validate("_coords")
@@ -173,38 +166,17 @@ class Coordinates(tl.HasTraits):
     @tl.validate("crs")
     def _validate_crs(self, d):
         val = d["value"]
-        pyproj.CRS(val)  # raises pyproj.CRSError if invalid
+        CRS = pyproj.CRS(val)  # raises pyproj.CRSError if invalid
+
+        # make sure CRS defines vertical units
+        if "alt" in self.udims and not CRS.is_vertical:
+            raise ValueError("Altitude dimension is defined, but CRS does not contain vertical unit")
+
         return val
 
     @tl.observe("crs")
     def _observe_crs(self, d):
         crs = d["new"]
-        self.set_trait("alt_units", get_vunits(crs))
-
-    @tl.validate("alt_units")
-    def _validate_alt_units(self, d):
-        val = d["value"]
-        if val is None:
-            return None
-
-        # check if the alt_units are valid by trying to set the vunits
-        try:
-            set_vunits(self.crs, val)
-        except pyproj.crs.CRSError:
-            raise ValueError("Invalid alt_units '%s', alt_units must be PROJ4 compliant distance units." % val)
-
-        # check if the alt_units contradict the crs
-        # this will only matter if a full proj4 string has been supplied
-        vunits = get_vunits(self.crs)
-        if vunits is not None and val != vunits:
-            raise ValueError("crs and alt_units mismatch, '%s' conflicts with crs '%s'" % (val, self.crs))
-
-        vunits = get_vunits(set_vunits(self.crs, val))
-        if vunits is None:
-            warnings.warn("alt_units ignored (crs '%s' does not support separate vunits)" % self.crs)
-            return None
-
-        return val
 
     # ------------------------------------------------------------------------------------------------------------------
     # Alternate constructors
@@ -235,7 +207,7 @@ class Coordinates(tl.HasTraits):
         return coords
 
     @classmethod
-    def grid(cls, dims=None, crs=None, alt_units=None, ctype=None, **kwargs):
+    def grid(cls, dims=None, crs=None, ctype=None, **kwargs):
         """
         Create a grid of coordinates.
 
@@ -265,10 +237,7 @@ class Coordinates(tl.HasTraits):
             List of dimension names, must match the provided keyword arguments. In Python 3.6 and above, the ``dims``
             argument is optional, and the dims will match the order of the provided keyword arguments.
         crs : str, optional
-            Coordinate reference system. Supports any PROJ4 compliant string (https://proj4.org/index.html).
-        alt_units : str, optional
-            Altitude units. Supports any `PROJ4 Distance Units
-            <https://proj4.org/operations/conversions/unitconvert.html#distance-units>`. Must not contradict the crs.
+            Coordinate reference system. Supports any PROJ4 or PROJ6 compliant string (https://proj.org).
         ctype : str, optional
             Default coordinates type. One of 'point', 'midpoint', 'left', 'right'.
 
@@ -283,10 +252,10 @@ class Coordinates(tl.HasTraits):
         """
 
         coords = cls._coords_from_dict(kwargs, order=dims)
-        return cls(coords, crs=crs, ctype=ctype, alt_units=alt_units)
+        return cls(coords, crs=crs, ctype=ctype)
 
     @classmethod
-    def points(cls, crs=None, alt_units=None, ctype=None, dims=None, **kwargs):
+    def points(cls, crs=None, ctype=None, dims=None, **kwargs):
         """
         Create a list of multidimensional coordinates.
 
@@ -319,10 +288,7 @@ class Coordinates(tl.HasTraits):
             List of dimension names, must match the provided keyword arguments. In Python 3.6 and above, the ``dims``
             argument is optional, and the dims will match the order of the provided keyword arguments.
         crs : str, optional
-            Coordinate reference system. Supports any PROJ4 compliant string (https://proj4.org/index.html).
-        alt_units : str, optional
-            Altitude units. Supports any `PROJ4 Distance Units
-            <https://proj4.org/operations/conversions/unitconvert.html#distance-units>`. Must not contradict the crs.
+            Coordinate reference system. Supports any PROJ4 or PROJ6 compliant string (https://proj.org/).
         ctype : str, optional
             Default coordinates type. One of 'point', 'midpoint', 'left', 'right'.
 
@@ -338,10 +304,10 @@ class Coordinates(tl.HasTraits):
 
         coords = cls._coords_from_dict(kwargs, order=dims)
         stacked = StackedCoordinates(coords)
-        return cls([stacked], crs=crs, ctype=ctype, alt_units=alt_units)
+        return cls([stacked], crs=crs, ctype=ctype)
 
     @classmethod
-    def from_xarray(cls, xcoord, crs=None, alt_units=None, ctype=None):
+    def from_xarray(cls, xcoord, crs=None, ctype=None):
         """
         Create podpac Coordinates from xarray coords.
 
@@ -350,10 +316,7 @@ class Coordinates(tl.HasTraits):
         xcoord : xarray.core.coordinates.DataArrayCoordinates
             xarray coords
         crs : str, optional
-            Coordinate reference system. Supports any PROJ4 compliant string (https://proj4.org/index.html).
-        alt_units : str, optional
-            Altitude units. Supports any `PROJ4 Distance Units
-            <https://proj4.org/operations/conversions/unitconvert.html#distance-units>`. Must not contradict the crs.
+            Coordinate reference system. Supports any PROJ4 or PROJ6 compliant string (https://proj.org/).
         ctype : str, optional
             Default coordinates type. One of 'point', 'midpoint', 'left', 'right'.
 
@@ -374,7 +337,7 @@ class Coordinates(tl.HasTraits):
                 c = StackedCoordinates.from_xarray(xcoord[dim])
             coords.append(c)
 
-        return cls(coords, crs=crs, alt_units=alt_units, ctype=ctype)
+        return cls(coords, crs=crs, ctype=ctype)
 
     @classmethod
     def from_json(cls, s):
@@ -767,21 +730,27 @@ class Coordinates(tl.HasTraits):
 
     @property
     def CRS(self):
-        if self.alt_units is not None and self.alt_units != get_vunits(self.crs):
-            crs = set_vunits(self.crs, self.alt_units)
-        else:
-            crs = self.crs
-
+        crs = self.crs
         return pyproj.CRS(crs)
 
+    # TODO: add a convience property for displaying altitude units for the CRS
+    # @property
+    # def alt_units(self):
+    #     CRS = self.CRS
+
+    #     if CRS.is_vertical:
+    #         alt_units = <unsure how to get this from CRS>
+    #     else:
+    #         raise ValueError("CRS does not contain vertical component")
+
+    #     return alt_units
+
     @property
     def properties(self):
         """:dict: Dictionary of the coordinate properties. """
 
         d = OrderedDict()
         d["crs"] = self.crs
-        if self.alt_units is not None and self.alt_units != get_vunits(self.crs):
-            d["alt_units"] = self.alt_units
         return d
 
     @property
@@ -799,7 +768,11 @@ class Coordinates(tl.HasTraits):
 
         d = OrderedDict()
         d["coords"] = [c.full_definition for c in self._coords.values()]
-        d["crs"] = self.CRS.to_proj4()
+        d[
+            "crs"
+        ] = (
+            self.CRS.to_wkt()
+        )  # use "wkt" which is suggested as best format: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems
         return d
 
     @property
@@ -1169,12 +1142,12 @@ class Coordinates(tl.HasTraits):
         else:
             return Coordinates([self._coords[dim] for dim in dims], **self.properties)
 
-    def transform(self, crs=None, alt_units=None):
+    def transform(self, crs=None):
         """
         Transform coordinate dimensions (`lat`, `lon`, `alt`) into a different coordinate reference system.
-        Uses PROJ4 syntax for coordinate reference systems and units.
+        Uses PROJ syntax for coordinate reference systems and units.
 
-        See `PROJ4 Documentation <https://proj4.org/usage/projections.html#cartographic-projection>`_ for
+        See `PROJ Documentation <https://proj.org/usage/projections.html#cartographic-projection>`_ for
         more information about creating PROJ4 strings. See `PROJ4 Distance Units
         <https://proj4.org/operations/conversions/unitconvert.html#distance-units>`_ for unit string
         references.
@@ -1212,41 +1185,11 @@ class Coordinates(tl.HasTraits):
                 lat: ArrayCoordinates1d(lat): Bounds[-1847545.541169525, -615848.513723175], N[21], ctype['midpoint']
                 lon: ArrayCoordinates1d(lon): Bounds[-614897.0725896168, 614897.0725896184], N[21], ctype['midpoint']
 
-        Transform coordinates with altitude::
-
-            # include alt units in proj4 string
-            c = Coordinates([[0, 1, 2], [0, 1, 2], [1, 2, 3]], dims=['lat', 'lon', 'alt'])
-            c.transform('+proj=tmerc +vunits=ft')
-
-            >> Coordinates
-                lat: ArrayCoordinates1d(lat): Bounds[594971.8894642257, 819117.0608407748], N[3], ctype['midpoint']
-                lon: ArrayCoordinates1d(lon): Bounds[29772096.71234478, 29995929.885877542], N[3], ctype['midpoint']
-                alt: ArrayCoordinates1d(alt): Bounds[3.280839895013123, 9.842519685039369], N[3], ctype['midpoint']
-
-
-            # specify alt units seperately
-            c.transform('EPSG:2193', alt_units='ft')
-
-            >> Coordinates
-                lat: ArrayCoordinates1d(lat): Bounds[594971.8894642257, 819117.0608407748], N[3], ctype['midpoint']
-                lon: ArrayCoordinates1d(lon): Bounds[29772096.71234478, 29995929.885877542], N[3], ctype['midpoint']
-                alt: ArrayCoordinates1d(alt): Bounds[3.280839895013123, 9.842519685039369], N[3], ctype['midpoint']
-
-            # using alt_units will save the property `crs` as a proj4 string:
-            ct = c.transform('EPSG:2193', alt_units='ft')
-            ct.crs
-
-            >> '+proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +units=m +no_defs +type=crs +vunits=ft'
-
         Parameters
         ----------
         crs : str, optional
             PROJ4 compatible coordinate reference system string.
             Defaults to the current `crs`
-        alt_units : str, optional
-            Override the alt units defined in `crs` string.
-            This is implemented to provide a shorthand for transforming alt units
-            without specifying the whole proj4 string.
 
         Returns
         -------
@@ -1259,23 +1202,22 @@ class Coordinates(tl.HasTraits):
             Coordinates must have both lat and lon dimensions if either is defined
         """
 
-        if crs is None and alt_units is None:
-            raise TypeError("transform requires crs and/or alt_units argument")
+        if crs is None:
+            raise TypeError("transform requires crs argument")
 
         input_crs = crs
-        input_alt_units = alt_units
 
         # use self.crs by default
         if crs is None:
             crs = self.crs
 
-        # combine crs and alt_units
-        if alt_units is not None:
-            crs = set_vunits(crs, alt_units)
-
         from_crs = self.CRS
         to_crs = pyproj.CRS(crs)
 
+        # make sure to CRS defines vertical units
+        if "alt" in self.udims and not to_crs.is_vertical:
+            raise ValueError("Altitude dimension is defined, but CRS to transform does not contain vertical unit")
+
         # no transform needed
         if from_crs == to_crs:
             return deepcopy(self)
@@ -1283,8 +1225,8 @@ class Coordinates(tl.HasTraits):
         cs = [c for c in self.values()]
 
         # if lat-lon transform is required, check dims and convert unstacked lat-lon coordinates if necessary
-        from_spatial = pyproj.CRS(rem_vunits(self.crs))
-        to_spatial = pyproj.CRS(rem_vunits(crs))
+        from_spatial = pyproj.CRS(self.crs)
+        to_spatial = pyproj.CRS(crs)
         if from_spatial != to_spatial:
             if "lat" in self.dims and "lon" in self.dims:
                 ilat = self.dims.index("lat")
@@ -1318,7 +1260,7 @@ class Coordinates(tl.HasTraits):
         # transform
         transformer = pyproj.Transformer.from_proj(from_crs, to_crs, always_xy=True)
         ts = [c._transform(transformer) for c in cs]
-        return Coordinates(ts, crs=input_crs, alt_units=input_alt_units)
+        return Coordinates(ts, crs=input_crs)
 
     # ------------------------------------------------------------------------------------------------------------------
     # Operators/Magic Methods
diff --git a/podpac/core/coordinates/test/test_coordinates.py b/podpac/core/coordinates/test/test_coordinates.py
index cb1022fc..945f9927 100644
--- a/podpac/core/coordinates/test/test_coordinates.py
+++ b/podpac/core/coordinates/test/test_coordinates.py
@@ -404,70 +404,45 @@ class TestCoordinateCreation(object):
         # default
         c = Coordinates([lat, lon])
         assert c.crs == podpac.settings["DEFAULT_CRS"]
-        assert c.alt_units is None
         assert set(c.properties.keys()) == {"crs"}
 
         # crs
         c = Coordinates([lat, lon], crs="EPSG:2193")
         assert c.crs == "EPSG:2193"
-        assert c.alt_units is None
         assert set(c.properties.keys()) == {"crs"}
 
         # proj4
         c = Coordinates([lat, lon], crs="EPSG:2193")
         assert c.crs == "EPSG:2193"
-        assert c.alt_units is None
         assert set(c.properties.keys()) == {"crs"}
 
         c = Coordinates([lat, lon], crs="+proj=merc +lat_ts=56.5 +ellps=GRS80")
         assert c.crs == "+proj=merc +lat_ts=56.5 +ellps=GRS80"
-        assert c.alt_units is None
         assert set(c.properties.keys()) == {"crs"}
 
         # with vunits
         c = Coordinates([lat, lon], crs="+proj=merc +lat_ts=56.5 +ellps=GRS80 +vunits=ft")
         assert c.crs == "+proj=merc +lat_ts=56.5 +ellps=GRS80 +vunits=ft"
-        assert c.alt_units == "ft"
-        assert set(c.properties.keys()) == {"crs"}  # no alt_units, it is in the crs
+        assert set(c.properties.keys()) == {"crs"}
 
         # invalid
         with pytest.raises(pyproj.crs.CRSError):
             Coordinates([lat, lon], crs="abcd")
 
-    def test_alt_units(self):
+    def test_crs_with_vertical_units(self):
 
         alt = ArrayCoordinates1d([0, 1, 2], name="alt")
 
-        # None
-        c = Coordinates([alt])
-        assert c.alt_units is None
+        c = Coordinates([alt], crs="+proj=merc +vunits=m")
         assert set(c.properties.keys()) == {"crs"}
 
-        c.set_trait("alt_units", None)
-
-        # proj4
-        c = Coordinates([alt], crs="EPSG:2193", alt_units="ft")
-        assert c.alt_units == "ft"
-        assert set(c.properties.keys()) == {"crs", "alt_units"}
-
         # with crs
-        c = Coordinates([alt], crs="EPSG:2193", alt_units="ft")
-        assert c.crs == "EPSG:2193"
-        assert c.alt_units == "ft"
-        assert set(c.properties.keys()) == {"crs", "alt_units"}
+        ct = c.transform("+proj=merc +vunits=ft")
+        assert round(ct["alt"].coordinates[-1]) == 7.0  # transformed to feet
 
         # invalid
-        with pytest.raises(ValueError, match="Invalid alt_units"):
-            Coordinates([alt], alt_units="feet")
-
-        # crs mismatch
-        with pytest.raises(ValueError, match="crs and alt_units mismatch"):
-            Coordinates([alt], crs="+proj=merc +vunits=ft", alt_units="m")
-
-        # ignore
-        with pytest.warns(UserWarning, match="alt_units ignored"):
-            c = Coordinates([alt], crs="EPSG:4326", alt_units="ft")
-        assert c.alt_units is None
+        with pytest.raises(ValueError):
+            Coordinates([alt], crs="EPSG:2193")
 
     def test_ctype(self):
         # assign
@@ -521,13 +496,12 @@ class TestCoordinatesSerialization(object):
         lon = ArrayCoordinates1d([0, 1, 2], "lon")
 
         # default
-        c = Coordinates([lat, lon], crs="EPSG:2193", alt_units="ft")
+        c = Coordinates([lat, lon], crs="EPSG:2193")
         d = c.definition
         json.dumps(d, cls=podpac.core.utils.JSONEncoder)
         c2 = Coordinates.from_definition(d)
         assert c2 == c
         assert c2.crs == "EPSG:2193"
-        assert c2.alt_units == "ft"
 
     def test_from_definition(self):
         d = {
@@ -1214,28 +1188,8 @@ class TestCoordinatesMethods(object):
         assert round(c1_trans["alt"].coordinates[1]) == 656.0
         assert c1_trans.crs == proj
 
-        # make sure vunits can be overwritten appropriately
-        c2_trans = c1_trans.transform("EPSG:2193", alt_units="m")
-        assert round(c2_trans["alt"].coordinates[0]) == 100.0
-        assert c2_trans.alt_units == "m"
-
-        # alt_units parameter
-        with pytest.warns(UserWarning, match="transformation of coordinate segment lengths not yet implemented"):
-            c_trans = c.transform("EPSG:2193", alt_units="us-ft")
-        assert round(c_trans["alt"].coordinates[1]) == 3.0
-        assert round(c_trans["alt"].coordinates[2]) == 7.0
-        assert c_trans.crs == "EPSG:2193"
-        assert c_trans.alt_units == "us-ft"
-
-        with pytest.warns(UserWarning, match="transformation of coordinate segment lengths not yet implemented"):
-            c_trans = c.transform("EPSG:2193", alt_units="km")
-        assert c_trans["alt"].coordinates[1] == 0.001
-        assert c_trans["alt"].coordinates[2] == 0.002
-        assert c_trans.crs == "EPSG:2193"
-        assert c_trans.alt_units == "km"
-
         # no parameter
-        with pytest.raises(TypeError, match="transform requires crs and/or alt_units argument"):
+        with pytest.raises(TypeError, match="transform requires crs argument"):
             c.transform()
 
     def test_select_single(self):
@@ -1520,56 +1474,6 @@ class TestCoordinatesSpecial(object):
         assert c1.hash != c2.hash
         assert c2.hash == deepcopy(c2).hash
 
-    def test_eq_ne_hash_alt_units(self):
-        lat = [0, 1, 2]
-        lon = [10, 20, 30]
-        c1 = Coordinates([lat, lon], dims=["lat", "lon"], crs="EPSG:2193")
-        c2 = Coordinates([lat, lon], dims=["lat", "lon"], crs="EPSG:2193", alt_units="ft")
-        c3 = Coordinates([lat, lon], dims=["lat", "lon"], crs="EPSG:2193", alt_units="m")
-
-        # eq
-        assert not c1 == c2
-        assert not c1 == c3
-        assert not c2 == c3
-
-        assert c1 == deepcopy(c1)
-        assert c2 == deepcopy(c2)
-        assert c3 == deepcopy(c3)
-
-        # ne (this only matters in python 2)
-        assert c1 != c2
-        assert c1 != c3
-        assert c2 != c3
-
-        assert not c1 != deepcopy(c1)
-        assert not c2 != deepcopy(c2)
-        assert not c3 != deepcopy(c3)
-
-        # hash
-        assert c1.hash != c2.hash
-        assert c1.hash != c3.hash
-        assert c2.hash != c3.hash
-
-        assert c1.hash == deepcopy(c1).hash
-        assert c2.hash == deepcopy(c2).hash
-        assert c3.hash == deepcopy(c3).hash
-
-    def test_eq_ne_hash_crs_alt_units(self):
-        lat = [0, 1, 2]
-        lon = [10, 20, 30]
-
-        # special case, these should be the same
-        c1 = Coordinates([lat, lon], dims=["lat", "lon"], crs="EPSG:2193", alt_units="ft")
-        c2 = Coordinates(
-            [lat, lon],
-            dims=["lat", "lon"],
-            crs=("+proj=tmerc +vunits=ft +ellps=GRS80 +lon_0=173 +k=0.9996" " +x_0=1600000 +y_0=10000000"),
-        )
-
-        assert c1 == c2
-        assert not c1 != c2
-        assert c1.hash == c2.hash
-
 
 class TestCoordinatesFunctions(object):
     def test_merge_dims(self):
diff --git a/podpac/core/coordinates/test/test_coordinates_utils.py b/podpac/core/coordinates/test/test_coordinates_utils.py
index 0c0e1f81..72b93168 100644
--- a/podpac/core/coordinates/test/test_coordinates_utils.py
+++ b/podpac/core/coordinates/test/test_coordinates_utils.py
@@ -8,7 +8,6 @@ from datetime import datetime
 from podpac.core.coordinates.utils import get_timedelta, get_timedelta_unit, make_timedelta_string
 from podpac.core.coordinates.utils import make_coord_value, make_coord_delta, make_coord_array, make_coord_delta_array
 from podpac.core.coordinates.utils import add_coord, divide_delta, divide_timedelta
-from podpac.core.coordinates.utils import get_vunits, set_vunits, rem_vunits
 
 
 def test_get_timedelta():
@@ -466,44 +465,3 @@ def test_divide_delta():
     assert divide_delta(np.timedelta64(1, "D"), 2) == np.timedelta64(12, "h")
     with pytest.raises(ValueError, match="Cannot divide timedelta .* evenly"):
         divide_delta(np.timedelta64(1, "D"), 17)
-
-
-def test_get_vunits():
-    assert get_vunits("+proj=merc +vunits=m") == "m"
-    assert get_vunits("+proj=merc +vunits=us-ft") == "us-ft"
-    assert get_vunits("+proj=merc +lat_ts=56.5 +vunits=ft +ellps=GRS80") == "ft"
-    assert get_vunits("+proj=merc +lat_ts=56.5 +ellps=GRS80") is None
-
-
-def test_set_vunits():
-    crs = set_vunits("+proj=merc +vunits=m", "ft")
-    assert get_vunits(crs) == "ft"
-
-    crs = set_vunits("+proj=merc +vunits=us-ft", "m")
-    assert get_vunits(crs) == "m"
-
-    crs = set_vunits("+proj=merc +lat_ts=56.5 +vunits=m +ellps=GRS80", "ft")
-    assert get_vunits(crs) == "ft"
-
-    crs = set_vunits("+proj=merc +lat_ts=56.5 +ellps=GRS80", "ft")
-    assert get_vunits(crs) == "ft"
-
-    crs = set_vunits("EPSG:2193", "ft")
-    assert get_vunits(crs) == "ft"
-
-    crs = set_vunits("EPSG:4326", "ft")
-    assert get_vunits(crs) is None
-
-
-def test_rem_vunits():
-    crs = rem_vunits("+proj=merc +vunits=m")
-    assert get_vunits(crs) is None
-
-    crs = rem_vunits("+proj=merc +vunits=us-ft")
-    assert get_vunits(crs) is None
-
-    crs = rem_vunits("+proj=merc +lat_ts=56.5 +vunits=ft +ellps=GRS80")
-    assert get_vunits(crs) is None
-
-    crs = rem_vunits("+proj=merc +lat_ts=56.5 +ellps=GRS80")
-    assert get_vunits(crs) is None
diff --git a/podpac/core/coordinates/utils.py b/podpac/core/coordinates/utils.py
index 1c15b665..27287756 100644
--- a/podpac/core/coordinates/utils.py
+++ b/podpac/core/coordinates/utils.py
@@ -463,74 +463,6 @@ class CoordinateType(tl.Enum):
         super(CoordinateType, self).__init__(["point", "left", "right", "midpoint"], *args, **kwargs)
 
 
-def get_vunits(crs):
-    """
-    Get vunits from a coordinate reference system string.
-
-    Arguments
-    ---------
-    crs : str
-        PROJ4 coordinate reference system.
-
-    Returns
-    -------
-    vunits : str
-        PROJ4 distance units for altitude, or None if no vunits present.
-    """
-
-    if "+vunits" not in crs:
-        return None
-
-    return re.search(r"(?<=\\+vunits=)[a-z\\-]+", crs).group(0)
-
-
-def set_vunits(crs, vunits):
-    """
-    Set the vunits of a coordinate reference system string. The vunits will be replaced or added, as needed.
-
-    Arguments
-    ---------
-    crs : str
-        PROJ4 coordinate reference system.
-    vunits : str
-        desired altitude units in PROJ4 distance units.
-
-    Returns
-    -------
-    crs : str
-        PROJ4 coordinate reference system string with the desired vunits.
-    """
-
-    if "+vunits" in crs:
-        crs = re.sub(r"(?<=\\+vunits=)[a-z\\-]+", vunits, crs)
-    else:
-        crs = pyproj.CRS(crs).to_proj4()  # convert EPSG-style strings
-        crs += " +vunits={}".format(vunits)
-
-    crs = pyproj.CRS(crs).to_proj4()  # standardize, this is optional
-
-    return crs
-
-
-def rem_vunits(crs):
-    """
-    Remove the vunits of a coordinate reference system string, if present.
-
-    Arguments
-    ---------
-    crs : str
-        PROJ4 coordinate reference system.
-
-    Returns
-    crs : str
-        PROJ4 coordinate referenc system without vunits.
-    """
-
-    if "+vunits" in crs:
-        crs = re.sub(r"\\+vunits=[a-z\\-]+", "", crs)
-    return crs
-
-
 def lower_precision_time_bounds(my_bounds, other_bounds, outer):
     """
     When given two bounds of np.datetime64, this function will convert both bounds to the lower-precision (in terms of