Skip to content

Commit

Permalink
Merge 3f3e30f into fcbb909
Browse files Browse the repository at this point in the history
  • Loading branch information
mpu-creare committed Apr 20, 2020
2 parents fcbb909 + 3f3e30f commit e3c65af
Show file tree
Hide file tree
Showing 9 changed files with 397 additions and 115 deletions.
19 changes: 18 additions & 1 deletion podpac/core/coordinates/array_coordinates1d.py
Original file line number Diff line number Diff line change
Expand Up @@ -90,7 +90,7 @@ def __init__(self, coordinates, name=None, ctype=None, segment_lengths=None):
else:
self._is_monotonic = True
self._is_descending = self.coordinates[1] < self.coordinates[0]
self._is_uniform = np.allclose(deltas, deltas[0])
self._is_uniform = np.allclose(deltas, deltas[0], atol=1e-7)

# set common properties
super(ArrayCoordinates1d, self).__init__(name=name, ctype=ctype, segment_lengths=segment_lengths)
Expand Down Expand Up @@ -233,6 +233,23 @@ def copy(self):
kwargs = self.properties
return ArrayCoordinates1d(self.coordinates, **kwargs)

def simplify(self):
""" Get the simplified/optimized representation of these coordinates.
Returns
-------
simplified : ArrayCoordinates1d, UniformCoordinates1d
UniformCoordinates1d if the coordinates are uniform, otherwise ArrayCoordinates1d
"""

from podpac.core.coordinates.uniform_coordinates1d import UniformCoordinates1d

if self.is_uniform:
# return UniformCoordinates1d(self.start, self.stop, self.step, **self.properties)
return UniformCoordinates1d(self.coordinates[0], self.coordinates[-1], size=self.size, **self.properties)

return self

# ------------------------------------------------------------------------------------------------------------------
# standard methods, array-like
# ------------------------------------------------------------------------------------------------------------------
Expand Down
207 changes: 161 additions & 46 deletions podpac/core/coordinates/coordinates.py
Original file line number Diff line number Diff line change
Expand Up @@ -27,11 +27,13 @@
from podpac.core.utils import OrderedDictTrait, _get_query_params_from_url, _get_param
from podpac.core.coordinates.base_coordinates import BaseCoordinates
from podpac.core.coordinates.coordinates1d import Coordinates1d
from podpac.core.coordinates.dependent_coordinates import ArrayCoordinatesNd
from podpac.core.coordinates.array_coordinates1d import ArrayCoordinates1d
from podpac.core.coordinates.uniform_coordinates1d import UniformCoordinates1d
from podpac.core.coordinates.stacked_coordinates import StackedCoordinates
from podpac.core.coordinates.dependent_coordinates import DependentCoordinates
from podpac.core.coordinates.rotated_coordinates import RotatedCoordinates
from podpac.core.coordinates.cfunctions import clinspace

# Optional dependencies
from lazy_import import lazy_module, lazy_class
Expand Down Expand Up @@ -787,8 +789,7 @@ def coords(self):

@property
def CRS(self):
crs = self.crs
return pyproj.CRS(crs)
return pyproj.CRS(self.crs)

# TODO: add a convience property for displaying altitude units for the CRS
# @property
Expand Down Expand Up @@ -1265,7 +1266,7 @@ def transpose(self, *dims, **kwargs):
else:
return Coordinates(coords, validate_crs=False, **self.properties)

def transform(self, crs=None):
def transform(self, crs):
"""
Transform coordinate dimensions (`lat`, `lon`, `alt`) into a different coordinate reference system.
Uses PROJ syntax for coordinate reference systems and units.
Expand Down Expand Up @@ -1310,9 +1311,8 @@ def transform(self, crs=None):
Parameters
----------
crs : str, optional
crs : str
PROJ4 compatible coordinate reference system string.
Defaults to the current `crs`
Returns
-------
Expand All @@ -1324,16 +1324,6 @@ def transform(self, crs=None):
ValueError
Coordinates must have both lat and lon dimensions if either is defined
"""

if crs is None:
raise TypeError("transform requires crs argument")

input_crs = crs

# use self.crs by default
if crs is None:
crs = self.crs

