Move segments from coordinates to datasource #395

jmilloy · 2020-04-17T17:37:48Z

Summary

coordinates ctype and segment_lengths are replaced by DataSource boundary
coordinates area_bounds is removed
rasterio reprojection bounds (and WCS request bounds) are calculated from the coordinates start, stop, and step
coordinates geotransform supports uniform ArrayCoordinates1d

Datasource coordinates boundary

The Datasource boundary is a dictionary of coordinate boundaries (one for each dimension). Point coordinates are omitted or None. The trait is validated and tested, but it is unused at this point and raises NotImplementedError in some cases. Note that all of these are coordinate offsets.

# simple centered boxes/ranges
boundary={'lat': 0.5, 'lon': 0.5}
boundary={'time': "1,D"}

# simple non-centered boxes/ranges
boundary={'lat': [0, 0.5], 'lon': [0, 0.5]}
boundary={'time': ["0,D", "1,D"]}

# polygons 
boundary={'lat': [0.0, -0.5, 0.0, 0.5], 'lon': [-0.5, 0.0, 0.5, 0.0]}

# array of boundaries, specified for each coordinate
boundary={"lat": [[-0.1, 0.4], [-0.2, 0.3], [-0.3, 0.2]], "lon": [[-0.1, 0.4], [-0.2, 0.3]]}

get_area_bounds

Area bounds are not currently used anywhere, but the logic to compute area bounds from coordinates and a boundary has been updated and may be useful. See get_area_bounds and associated tests. It supports all of the above boundary types.

The area_bounds property is removed as well (it depends on these). Note that ctype and segment_lengths are not used in any tests or core podpac code. The midpoint default for uniform coordinates is expected for some tests the rely on interpolation, which now fail.

coveralls · 2020-04-17T17:56:38Z

Coverage decreased (-0.04%) to 91.302% when pulling 0ca7e5c on feature/move-segments into 0fcf40f on develop.

…_bounds method, and uses both in interpolation.

…n docs and docstrings.

…e indexes non-uniform boundaries during evaluation as needed, but interpolation raises an exception for non-uniform or non-centered boundaries.

…s when necessary.

mpu-creare · 2020-04-19T02:41:48Z

#396 touches segments / ctypes in a few minor places. I think we should merge that FIRST, as this one probably has a lot of changes / updates.

jmilloy · 2020-04-19T03:29:42Z

I admit that I only have a cursory understanding, but I think that this is revealing to me the extent to which we are mis-using rasterio reprojection to do interpolation and getting lucky that our defaults usually work.

Basically, we infer the rasterio transform bounds from the coordinates within the coordinates object via the inferred segment_lengths. If the segment lengths cannot be inferred, or if the user supplies their own segment_lengths, the rasterio "interpolation" will be incorrect or raise an exception. (This is why you get an exception when you request coordinates with only one value in either dimension: we cannot infer the segment_lengths, so we silently fall back on the bounds which make an invalid geotransform because the high and low bonud is the same!).

I think I can make it work by inferring the area_bounds in the rasterio interpolate method, which will make everything a lot clearer. It will be clear what's actually going on, and I will be able to raise the right exceptions at that time.

Note that nothing else uses the soon-to-be-removed area_bounds. I believe that rasterio is unable to handle anything but centered rectangular coordinates, and that non-default coordinate boundaries (e.g. ctype left or non-uniform segment_lengths) must use a different interpolator.

mpu-creare · 2020-04-19T12:03:15Z

See if the new Coordinates.geotransform will serve the needed purpose in the rasterio interpolator.

Perhaps we can return False for the rasterio interpolator in cases where area_bounds cannot be inferred.

jmilloy · 2020-04-19T12:41:04Z

Perhaps we can return False for the rasterio interpolator in cases where area_bounds cannot be inferred.

Yes, the can_interpolate method is doing a better job that I realized.

jmilloy · 2020-04-19T12:47:45Z

See if the new Coordinates.geotransform will serve the needed purpose in the rasterio interpolator.

