Skip to content

Commit

Permalink
Merge b6a3e57 into 0a7ba0b
Browse files Browse the repository at this point in the history
  • Loading branch information
@mlshapiro
mlshapiro committed Dec 31, 2019
2 parents 0a7ba0b + b6a3e57 commit a410641
Show file tree
Hide file tree
Showing 4 changed files with 56 additions and 320 deletions.
152 changes: 47 additions & 105 deletions podpac/core/coordinates/coordinates.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,6 @@
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
Expand Down Expand Up @@ -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.
Expand All @@ -87,10 +85,7 @@ def __init__(self, coords, dims=None, crs=None, alt_units=None, ctype=None):
* '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'.
"""
Expand Down Expand Up @@ -144,8 +139,6 @@ def __init__(self, coords, dims=None, crs=None, alt_units=None, ctype=None):
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")
Expand Down Expand Up @@ -173,38 +166,17 @@ def _default_crs(self):
@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
Expand Down Expand Up @@ -235,7 +207,7 @@ def _coords_from_dict(d, order=None):
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.
Expand Down Expand Up @@ -265,10 +237,7 @@ def grid(cls, dims=None, crs=None, alt_units=None, ctype=None, **kwargs):
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'.
Expand All @@ -283,10 +252,10 @@ def grid(cls, dims=None, crs=None, alt_units=None, ctype=None, **kwargs):
"""

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.
Expand Down Expand Up @@ -319,10 +288,7 @@ def points(cls, crs=None, alt_units=None, ctype=None, dims=None, **kwargs):
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'.
Expand All @@ -338,10 +304,10 @@ def points(cls, crs=None, alt_units=None, ctype=None, dims=None, **kwargs):

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.
Expand All @@ -350,10 +316,7 @@ def from_xarray(cls, xcoord, crs=None, alt_units=None, ctype=None):
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'.
Expand All @@ -374,7 +337,7 @@ def from_xarray(cls, xcoord, crs=None, alt_units=None, ctype=None):
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):
Expand Down Expand Up @@ -767,21 +730,27 @@ def coords(self):

@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
Expand All @@ -799,7 +768,11 @@ def full_definition(self):

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
Expand Down Expand Up @@ -1169,12 +1142,12 @@ def transpose(self, *dims, **kwargs):
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.
Expand Down Expand Up @@ -1212,41 +1185,11 @@ def transform(self, crs=None, alt_units=None):
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
-------
Expand All @@ -1259,32 +1202,31 @@ def transform(self, crs=None, alt_units=None):
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)

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")
Expand Down Expand Up @@ -1318,7 +1260,7 @@ def transform(self, crs=None, alt_units=None):
# 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
Expand Down

0 comments on commit a410641

Please sign in to comment.