busstoptaktik / geodesy

Rust geodesy
Apache License 2.0
66 stars 6 forks source link

Should RG accept ellipse params as well as an ellps=<name>? #61

Closed Rennzie closed 10 months ago

Rennzie commented 10 months ago

While testing various output from projinfo I've come across an example where the ellipse is described in terms of rf and a parameters instead of ellps=<name>. For example a transformation between EPSG:2001 and EPSG:3857 looks like:

+proj=pipeline
  +step +inv +proj=tmerc +lat_0=0 +lon_0=-62 +k=0.9995 +x_0=400000 +y_0=0 +a=6378249.145 +rf=293.465
  +step +inv +proj=longlat +a=6378249.145 +rf=293.465 # Note: I've implemented longlat as a noop in geodesy-wasm.
  +step +proj=push +v_3
  +step +proj=cart +a=6378249.145 +rf=293.465
  +step +proj=helmert +x=-255 +y=-15 +z=71
  +step +inv +proj=cart +ellps=WGS84
  +step +proj=pop +v_3
  +step +proj=webmerc +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84

# With the equivalent RG pipeline being:

  | tmerc inv lat_0=0 lon_0=-62 k=0.9995 x_0=400000 y_0=0  a=6378249.145 rf=293.465
  | noop inv a=6378249.145 rf=293.465
  | push v_3
  | cart a=6378249.145 rf=293.465
  | helmert x=-255 y=-15 z=71
  | cart inv ellps=WGS84
  | pop v_3
  | webmerc lat_0=0 lon_0=0 x_0=0 y_0=0 ellps=WGS84

When converting the same coordinate with PROJ and RG I get differences in the order of 10's of meters. The a=6378249.145 rf=293.465 parameters describe the clrk80 ellipse, so if I update both pipelines to use ellps=clrk80 instead of using a and rf the results are the same.

The question is, should RG be able to handle ellipse input via the rf and a parameters? If yes, I'm happy to attempt a PR with some guidance on where to do this and what else needs to be covered.

busstoptaktik commented 10 months ago

I decided skipping the "anything allowed" approach of PROJ in favor of the (I believe) "only a and rf" dogma of EPSG, since "if it does not appear in practice, no reason to support it" - or at least that was my bad excuse. Really I just had to shoehorn and cut down to be able to get to a 80% solution in 0.1% of the time.

The workaround is to use the "a,rf" syntax for ellps, i.e. ellps=6378249.145,293.46 in your case. This is handled by Ellipsoid::named().

The all-encompassing solution would probably either be rather intrusive, probably starting somewhere like here, or be based on additional pattern recognition in the parse_proj() function, scanning for all valid PROJ ellipsoid definition tokens, and translating them into a ellps=a,rf pair.

Note, however, that PROJ's interpretation of redundant ellipsoid specifiers may be surprising - e.g. the combination ellps=GRS80 a=1 results in the instantiation of an ellipsoid the shape of GRS80 but with a size of unity. If you decide to give it a try, also consider doing a bit of work documenting known (and potentially unknown) differences: parse_proj is horribly underdocumented (totally my fault)

Rennzie commented 10 months ago

Ok thanks for the guidance @busstoptaktik. I'll have a think and give it a whirl as time permits.