I haven't looked closely, but it is giving be the exception when I don't think it should.

self.coordinates.geotransform
*** TypeError: Only 2-D coordinates that are uniform or rotated have a GDAL transform. These coordinates Coordinates (EPSG:4326)
	lat: ArrayCoordinates1d(lat): Bounds[1.0, 3.0], N[3], ctype['midpoint']
	lon: ArrayCoordinates1d(lon): Bounds[1.0, 3.0], N[3], ctype['midpoint'] do not.

jmilloy · 2020-04-19T13:08:32Z

Ahh, you check isinstance(self._coord['lat'], UniformCoordinats1d) when you could check self._coord[lat].is_uniform.

…g the coordinates start, stop and step. * The boundary does not currently need to be passed through to the interpolator. * Adds start, stop, and step properties to the ArrayCoordinates1d (only valid when the coordinates are uniform). * Coordinates geotransform property supports uniform ArrayCoordinates1d. * Singleton ArrayCoordinates1d objects are not uniform. * WCS node also computes bounds using start, stop, and step. * Datasource boundary defaults are removed, validation and tests added.

podpac/core/coordinates/coordinates.py

podpac/core/coordinates/coordinates1d.py

podpac/core/coordinates/test/test_array_coordinates1d.py

jmilloy · 2020-04-19T16:30:28Z

podpac/core/data/datasource.py

+        # get indexed boundary
+        self._requested_source_boundary = self.get_boundary(self._requested_source_coordinates_index)
+


This call (and the get_boundary method) is only necessary for arrays of non-uniform boundaries. And since we aren't actually using the boundary anywhere in the interpolation, it isn't necessary at all. But I didn't want to forget about it. Otherwise, the datasource eval method is fundamentally unchanged.

jmilloy · 2020-04-19T16:32:03Z

podpac/core/interpolation/interpolators.py

+            west = c["lon"].bounds[0] - c["lon"].step / 2.0
+            east = c["lon"].bounds[1] + c["lon"].step / 2.0
+            south = c["lat"].bounds[0] - c["lat"].step / 2.0
+            north = c["lat"].bounds[1] + c["lat"].step / 2.0
+            cols = c["lon"].size
+            rows = c["lat"].size


@mpu-creare This is what we would replace c.geotransform if possible, but I couldn't figure out how exactly.

This is fine. As a one-liner, it would have been

return rasterio.Affine.from_gdal(*c.geotransform)

But this is pretty clear in case of future debugging...

jmilloy · 2020-04-19T16:36:22Z

