fatiando
Coordinate conversions proposal
I propose that we introduce a new set of unified coordinate conversion routines for translating between geodetic, spherical, ellipsoidal-harmonic and cartesian coordinates. These would partially replace those that exist, and would introduce a more uniform logic for Boule. All coordinate transformations would require knowledge of the ellipsoid parameters and are hence appropriate for this project.
geodetic, spherical, and ellipsoidal-harmonic conversions
These routines would convert between "latitude" and "height" above the ellipsoid. Latitude would be geodetic, spherical or reduced (for ellipsoidal-harmonic coordinates). Height would be geodetic height (ellipsoidal height) or spherical height (height above the ellipsoid along the radius vector direction). Ellipsoidal-harmonic coordinate don't use a height, but instead use "u", which is the ellipsoid semiminor axis that passes through the point.
None of these routines require longitude, as the longitude is the same for all coordinate systems:
spherical_latitude, spherical_height = geodetic_to_spherical(latitude, height)
latitude, height = spherical_to_geodetic(spherical_latitude, spherical_height)
reduced_latitude, u = geodetic_to_ellipsoidal_harmonic(latitude, height)
latitude, height = ellipsoidal_harmonic_to_geodetic(reduced_latitude, u)
reduced_latitude, u = spherical_to_ellipsoidal_harmonic(spherical_latitude, spherical_height)
spherical_latitude, spherical_height = ellipsoidal_harmonic_to_spherical(reduced_latitude, u)
One could potentially simplify these routines by considering the height as an optional parameter with a default of zero, and the parameter u as an optional parameter with a default of semiminor_axis.
Cartesian conversions
Conversion to and from cartesian coordinates require knowledge of the longitude.
x, y, z = geodetic_to_cartesian(longitude, latitude, height)
x, y, z = spherical_to_cartesian(longitude, spherical_latitude, spherical_height)
x, y, z = ellipsoidal_harmonic_to_cartesian(longitude, reduced_latitude, u)
longitude, latitude, height = cartesian_to_geodetic(x, y, z)
longitude, spherical_latitude, spherical_height = cartesian_to_spherical(x, y, z)
longtidue, reduced_latitude, u = cartesian_to_ellipsoidal_harmonic(x, y, z)
To simplify these functions, one could consider the height and u parameters as optional.
Comments
- There was a discussion earlier about deprecating some coordinate transforms and moving these to pymap3d (https://github.com/fatiando/boule/pull/126). What is being proposed here is somewhat different: In that discussion, the spherical coordinates were latitude and radius, which didn't depend on knowledge of the ellipsoid (see als #165). Here they are latitude and spherical height, which do require knowledge of the ellipsoid
- At the present time, there is no easy way to convert from geodetic to spherical latitude, which is probably the most commonly used conversion. This was because you needed to input the latitude and radius, and radius needed to be computed by another function. What is proposed here rectifies this problem.
- Even if some of these functions could be implemented in pymap3d instead of here, pymap3d does not have all of the ellipsoids that we have. This would complicate any such conversion. Furthermore, as far as i can tell, pymap3d doesn't really implement the above routines when height above the ellipsoid is non-zero.
- This would be a breaking change.
Are you willing to help implement and maintain this feature? Yes, but I am not sure when I will get around to it. Help (even for just 1 of the above) would be appreciated!
@MarkWieczorek I've changed my mind a bit about minimizing overlap with pymap3d. I guess it's not a problem to repeat a few things for ease of use and consistency of our API. Sorry about pushing against this earlier.
Regarding the proposal: the main think I don't understand is why use a spherical_height instead of a radius? It's not really a coordinate that comes up in most equations and would be a pain to determine. I'd propose to just use the radius.
A few things I don't quite like but I'd like to know what this looks like from your experience:
Omitting longitude: I know it's not needed but when dealing with data we often already have it. For example:
point = (0, 14, 1500) # lon, lat, height
point_sph = WGS84.geodetic_to_spherical(*point)
otherwise this would look like:
point = (0, 14, 1500) # lon, lat, height
lat_sph, radius = WGS84.geodetic_to_spherical(*point[1:])
point_sph = (point[0], lat_sph, radius)
At least I'm not usually carrying coordinates in individual arrays since it's a pain to have to name them all. But this is how I code these things. I'm curious to know if you usually do things differently and passing longitude is a chore.
Defaulting to height/u of zero: I'm not a fan of having the default coordinate like this since it feels arbitrary. Kind of like default latitude to the equator. It's easy enough to pass in 0 if needed and it makes the calculation more explicit.
At the present time, there is no easy way to convert from geodetic to spherical latitude, which is probably the most commonly used conversion. This was because you needed to input the latitude and radius, and radius needed to be computed by another function. What is proposed here rectifies this problem.
This I also don't understand. There is no way to convert latitudes without specifying a radius or height since the conversion depends on those parameters. It would only make sense if the you mean converting for points on the surface of the ellipsoid. Is that it?
Regarding the proposal: the main think I don't understand is why use a spherical_height instead of a radius? It's not really a coordinate that comes up in most equations and would be a pain to determine. I'd propose to just use the radius.
Either way would be fine for me. Specifying a spherical height would make this analogous to specifying a geodetic height. That's the only reason I suggested this.
At least I'm not usually carrying coordinates in individual arrays since it's a pain to have to name them all. But this is how I code these things. I'm curious to know if you usually do things differently and passing longitude is a chore.
I do it both ways. For me its common to just have individual arrays of latitude and longitude. One option would be to allow to input either a coordinate (array of three values) or to input individual arrays for latitude, longitude, height, etc. For most of my work, the height is always zero and longitude isn't important. This might be the most "useful" way of doing this, but then we would probably have to have some complicated check that the input keywords are consistent.
One more thing: The reason that I use individual arrays of latitude and longitude, is that I often work with global grids of a specified size, and the lat and lon values are just the numerical values of the two coordinates.