judgej / Proj4

Experiments in a Proj4 reboot
MIT License
3 stars 0 forks source link

Geodetic conversion #11

Closed SamMousa closed 8 years ago

SamMousa commented 8 years ago

According to most documentation, a geodetic projection is just a transformation that consists of 3 operations:

  1. Translation
  2. Rotation
  3. Scaling

These can all be captured in 7 parameters used in the transformation matrix. https://en.wikipedia.org/wiki/Transformation_matrix#Affine_transformations

When converting coordinates from one geodetic coordinate system to another we can first construct a transformation matrix and then apply that to the input coordinates.

Assuming the above is correct:

  1. There should be a class that represents a transformation matrix.
  2. There should be functionality to combine these transformation matrices. (I think that is what the forward / inverse functions are (partially) doing right now.

Advantage of doing it this way is you can increase precision and improve performance by doing less calculations.

This page lists several methods to convert from geodetic to geodetic without going through geocentric. https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#Datum_transformations

judgej commented 8 years ago

I would like to keep the ability to shift the datum in geocentric space, even if shifting in geodetic space is supported too. I believe the former was chosen because it was more accurate and faster for computing, but that may have been the case thirty years ago - I'm not sure how they compare for speed and accuracy in today's PHP.

Fewer calculations do not necessarily mean more accuracy or precision. It really depends on the intermediate steps and what precision each of those must have for a given final precision. In this case, I don't know the exact details, but these things are never simple.

I'd be happy with both, then they can be used as needed depending on circumstances.

judgej commented 8 years ago

That datum transformation diagram is great. We should be able to support each of those paths. But unlike Proj4, keep the datum details in the points, so they always retain context.

SamMousa commented 8 years ago

Okay, so when taking into account the fact that we don't want to lose source information (since that is assumed to be most precise).

I think we should have:

  1. Abstract coordinate class that uses X,Y, Z
  2. Abstract coordinate class that uses lat, long, h
  3. Point that is a tuple of Coordinate and Datum
  4. Concrete classes GeodeticCoordinate and GeocentricCoordinate both implementing CoordinateInterface
  5. Interface CoordinateInterface which has functions getX(), getY(), getZ(), getLon(), getLat(), getH(); returning either calculated or source values (depending on the type of the source data).

Question: to convert from geodetic to geocentric do we need to know the datum? If I give you X, Y, Z can you give me Lat, Lon, H for the same (but unknown to you) datum?

Then we add:

  1. MapperInterface, which exposes function CoordinateInterface map(CoordinateInterface $c) which maps a coordinate from a source datum. This is merely the exposed API, a simple implemention would be an IdentityMapper which wouldn't need any parameters. Or an ABMapper which might take 2 datums in its constructor.
  2. Optionally we could add a function mapPoint() to this interface and implement it in a trait, this would then return a Point with the datum set to the target datum. If the point is not in the same source datum as the mapper then an exception is thrown. (Basically from an implementors point of view you don't care about points but about coordinates, from a user point of view you care about mapping points and having data that has context).

The concrete mappers could:

  1. Use a fixed route ([Geodetic] -> ECEF -> Helmert -> ECEF).
  2. The mapper should always do as little as possible, if it converts a geodetic coordinate by converting it to ECEF and then applying Helmert, it should return a GeocentricCoordinate instance, since the conversion to geodetic will be done when it is needed in the interface.

So a mapper might look like this (pseudo code):

public function map(CoordinateInterface $c) {
    // We support Molodensky
    if ($c instanceof GeodeticCoordinate) {
        return $this->molodensky($c);
    } else {
        // Source is geocentric.
       return $this->helmert($c);
    }
}

Or like this, in case it doesn't support Molodensky:

public function map(CoordinateInterface $c) {
    // We don't support Molodensky
    // Even though we might have a geodetic coordinate as input; the concrete class must implement the `getX()`, `getY()` and `getZ()`.
    // We return a geocentric coordinate (because that is what helmert gives us and we don't want to do extra conversions. (And if you need the geodetic coordinate values you can them via the exact same interface (`getLan()`, `getLon()` and `getH()`)
   return $this->helmert($c);
}
judgej commented 8 years ago

Question: to convert from geodetic to geocentric do we need to know the datum? If I give you X, Y, Z can you give me Lat, Lon, H for the same (but unknown to you) datum?

Yes, the same datum carries forward.

SamMousa commented 8 years ago

And there is always just these 2 pairs of 3 attributes (whcih may appear in different units) or are there more coordinate types?

SamMousa commented 8 years ago

https://github.com/SamMousa/Proj4/commit/7c7de77390b94d51d090e4545c1b9dd1082a5665

Have implemented the basic structure, any comments?

SamMousa commented 8 years ago

Yes, the same datum carries forward.

But can you convert without knowing the datum? As far as I can tell from the current source code you need to know the radius of the ellipsoid, meaning you cannot give me x y z when given just lat lon h.

SamMousa commented 8 years ago

Closing this as I've implemented a generic coordinate that supports both geodetic and geocentric values given an ellipsoid.