@mpu-creare This is all set. Note that the interpolation tests are completed unchanged, I just broke them up into three files. Let me know if you to change anything in the boundary spec (it is similar but I think better to what we discussed in #395).

mpu-creare · 2020-04-20T12:38:24Z

Ahh, you check isinstance(self._coord['lat'], UniformCoordinats1d) when you could check self._coord[lat].is_uniform.

OH! That's good to know! I might have to fix that in #396 as well.

mpu-creare

Look pretty good. A few questions (see below).
Unit tests are failing.
I also want to play with it a little bit on my end.

podpac/core/compositor/tile_compositor.py

mpu-creare · 2020-04-20T15:39:22Z

podpac/core/coordinates/coordinates1d.py

+        ---------
+        boundary : float, timedelta, array, None
+            Boundary offsets in this dimension.
+            * For a centered uniform boundary, use a single positive float or timedelta offset


should this be step / 2 ?

otherwise in the offsets below you have to divide boundary by 2

I get it down. I'll make the documentation a little clearer.

mpu-creare · 2020-04-20T15:44:04Z

podpac/core/coordinates/test/test_array_coordinates1d.py

-        # uniform segment_lengths
-        c = ArrayCoordinates1d([1, 2, 4, 5], segment_lengths=0.5)
+        # polygon (i.e. there would be corresponding offets for another dimension)
+        area_bounds = c.get_area_bounds([-0.2, -0.5, 0.7, 0.5])


I don't quite understand this one...

I'm not sure if we will ever use this one, but it is very general.

For a centered square, you could use boundary = {'lat': [-1.0, -1.0, 1.0, 1.0], 'lon': [-1.0, 1.0, -1.0, 1.0]}. So for lat, you get boundary offsets [-1.0, -1.0, 1.0, 1.0], which appear a bit meaningless at that point.

For area bounds, all you care about is the min and max offset, and currently Coordinates1d get_area_bounds does that conversion (from boundary offsets to area bounds offsets). I guess it could be more clear if full Coordinates get_area_bounds method did the conversion, e.g:

class Coordinates: def get_area_bounds(self, boundary): offsets = {} for dim, value in boundary.items(): if np.ndarray(vaue).ndim == 0: value = make_coord_delta(value) offsets[dim] = [-value, value] elif np.ndarray(value).ndim == 1: value = make_coord_delta_array(value) offsets[dim] = [min(value), max(value)] else: value = np.array([make_coord_delta_array(elem) for elem in value]) offsets[dim] = np.array(np.min(value, axis=1), np.max(value, axis=1)).T return {dim: self[dim].get_area_bounds(offsets.get(dim)) for dim in self.udims}

And then Coordinates1d.get_area_bounds would expect offsets, and be a bit clearer:

class Coordinates1d: def get_area_bounds(self, offsets): # point coordinates if offsets is None: return self.bounds # empty coordinates if self.size == 0: return self.bounds if np.array(offsets).ndim == 1: # uniform boundary lo_offset, hi_offset = offsets else: # non-uniform boundary L, H = self.argbounds lo_offset = offsets[L][0] hi_offset = offsets[H][1] lo, hi = self.bounds lo = add_coord(lo, lo_offset) hi = add_coord(hi, hi_offset) return lo, hi

That is a little bit less efficient, but maybe worth it.

I'm not sure when/where we will need the area_bounds though.

Ah, this is great for future. Nice. Good spec.

podpac/core/coordinates/test/test_array_coordinates1d.py

podpac/core/data/datasource.py

mpu-creare · 2020-04-20T15:54:51Z

podpac/core/interpolation/interpolators.py

+            west = c["lon"].bounds[0] - c["lon"].step / 2.0
+            east = c["lon"].bounds[1] + c["lon"].step / 2.0
+            south = c["lat"].bounds[0] - c["lat"].step / 2.0
+            north = c["lat"].bounds[1] + c["lat"].step / 2.0
+            cols = c["lon"].size
+            rows = c["lat"].size


This is fine. As a one-liner, it would have been

return rasterio.Affine.from_gdal(*c.geotransform)

But this is pretty clear in case of future debugging...

jmilloy · 2020-04-20T16:50:28Z

return rasterio.Affine.from_gdal(*c.geotransform)

I think you can see that I tried that and left it commented it out. It is not the same and doesn't work.

…was renamed).

jmilloy · 2020-04-20T16:55:28Z

Okay fixed that typo in the tests.

mpu-creare · 2020-04-20T17:17:35Z

return rasterio.Affine.from_gdal(*c.geotransform)

I think you can see that I tried that and left it commented it out. It is not the same and doesn't work.

Yeah, I'll need to look into that because it should work. They are doing the same thing, or should be. This means one of the two are wrong.

mpu-creare · 2020-04-20T17:39:43Z

return rasterio.Affine.from_gdal(*c.geotransform)

I think you can see that I tried that and left it commented it out. It is not the same and doesn't work.

Yeah, I'll need to look into that because it should work. They are doing the same thing, or should be. This means one of the two are wrong.

Yeah, I see the issue. The coordinates.geotransform "naturally" handles the order of the dimensions (is_descending). For example, these give the same results:

coords_src = Coordinates([clinspace(10, 0, 3), clinspace(0, 10, 5)], dims=["lat", "lon"])
rasterio.Affine.from_gdal(*coords_src.geotransform)
Affine(2.5, 0.0, -1.25,
       0.0, -5.0, 12.5)
