[PROJ] [EXTERNAL] Re: Proj 6 API questions

[Kristian Evers]

In response to Jochem, Nyall wrote:

> > > To add a 3rd opinion in the debate: I think the ellipsoid should should
> > > NOT automatically be set to the ellipsoid of a CRS used by the user, but
> > > it should by default kept set to GRS80, unless the user explicitly
> > > manually selects another ellipsoid. GRS80 is scientifically the best
> > > ellipsoid, all others are outdated (not the most accurate estimation of
> > > the size of the earth), unconventional (seldomly used) or wrong (e.g. due
> > > to rounding of parameters). The only exception would be extra-terrestrial
> > > CRSs: there GRS80 make no sense and the scientifically best ellipsoid of
> > > the planet in question should be selected by default. When a user is
> > > using a national CRS based on a old ellipsoid, (s)he mostly doesn't want
> > > distances and areas based on this outdated ellispoid but distances and
> > > areas computed in the best way possible, which is GRS80.

> > I'm very much in favour of this approach!

I can understand that you are in favour of this approach; It is easy and simple.
Unfortunately, geodesy is rarely easy and simple. If you want your calculations
to be correct you need to accept this fact. Even though the solution is not
always straight forward. To explain this we need to brush up our understanding
of what a coordinate reference system is. According to the EPSG guidance note
7-1 a CRS is comprised of a datum and a coordinate system. The datum is where
we define the ellipsoid (and often times other parameters) and the coordinate
system is broadly speaking the projection. This division works really well with
modern data where we can easily separate the two because we have access
to global positioning systems. For legacy coordinate reference systems
the datum and coordinate system are often times very tightly tied together
because physical measurements of coordinates rely on both the datum and
the coordinate system. This has many unfortunate effects, for instance that
the scale factor of the used projection is affecting distances measures etc.

The whole point of back-projecting coordinates to the ellipsoid is to be able
to determine distances without being affected by distortions inherent to the
used projection. If we back-project linear coordinates to a different ellipsoid
than the one a CRS is based on we inherently introduce scale differences
again. So this should be avoided if the aim is to provide as exact distance
and area measurements as possible.

On the surface this can seem like a problem from the past that we in modern
days do not have to deal with, but that is unfortunately not the case. There
are many examples of modern day use of legacy CRS's, for example cadastral
work in many countries and in territorial disputes between national states
where the original data is usually only available in older CRS's. Both scenarios
are involve lawyers and incorrect area calculations can end up having severe
consequences for those involved. That is my reason for pointing this out - It
would be nice if my GIS tools of choice dealt with this correctly so I don't have
to worry about it.

So, to handle this correctly I would suggest a two-step approach:

1. back-project the projected coordinates to geodetic coordinates using the
    Ellipsoid parameters of the CRS datum.
2. Perform geodesic calculations using the ellipsoid parameters of the CRS datum

This is best done on a per dataset/layer basis, using the CRS that the data is
registered as in the metadata of the file.

If I undertand the documentation correctly I believe that using a PJ
returned from proj_crs_get_coordoperation in the inverse direction should be
able to transform projected coordinates to geodetic coordinates using the
proper ellipsoid.

/Kristian


-----Oprindelig meddelelse-----
Fra: PROJ <proj-bounces at lists.osgeo.org> På vegne af Nyall Dawson
Sendt: 27. marts 2019 02:28
Til: Even Rouault <even.rouault at spatialys.com>
Cc: PROJ <proj at lists.osgeo.org>
Emne: Re: [PROJ] [EXTERNAL] Re: Proj 6 API questions

On Tue, 26 Mar 2019 at 20:57, Even Rouault <even.rouault at spatialys.com> wrote:
>
> On mardi 26 mars 2019 19:56:22 CET Nyall Dawson wrote:
> > On Tue, 26 Mar 2019 at 17:58, Lesparre, Jochem
> >
> > <Jochem.Lesparre at kadaster.nl> wrote:
> > > > I honestly don't have any plans to change this behavior. I cannot see
> > > > a way to remove this choice from users, as QGIS projects don't have a
> > > > single CRS, but they DO need a single ellipsoid. The best we can do is
> > > > make sure that a good default selection is given to users, but we
> > > > cannot totally remove this choice.
> > >
>
> The question raised by Kristian is still valid. How do you convert from
> project coordinates to coordinates on the measurement ellipsoid ?
> - is is a no-op conversion, that is
>     lon,lat(measurement_ellipsoid) = lon,lat(project_crs)  ?
>   That would be wrong.
> - or some form of coordinate transformation from the project CRS to the
>   measurement ellipsoid, but then you'll get no datum shifts. That said, for
>   measurement, that probably doesn't make a lot of difference in the resulting
>   computation (assuming that the shift would be almost homogeneous on the area
>   of interest)

See earlier in this thread -- it's the second, we reproject from the
layer coordinates to a "+proj=longlat +a=%1 +b=%2 +no_defs" or
"+proj=longlat +a=%1 +rf=%2 +no_defs" definition.

>
> > I wonder if this is something which belongs in the proj db itself.
> > I.e. some API to query the recommended ellipsoid for a given planetary
> > body.
>
> Before the API, we would need the information. The database doesn't store a
> concept of recommended ellipsoid. Could be a new column in 'ellipsoid' table
> with a check that there's a single recommended ellipsoid for a celestial_body,
> and that for each body, there's a recommended ellipsoid (commit.sql can be
> used to add consistency checks that don't fall in usual SQL constraints)
> (thinking out loud. haven't given that more thought)

That's exactly what I had in mind.

What's semi_major_axis in the celestial_body table used for? To me
this is basically a default ellipsoid definition for those bodies, but
bypassing the ellipsoid table. Maybe it makes sense to remove this
column?

Nyall


>
> Even
> --
> Spatialys - Geospatial professional services
> http://www.spatialys.com
_______________________________________________
PROJ mailing list
PROJ at lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/proj