from_crs = self.CRS
to_crs = pyproj.CRS(crs)

Expand All @@ -1345,45 +1335,170 @@ def transform(self, crs=None):
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(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")
ilon = self.dims.index("lon")
if ilat == ilon - 1:
c1, c2 = self["lat"], self["lon"]
elif ilon == ilat - 1:
c1, c2 = self["lon"], self["lat"]
else:
raise ValueError("Cannot transform coordinates with nonadjacent lat and lon, transpose first")
transformer = pyproj.Transformer.from_proj(from_crs, to_crs, always_xy=True)

# Let's test to see if the transformed coordinates can be simplified to either UniformCoordinates1D or
# ArrayCoordinates1D
def simplified_transform():
""" I'm wrapping this in a function so that I can escape early if the coordinates prove dependent
"""
# TODO: ADD ALTITUDE
TOL = 1e-7 # assume float 32 coordinate representation

def make_uniform(coords, size, name):
tc = None
dcoords = np.diff(coords)
if np.all(np.abs(dcoords[1:] - dcoords[0]) < (np.abs(dcoords[0]) * TOL)):
# We've actually already done the needed transformation, so we can just go ahead and create the
# coordinates
tc = clinspace(coords[0], coords[-1], size, name)
return tc

lat = clinspace(self["lat"].coordinates[0], self["lat"].coordinates[-1], 5)
lon = clinspace(self["lon"].coordinates[0], self["lon"].coordinates[-1], 5)
if "alt" in self.dims:
alt = clinspace(self["alt"].coordinates[0], self["alt"].coordinates[-1], 5)
c = DependentCoordinates(
np.meshgrid(c1.coordinates, c2.coordinates, indexing="ij"),
dims=[c1.name, c2.name],
ctypes=[c1.ctype, c2.ctype],
segment_lengths=[c1.segment_lengths, c2.segment_lengths],
np.meshgrid(lat.coordinates, lon.coordinates, alt.coordinates, indexing="ij"),
dims=["lat", "lon", "alt"],
ctypes=[self["lat"].ctype, self["lon"].ctype, self["alt"].ctype],
segment_lengths=[
self["lat"].segment_lengths,
self["lon"].segment_lengths,
self["alt"].segment_lengths,
],
)
else:
c = DependentCoordinates(
np.meshgrid(lat.coordinates, lon.coordinates, indexing="ij"),
dims=["lat", "lon"],
ctypes=[self["lat"].ctype, self["lon"].ctype],
segment_lengths=[self["lat"].segment_lengths, self["lon"].segment_lengths],
)
t = c._transform(transformer)
if not isinstance(t, list):
# Then we are not independent -- this was checked in the Dependent Coordinates
return

# Great, we CAN simplify the transformed coordinates. Now let's test if these are uniform or not
# Note, the dependent coordinates may already give UniformCoordinates, but these are not the correct size

# Lat uniform test and transform
if isinstance(t[0], UniformCoordinates1d):
t_lat = make_uniform(t[0].coordinates, self["lat"].size, "lat")
else:
# Need to compute the non-uniform transformed coordinates
temp_coords = Coordinates(
[[self["lat"].coordinates, self["lat"].coordinates * 0 + self["lon"].coordinates.mean()]],
["lat_lon"],
crs=self.crs,
)
t_lat = temp_coords.transform(crs)["lat"]

# There is a small possibility here that these are again uniform (i.e. non-inform array input gives
# uniform array output). So let's check that...
t_lat_2 = make_uniform(t_lat.coordinates, self["lat"].size, "lat")
if t_lat_2 is not None:
t_lat = t_lat_2

# Lon uniform test and transform
if isinstance(t[1], UniformCoordinates1d):
t_lon = make_uniform(t[1].coordinates, self["lon"].size, "lon")
else:
# Need to compute the non-uniform coordinates
temp_coords = Coordinates(
[[self["lon"].coordinates, self["lon"].coordinates * 0 + self["lat"].coordinates.mean()]],
["lon_lat"],
crs=self.crs,
)
t_lon = temp_coords.transform(crs)["lon"]

# There is a small possibility here that these are again uniform (i.e. non-inform array input gives
# uniform array output). So let's check that...
t_lon_2 = make_uniform(t_lon.coordinates, self["lon"].size, "lon")
if t_lon_2 is not None:
t_lon = t_lon_2

if len(t) < 3:
return (t_lat, t_lon)

# If we also have altitude, we haven't returned yet, so compute its transform
if isinstance(t[2], UniformCoordinates1d):
t_alt = make_uniform(t[2].coordinates, self["alt"].size, "alt")
else:
# Need to compute the non-uniform coordinates
temp_coords = Coordinates(
[
[
self["alt"].coordinates,
self["alt"].coordinates * 0 + self["lat"].coordinates.mean(),
self["alt"].coordinates * 0 + self["lon"].coordinates.mean(),
]
],
["alt_lon_lat"],
crs=self.crs,
)
t_alt = temp_coords.transform(crs)["alt"]

# There is a small possibility here that these are again uniform (i.e. non-inform array input gives
# uniform array output). So let's check that...
t_alt_2 = make_uniform(t_alt.coordinates, self["alt"].size, "alt")
if t_alt_2 is not None:
t_alt = t_alt_2

# replace 'lat' and 'lon' entries with single 'lat,lon' entry
i = min(ilat, ilon)
cs.pop(i)
cs.pop(i)
cs.insert(i, c)
return (t_lat, t_lon, t_alt)

elif "lat" in self.dims:
raise ValueError("Cannot transform lat coordinates without lon coordinates")
# Collect the individual coordinates
cs = [c for c in self.values()]

# if lat-lon transform is required, check dims and convert unstacked lat-lon coordinates if necessary
# Also, check to see if we can short-cut the transformation because the results are uniform coordinates
if "lat" in self.dims and "lon" in self.dims:
ilat = self.dims.index("lat")
ilon = self.dims.index("lon")
if ilat == ilon - 1:
c1, c2 = self["lat"], self["lon"]
elif ilon == ilat - 1:
c1, c2 = self["lon"], self["lat"]
else:
raise ValueError("Cannot transform coordinates with nonadjacent lat and lon, transpose first")

tc = simplified_transform()
if tc:
cs[ilat] = tc[0]
cs[ilon] = tc[1]
if "alt" in self.dims:
cs[self.dims.index("alt")] = tc[2]
return Coordinates(cs, crs=crs, validate_crs=False)

c = DependentCoordinates(
np.meshgrid(c1.coordinates, c2.coordinates, indexing="ij"),
dims=[c1.name, c2.name],
ctypes=[c1.ctype, c2.ctype],
segment_lengths=[c1.segment_lengths, c2.segment_lengths],
)

elif "lon" in self.dims:
raise ValueError("Cannot transform lon coordinates without lat coordinates")
# replace 'lat' and 'lon' entries with single 'lat,lon' entry
i = min(ilat, ilon)
cs.pop(i)
cs.pop(i)
cs.insert(i, c)

elif "lat" in self.dims:
raise ValueError("Cannot transform lat coordinates without lon coordinates")

elif "lon" in self.dims:
raise ValueError("Cannot transform lon coordinates without lat coordinates")

# 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, validate_crs=False)
ts = []
for c in cs:
tc = c._transform(transformer)
if isinstance(tc, list):
ts.extend(tc)
else:
ts.append(tc)
return Coordinates(ts, crs=crs, validate_crs=False)

# ------------------------------------------------------------------------------------------------------------------
# Operators/Magic Methods
Expand Down
11 changes: 11 additions & 0 deletions podpac/core/coordinates/coordinates1d.py
Original file line number Diff line number Diff line change
Expand Up @@ -293,6 +293,17 @@ def copy(self):

raise NotImplementedError

def simplify(self):
""" Get the simplified/optimized representation of these coordinates.
Returns
-------
simplified : Coordinates1d
simplified version of the coordinates
"""

raise NotImplementedError

def _select_empty(self, return_indices):
I = []
if return_indices:
Expand Down
21 changes: 19 additions & 2 deletions podpac/core/coordinates/dependent_coordinates.py
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@
from podpac.core.coordinates.coordinates1d import Coordinates1d
from podpac.core.coordinates.array_coordinates1d import ArrayCoordinates1d
from podpac.core.coordinates.stacked_coordinates import StackedCoordinates
from podpac.core.coordinates.cfunctions import clinspace


class DependentCoordinates(BaseCoordinates):
Expand Down Expand Up @@ -477,7 +478,9 @@ def _transform(self, transformer):
warnings.warn("transformation of coordinate segment lengths not yet implemented")
if self.ctypes[ialt] != "point":
_, _, tsl = transformer.transform(0, 0, self.segment_lengths[ialt])
properties["segment_lengths"][ialt] = tsl
new_prop = list(properties["segment_lengths"])
new_prop[ialt] = tsl
properties["segment_lengths"] = tuple(new_prop)

elif "lat" in self.dims and "lon" in self.dims:
ilat = self.dims.index("lat")
Expand Down Expand Up @@ -510,7 +513,21 @@ def _transform(self, transformer):
_, _, tsl = transformer.transform(0, 0, self.segment_lengths[ialt])
properties["segment_lengths"][ialt] = tsl

return DependentCoordinates(coords, **properties)
return DependentCoordinates(coords, **properties).simplify()

def simplify(self):
coords = [c.copy() for c in self.coordinates]
slc_start = [slice(0, 1) for d in self.dims]

for dim in self.dims:
i = self.dims.index(dim)
slc = slc_start.copy()
slc[i] = slice(None)
if not np.allclose(coords[i][tuple(slc)], coords[i], atol=1e-7):
return self
coords[i] = ArrayCoordinates1d(coords[i][tuple(slc)].squeeze(), name=dim).simplify()

return coords

def transpose(self, *dims, **kwargs):
"""
Expand Down
18 changes: 12 additions & 6 deletions podpac/core/coordinates/stacked_coordinates.py
Original file line number Diff line number Diff line change
Expand Up @@ -430,9 +430,11 @@ def _transform(self, transformer):
lon = coords[ilon]
alt = coords[ialt]
tlon, tlat, talt = transformer.transform(lon.coordinates, lat.coordinates, alt.coordinates)
coords[ilat].set_trait("coordinates", tlat)
coords[ilon].set_trait("coordinates", tlon)
coords[ialt].set_trait("coordinates", talt)

# Could check for uniform coordinates here, but a rare case.
coords[ilat] = ArrayCoordinates1d(tlat, "lat", ctype=lat.ctype, segment_lengths=lat.segment_lengths)
coords[ilon] = ArrayCoordinates1d(tlon, "lon", ctype=lon.ctype, segment_lengths=lon.segment_lengths)
coords[ialt] = ArrayCoordinates1d(talt, "alt", ctype=alt.ctype, segment_lengths=lon.segment_lengths)

# segment lengths
# TODO can we use the proj4 '+units' here, at least sometimes?
Expand All @@ -453,8 +455,10 @@ def _transform(self, transformer):
lat = coords[ilat]
lon = coords[ilon]
tlon, tlat = transformer.transform(lon.coordinates, lat.coordinates)
coords[ilat].set_trait("coordinates", tlat)
coords[ilon].set_trait("coordinates", tlon)

# Could check for uniform coordinates here, but a rare case.
coords[ilat] = ArrayCoordinates1d(tlat, "lat", ctype=lat.ctype, segment_lengths=lat.segment_lengths)
coords[ilon] = ArrayCoordinates1d(tlon, "lon", ctype=lon.ctype, segment_lengths=lon.segment_lengths)

# segment lengths
# TODO can we use proj4 '+units' here, at least sometimes?
Expand All @@ -469,7 +473,9 @@ def _transform(self, transformer):
# coordinates
alt = coords[ialt]
_, _, talt = transformer.transform(np.zeros(self.size), np.zeros(self.size), alt.coordinates)
coords[ialt].set_trait("coordinates", talt)

# Could check for uniform coordinates here, but a rare case.
coords[ialt] = ArrayCoordinates1d(talt, "alt", ctype=alt.ctype, segment_lengths=lon.segment_lengths)

# segment lengths
if alt.ctype != "point" and "segment_lengths" in lon.properties:
Expand Down
Loading

0 comments on commit e3c65af

Please sign in to comment.