coords_src = Coordinates([clinspace(0, 10, 3), clinspace(0, 10, 5)], dims=["lat", "lon"])
get_rasterio_transform(coords_src)
Affine(2.5, 0.0, -1.25,
       0.0, -5.0, 12.5)

But this does not give the correct answer (since step is < 0 and is not accounted for in get_rasterio_transform

coords_src = Coordinates([clinspace(10, 0, 3), clinspace(0, 10, 5)], dims=["lat", "lon"])
get_rasterio_transform(coords_src)
Affine(2.5, 0.0, -1.25,
       0.0, -1.6666666666666667, 7.5)

So there is a small bug there... I'll fix it.

jmilloy · 2020-04-20T17:49:47Z

return rasterio.Affine.from_gdal(*c.geotransform)

I think you can see that I tried that and left it commented it out. It is not the same and doesn't work.

Yeah, I'll need to look into that because it should work. They are doing the same thing, or should be. This means one of the two are wrong.

Yeah, I see the issue. The coordinates.geotransform "naturally" handles the order of the dimensions (is_descending). For example, these give the same results:
coords_src = Coordinates([clinspace(10, 0, 3), clinspace(0, 10, 5)], dims=["lat", "lon"])
rasterio.Affine.from_gdal(*coords_src.geotransform)
Affine(2.5, 0.0, -1.25,
       0.0, -5.0, 12.5)
coords_src = Coordinates([clinspace(0, 10, 3), clinspace(0, 10, 5)], dims=["lat", "lon"])
get_rasterio_transform(coords_src)
Affine(2.5, 0.0, -1.25,
       0.0, -5.0, 12.5)
But this does not give the correct answer (since step is < 0 and is not accounted for in get_rasterio_transform
coords_src = Coordinates([clinspace(10, 0, 3), clinspace(0, 10, 5)], dims=["lat", "lon"])
get_rasterio_transform(coords_src)
Affine(2.5, 0.0, -1.25,
       0.0, -1.6666666666666667, 7.5)
So there is a small bug there... I'll fix it.

Aha! Thanks.

jmilloy · 2020-04-20T17:52:34Z

Well, instead of fixing the bug in get_rasterio_transform, I would prefer to use geotransform here if that works. I feel like when I tried it, the tests didn't pass... maybe even threw an exception, but I could be wrong.

mpu-creare · 2020-04-20T18:08:06Z

Well, instead of fixing the bug in get_rasterio_transform, I would prefer to use geotransform here if that works. I feel like when I tried it, the tests didn't pass... maybe even threw an exception, but I could be wrong.

Yep, that's the plan! I think tests are now passing...

…Also updating unit tests to catch an error in case of descening 'lat' coordinates.

jmilloy added 2 commits April 17, 2020 10:21

Merge branch 'develop' into feature/move-segments

90ab00a

jmilloy self-assigned this Apr 17, 2020

jmilloy added 5 commits April 17, 2020 17:20

WIP: Adds DataSource coordinates boundary trait, Coordinates get_area…

fec45fe

…_bounds method, and uses both in interpolation.

WIP: Clean up references to ctype, segment lengths, and area_bounds i…

54b4e04

…n docs and docstrings.

WIP: Partial support for non-uniform coordinate boundaries. DataSourc…

5da4cb8

…e indexes non-uniform boundaries during evaluation as needed, but interpolation raises an exception for non-uniform or non-centered boundaries.

WIP: Default boundary for uniform coordinates.

70ecc41

Pass datasource boundary through to select in order to use area_bound…

f00936a

…s when necessary.

jmilloy added 2 commits April 19, 2020 11:42

WIP: No longer need to pass the boundary into intersect/select.

36d3418

jmilloy changed the title ~~WIP: move segments from coordinates to datasource~~ Move segments from coordinates to datasource Apr 19, 2020

jmilloy requested a review from mpu-creare April 19, 2020 16:16

WIP: Tests uncommented.

d0e39fe

jmilloy commented Apr 19, 2020

View reviewed changes

podpac/core/coordinates/coordinates.py Show resolved Hide resolved

WIP: cleanup stray TODO.

ccff81b

jmilloy commented Apr 19, 2020

View reviewed changes

podpac/core/coordinates/coordinates1d.py Outdated Show resolved Hide resolved

jmilloy commented Apr 19, 2020

View reviewed changes

podpac/core/coordinates/test/test_array_coordinates1d.py Show resolved Hide resolved

jmilloy commented Apr 19, 2020

View reviewed changes

jmilloy added 2 commits April 19, 2020 12:33

Merge branch 'develop' into feature/move-segments

4f291ce

WIP: Make the Datasource boundary indexing method private.

b72ad85

jmilloy mentioned this pull request Apr 20, 2020

Feature/uniform coords transform to uniform #396

Merged

mpu-creare requested changes Apr 20, 2020

View reviewed changes

BUGFIX: get_boundary tests uptade to use _get_boundary method (which …

89ad856

…was renamed).

FIX: Fixing rasterion get_transform to use more general machinery. …

0ca7e5c

…Also updating unit tests to catch an error in case of descening 'lat' coordinates.

mpu-creare approved these changes Apr 20, 2020

View reviewed changes

mpu-creare merged commit bc7c1d9 into develop Apr 20, 2020

mpu-creare deleted the feature/move-segments branch April 20, 2020 18:31

mpu-creare mentioned this pull request Apr 20, 2020

Move ctype/segments from Coordinates to DataSource #357

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Move segments from coordinates to datasource #395

Move segments from coordinates to datasource #395

jmilloy commented Apr 17, 2020 •

edited

Loading

coveralls commented Apr 17, 2020 •

edited

Loading

mpu-creare commented Apr 19, 2020

jmilloy commented Apr 19, 2020

mpu-creare commented Apr 19, 2020

jmilloy commented Apr 19, 2020

jmilloy commented Apr 19, 2020

jmilloy commented Apr 19, 2020

jmilloy Apr 19, 2020

jmilloy Apr 19, 2020

mpu-creare Apr 20, 2020

jmilloy commented Apr 19, 2020

mpu-creare commented Apr 20, 2020

mpu-creare left a comment

mpu-creare Apr 20, 2020

mpu-creare Apr 20, 2020

mpu-creare Apr 20, 2020

mpu-creare Apr 20, 2020

jmilloy Apr 20, 2020

jmilloy Apr 20, 2020

mpu-creare Apr 20, 2020

mpu-creare Apr 20, 2020

jmilloy commented Apr 20, 2020

jmilloy commented Apr 20, 2020

mpu-creare commented Apr 20, 2020

mpu-creare commented Apr 20, 2020

jmilloy commented Apr 20, 2020

jmilloy commented Apr 20, 2020

mpu-creare commented Apr 20, 2020

		# get indexed boundary
		self._requested_source_boundary = self.get_boundary(self._requested_source_coordinates_index)

Move segments from coordinates to datasource #395

Move segments from coordinates to datasource #395

Conversation

jmilloy commented Apr 17, 2020 • edited Loading

Summary

Datasource coordinates boundary

get_area_bounds

coveralls commented Apr 17, 2020 • edited Loading

mpu-creare commented Apr 19, 2020

jmilloy commented Apr 19, 2020

mpu-creare commented Apr 19, 2020

jmilloy commented Apr 19, 2020

jmilloy commented Apr 19, 2020

jmilloy commented Apr 19, 2020

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jmilloy commented Apr 19, 2020

mpu-creare commented Apr 20, 2020

mpu-creare left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jmilloy commented Apr 20, 2020

jmilloy commented Apr 20, 2020

mpu-creare commented Apr 20, 2020

mpu-creare commented Apr 20, 2020

jmilloy commented Apr 20, 2020

jmilloy commented Apr 20, 2020

mpu-creare commented Apr 20, 2020

jmilloy commented Apr 17, 2020 •

edited

Loading

coveralls commented Apr 17, 2020 •

edited

Loading