Repository: fatiando/verde
Pull request #231 (2020-02-26 18:33:28)

Fix bug with profile distances being unprojected
`BaseGridder.profile` has the option to apply a projection to the
coordinates before predicting so we can pass geographic coordinates to
cartesian gridders. In these cases, the distance along the profile is
calculated by `profile_coordinates` with the unprojected coordinates (in
the geographic case it would be degrees). Returning unprojected
distances is not very useful and instead we should have been returning
projected distances. They need to be calculated from the projected
coordinates (out of `profile_coordinates`) because the points are evenly
spaced in unprojected coordinates.





**Reminders**:

- [ ] Run `make format` and `make check` to make sure the code follows the style guide.
- [ ] Add tests for new features or tests that would have caught the bug that you're fixing.
- [ ] Add new public functions/methods/classes to `doc/api/index.rst` and the base `__init__.py` file for the package.
- [ ] Write detailed docstrings for all functions/classes/methods. It often helps to design better code if you write the docstrings first.
- [ ] If adding new functionality, add an example to the docstring, gallery, and/or tutorials.
- [ ] Add your full name, affiliation, and [ORCID](https://orcid.org) (optional) to the `AUTHORS.md` file (if you haven't already) in case you'd like to be listed as an author on the [Zenodo](https://zenodo.org/communities/fatiando) archive of the next release.


leouieda: @santisoler @jessepisel I could use another set of eyes on this one. Using the example of providing geographic points to `profile` and a projection:

1. Generate profile points in degrees before projecting (so they are evenly spaced in degrees).
2. Project the coordinates and predict.
3. Re-calculate distances between the projected coordinates so we don't return "distances" in degrees.

Does this seem reasonable? I found this bug while extracting a profile for modelling in 2D Cartesian coordinates and getting distances in degrees (which are kind of useless). Then again, maybe I should have just projected the coordinates of the points before passing them to `profile`. What do you think?
leouieda: Actually, it might make more sense to project the points before creating the profile. The way we were generating the profile coordinates is probably wrong if given geographic coordinates. 
santisoler: IMHO, I think we should avoid projecting inside the method and ask already projected coordinates.
This would make the workflow more transparent and it will increase users' code just by one line (projecting the coordinates, the projection has to be defined on both cases).
Worth remembering that projections may introduce errors if the projected region is too big, so I think we should avoid introducing such source of error inside the `profile` method.

If we would like to create profiles on geographic coordinates maybe we could write a different method for that, which don't involve projecting coordinates but computing equidistant points located on a sphere. But that's for another PR of course.
santisoler: It worth noting that this workflow would work if the data region is not too big. Otherwise we will be introducing errors due to projection. A solution for this would be to create profiles directly on geographic (or spherical) coordinates, where equidistant points can be defined as point located at constant arc-length distances. But this should be tackled on a different PR.
ahartikainen: What library are you using to calculate distances? Would it be a good idea to calculate distances with `pyproj` where it calculates the geographic distance?
leouieda: @santisoler yeah I agree that this is not the best workflow. The main problem is that the gridder isn't aware of coordinate systems. They are all Cartesian at the moment as well so the problem with large areas and profiles is present for all our gridding. 

@ahartikainen we're not really using any libraries for the reasons outlines above. Right now, we just assume that everything is Cartesian. We could use pyproj for great circle distances but that might involve having `profile` be aware of Cartesian vs geographic (which it currently isn't).

It seems more and more to be that we can't use the same abstraction for `profile` and `grid` (pass in the projection and don't ask questions). We might have to add an option to inform profile that we're talking about geographic coordinates.

The ideal case for geographic profiles would be:

1. Make the profile evenly spaced in actual physical coordinates (not degrees).
2. Return distances as great circle distances instead of "degrees".

When extracting profiles, we're usually plotting value x distance so it's more useful to have things work this way most times.

leouieda: A few thoughts after discussion with @santisoler:

* The Verde gridders are currently all Cartesian by nature.
* Even if they take a `projection` argument, the actual prediction is still happening in Cartesian coordinates.
* So it makes little sense to calculate distances on the sphere and then flatten it out for interpolation.
* In case we have spherical gridders in the future (#173), then it would make sense for those gridders to calculate distances on the sphere by default. That would require changes to how we generate the profile coordinates and it will likely have to be different than for the Cartesian case. We can deal with that when the need arises. 

I think the most sensible approach right now is:

* Project the profile points before passing them to `profile_coordinates`. That way, the profile coordinates are generated in Cartesian coordinates and so are the distances.
* Do the inverse projection of the profile coordinates before sending them back to the user. That way the profile location can be plotted on a map.

This would make the profile point generation be Cartesian (as it should for the current gridders) and solve the problem of returning "distances" as angles.

@ahartikainen @santisoler what do you think?
leouieda: I implemented projecting the input points and then projecting back the Cartesian coordinates. This way, the distances returned are Cartesian (projected) and evenly spaced in the projection domain. I also added a section in the projections tutorial about profiles:

![2020-03-17-111057_3200x1800_scrot](https://user-images.githubusercontent.com/290082/76851566-6121f600-6841-11ea-93d3-640935debed0.jpg)

My main question right now is that this breaks compatibility slightly since the returned points and distances are going to be different from before. But I consider the old behaviour to be a bug so it might not be too bad.

If we consider this to break the API compatibility, I'd be in favour of holding it until Verde 2.0. If not, then we can merge as include it in 1.4.0, which I want to release once #237 is done.

What do you think?

cc @santisoler @ahartikainen @jessepisel 

santisoler: This is looking great @leouieda!

I agree to merge this. Although it somehow breaks backward compatibility, I agree the previous implementation was bugged, so I think we can prioritize accurate behavior over backward compatibility. Maybe we could post a comment about this on Slack after we merge this so most users can become aware of this. In fact we could ask their opinion before merging.


leouieda: > Although it somehow breaks backward compatibility, I agree the previous implementation was bugged, so I think we can prioritize accurate behavior over backward compatibility. 

Yep, and also change the values of coordinates but not the number/types of arguments returned or any API so we should be fine. This is also not a very well publicized method and likely not very widely used.

Another option is to include a warning about the change which we can remove in 2.0.
leouieda: Added more details to the docstring and a warning (docs only) about the change.

Patch:
diff --git a/tutorials/projections.py b/tutorials/projections.py
index 16edce9..0d6f34b 100644
--- a/tutorials/projections.py
+++ b/tutorials/projections.py
@@ -132,5 +132,66 @@ pc = grid_geo.bathymetry.plot.pcolormesh(
 )
 plt.colorbar(pc).set_label("meters")
 vd.datasets.setup_baja_bathymetry_map(ax, land=None)
+plt.show()
+
+########################################################################################
+# Profiles
+# --------
+#
+# For profiles, things are a bit different. The projection is applied to the
+# input points before coordinates are generated. So the profile will be evenly
+# spaced in *projected coordinates*, not geographic coordinates. This is to
+# avoid issues with calculating distances on a sphere.
+#
+# The coordinates returned by the ``profile`` method will be in geographic
+# coordinates, so projections given to ``profile`` must take an ``inverse``
+# argument so we can undo the projection.
+#
+# The main utility of using a projection with ``profile`` is being able to pass
+# in points in geographic coordinates and get coordinates back in that same
+# system (making it easier to plot on a map).
+#
+# To generate a profile cutting across our bathymetry data, we can use
+# longitude and latitude points taken from the map above).
+
+start = (-114.5, 24.7)
+end = (-110, 20.5)
+profile = spline.profile(
+    point1=start,
+    point2=end,
+    size=200,
+    projection=projection,
+    dims=("latitude", "longitude"),
+    data_names=["bathymetry"],
+)
+print(profile)
+
+########################################################################################
+# Plot the profile location on our geographic grid from above.
+
+plt.figure(figsize=(7, 6))
+ax = plt.axes(projection=ccrs.Mercator())
+ax.set_title("Profile location")
+pc = grid_geo.bathymetry.plot.pcolormesh(
+    ax=ax, transform=ccrs.PlateCarree(), vmax=0, zorder=-1, add_colorbar=False
+)
+plt.colorbar(pc).set_label("meters")
+ax.plot(profile.longitude, profile.latitude, "-k", transform=ccrs.PlateCarree())
+ax.text(start[0], start[1], "A", transform=ccrs.PlateCarree())
+ax.text(end[0], end[1], "B", transform=ccrs.PlateCarree())
+vd.datasets.setup_baja_bathymetry_map(ax, land=None)
+plt.show()
+
+########################################################################################
+# And finally plot the profile.
+
+plt.figure(figsize=(8, 3))
+ax = plt.axes()
+ax.set_title("Profile of bathymetry (A-B)")
+ax.plot(profile.distance, profile.bathymetry, "-k")
+ax.set_xlabel("Distance (m)")
+ax.set_ylabel("Bathymetry (m)")
+ax.set_xlim(profile.distance.min(), profile.distance.max())
+ax.grid()
 plt.tight_layout()
 plt.show()
diff --git a/verde/base/base_classes.py b/verde/base/base_classes.py
index bc581c8..cc988d0 100644
--- a/verde/base/base_classes.py
+++ b/verde/base/base_classes.py
@@ -409,6 +409,30 @@ class BaseGridder(BaseEstimator):
         Includes the calculated Cartesian distance from *point1* for each data
         point in the profile.
 
+        To specify *point1* and *point2* in a coordinate system that would
+        require projection to Cartesian (geographic longitude and latitude, for
+        example), use the ``projection`` argument. With this option, the input
+        points will be projected using the given projection function prior to
+        computations. The generated Cartesian profile coordinates will be
+        projected back to the original coordinate system. **Note that the
+        profile points are evenly spaced in projected coordinates, not the
+        original system (e.g., geographic)**.
+
+        .. warning::
+
+            **The profile calculation method with a projection has changed in
+            Verde 1.4.0**. Previous versions generated coordinates (assuming
+            they were Cartesian) and projected them afterwards. This led to
+            "distances" being incorrectly handled and returned in unprojected
+            coordinates. For example, if ``projection`` is from geographic to
+            Mercator, the distances would be "angles" (incorrectly calculated
+            as if they were Cartesian). After 1.4.0, *point1* and *point2* are
+            projected prior to generating coordinates for the profile,
+            guaranteeing that distances are properly handled in a Cartesian
+            system. **With this change, the profile points are now evenly
+            spaced in projected coordinates and the distances are returned in
+            projected coordinates as well.**
+
         Parameters
         ----------
         point1 : tuple
@@ -433,13 +457,16 @@ class BaseGridder(BaseEstimator):
             ``['east_component', 'north_component', 'vertical_component']`` for
             3D vector data.
         projection : callable or None
-            If not None, then should be a callable object
-            ``projection(easting, northing) -> (proj_easting, proj_northing)``
-            that takes in easting and northing coordinate arrays and returns
-            projected northing and easting coordinate arrays. This function
-            will be used to project the generated profile coordinates before
-            passing them into ``predict``. For example, you can use this to
-            generate a geographic profile from a Cartesian gridder.
+            If not None, then should be a callable object ``projection(easting,
+            northing, inverse=False) -> (proj_easting, proj_northing)`` that
+            takes in easting and northing coordinate arrays and returns
+            projected northing and easting coordinate arrays. Should also take
+            an optional keyword argument ``inverse`` (default to False) that if
+            True will calculate the inverse transform instead. This function
+            will be used to project the profile end points before generating
+            coordinates and passing them into ``predict``. It will also be used
+            to undo the projection of the coordinates before returning the
+            results.
 
         Returns
         -------
@@ -448,11 +475,19 @@ class BaseGridder(BaseEstimator):
 
         """
         dims = self._get_dims(dims)
+        # Project the input points to generate the profile in Cartesian
+        # coordinates (the distance calculation doesn't make sense in
+        # geographic coordinates since we don't do actual distances on a
+        # sphere).
+        if projection is not None:
+            point1 = projection(*point1)
+            point2 = projection(*point2)
         coordinates, distances = profile_coordinates(point1, point2, size, **kwargs)
-        if projection is None:
-            data = check_data(self.predict(coordinates))
-        else:
-            data = check_data(self.predict(projection(*coordinates)))
+        data = check_data(self.predict(coordinates))
+        # Project the coordinates back to have geographic coordinates for the
+        # profile but Cartesian distances.
+        if projection is not None:
+            coordinates = projection(*coordinates, inverse=True)
         data_names = get_data_names(data, data_names)
         columns = [
             (dims[0], coordinates[1]),


Test patch:
diff --git a/verde/tests/test_base.py b/verde/tests/test_base.py
index c95eef9..4f2f61d 100644
--- a/verde/tests/test_base.py
+++ b/verde/tests/test_base.py
@@ -128,6 +128,7 @@ def test_basegridder():
     coordinates_true = grid_coordinates(region, shape)
     data_true = angular * coordinates_true[0] + linear
     grid = grd.grid(region, shape)
+    prof = grd.profile((0, -10), (10, -10), 30)
 
     npt.assert_allclose(grd.coefs_, [linear, angular])
     npt.assert_allclose(grid.scalars.values, data_true)
@@ -135,40 +136,68 @@ def test_basegridder():
     npt.assert_allclose(grid.northing.values, coordinates_true[1][:, 0])
     npt.assert_allclose(grd.scatter(region, 1000, random_state=0).scalars, data)
     npt.assert_allclose(
-        grd.profile((0, 0), (10, 0), 30).scalars,
-        angular * coordinates_true[0][0, :] + linear,
+        prof.scalars, angular * coordinates_true[0][0, :] + linear,
     )
+    npt.assert_allclose(prof.easting, coordinates_true[0][0, :])
+    npt.assert_allclose(prof.northing, coordinates_true[1][0, :])
+    npt.assert_allclose(prof.distance, coordinates_true[0][0, :])
 
 
 def test_basegridder_projection():
     "Test basic functionality of BaseGridder when passing in a projection"
 
+    # Lets say the projection is doubling the coordinates
+    def proj(lon, lat, inverse=False):
+        "Project from the new coordinates to the original"
+        if inverse:
+            return (lon / 2, lat / 2)
+        return (lon * 2, lat * 2)
+
+    # Values in "geographic" coordinates
     region = (0, 10, -10, -5)
-    shape = (50, 30)
+    shape = (51, 31)
     angular, linear = 2, 100
     coordinates = scatter_points(region, 1000, random_state=0)
     data = angular * coordinates[0] + linear
+    # Project before passing to our Cartesian gridder
+    grd = PolyGridder().fit(proj(coordinates[0], coordinates[1]), data)
+
+    # Check the estimated coefficients
+    # The grid is estimated in projected coordinates (which are twice as large)
+    # so the rate of change (angular) should be half to get the same values.
+    npt.assert_allclose(grd.coefs_, [linear, angular / 2])
+
+    # The actual values for a grid
     coordinates_true = grid_coordinates(region, shape)
     data_true = angular * coordinates_true[0] + linear
-    grd = PolyGridder().fit(coordinates, data)
 
-    # Lets say we want to specify the region for a grid using a coordinate
-    # system that is lon/2, lat/2.
-    def proj(lon, lat):
-        "Project from the new coordinates to the original"
-        return (lon * 2, lat * 2)
-
-    proj_region = [i / 2 for i in region]
-    grid = grd.grid(proj_region, shape, projection=proj)
-    scat = grd.scatter(proj_region, 1000, random_state=0, projection=proj)
-    prof = grd.profile((0, 0), (5, 0), 30, projection=proj)
+    # Check the scatter
+    scat = grd.scatter(region, 1000, random_state=0, projection=proj)
+    npt.assert_allclose(scat.scalars, data)
+    npt.assert_allclose(scat.easting, coordinates[0])
+    npt.assert_allclose(scat.northing, coordinates[1])
 
-    npt.assert_allclose(grd.coefs_, [linear, angular])
+    # Check the grid
+    grid = grd.grid(region, shape, projection=proj)
     npt.assert_allclose(grid.scalars.values, data_true)
-    npt.assert_allclose(grid.easting.values, coordinates_true[0][0, :] / 2)
-    npt.assert_allclose(grid.northing.values, coordinates_true[1][:, 0] / 2)
-    npt.assert_allclose(scat.scalars, data)
-    npt.assert_allclose(prof.scalars, angular * coordinates_true[0][0, :] + linear)
+    npt.assert_allclose(grid.easting.values, coordinates_true[0][0, :])
+    npt.assert_allclose(grid.northing.values, coordinates_true[1][:, 0])
+
+    # Check the profile
+    prof = grd.profile(
+        (region[0], region[-1]), (region[1], region[-1]), shape[1], projection=proj
+    )
+    npt.assert_allclose(prof.scalars, data_true[-1, :])
+    # Coordinates should still be evenly spaced since the projection is a
+    # multiplication.
+    npt.assert_allclose(prof.easting, coordinates_true[0][0, :])
+    npt.assert_allclose(prof.northing, coordinates_true[1][-1, :])
+    # Distance should still be in the projected coordinates. If the projection
+    # is from geographic, we shouldn't be returning distances in degrees but in
+    # projected meters. The distances will be evenly spaced in unprojected
+    # coordinates.
+    distance_true = np.linspace(region[0] * 2, region[1] * 2, shape[1])
+    npt.assert_allclose(prof.distance, distance_true)
 
 
 def test_check_fit_input():