NetTopologySuite / ProjNet4GeoAPI

.NET Spatial Reference and Projection Engine
GNU Lesser General Public License v2.1
275 stars 83 forks source link

Add EPSG 4979 Support? #79

Open jhaverkost opened 4 years ago

jhaverkost commented 4 years ago

EPSG 4979 is the 3D equivalent of EPSG 4326 with ellipsoidal height as the third dimension.

Attempting to create this coordinate system from the WKT results in an ArgumentException.

Code:

CoordinateSystemFactory csFact = new CoordinateSystemFactory();
coordinateSystem = csFact.CreateFromWkt("GEOGCS[\"WGS 84\",DATUM[\"World Geodetic System 1984\",SPHEROID[\"WGS 84\",6378137.0,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0.0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.017453292519943295],AXIS[\"Geodetic latitude\",NORTH],AXIS[\"Geodetic longitude\",EAST],AXIS[\"Ellipsoidal height\",UP],AUTHORITY[\"EPSG\",\"4979\"]]");

Error:

ArgumentException: Axis info should contain two axes for horizontal coordinate systems
ProjNet.CoordinateSystems.HorizontalCoordinateSystem..ctor (ProjNet.CoordinateSystems.HorizontalDatum datum, System.Collections.Generic.List`1[T] axisInfo, System.String name, System.String authority, System.Int64 code, System.String alias, System.String remarks, System.String abbreviation) (at <45d4f050441b42ac9796be647bb4c7ef>:0)
ProjNet.CoordinateSystems.GeographicCoordinateSystem..ctor (ProjNet.CoordinateSystems.AngularUnit angularUnit, ProjNet.CoordinateSystems.HorizontalDatum horizontalDatum, ProjNet.CoordinateSystems.PrimeMeridian primeMeridian, System.Collections.Generic.List`1[T] axisInfo, System.String name, System.String authority, System.Int64 authorityCode, System.String alias, System.String abbreviation, System.String remarks) (at <45d4f050441b42ac9796be647bb4c7ef>:0)
ProjNet.IO.CoordinateSystems.CoordinateSystemWktReader.ReadGeographicCoordinateSystem (ProjNet.IO.CoordinateSystems.WktStreamTokenizer tokenizer) (at <45d4f050441b42ac9796be647bb4c7ef>:0)
ProjNet.IO.CoordinateSystems.CoordinateSystemWktReader.ReadCoordinateSystem (System.String coordinateSystem, ProjNet.IO.CoordinateSystems.WktStreamTokenizer tokenizer) (at <45d4f050441b42ac9796be647bb4c7ef>:0)
ProjNet.IO.CoordinateSystems.CoordinateSystemWktReader.Parse (System.String wkt) (at <45d4f050441b42ac9796be647bb4c7ef>:0)
ProjNet.CoordinateSystems.CoordinateSystemFactory.CreateFromWkt (System.String WKT) (at <45d4f050441b42ac9796be647bb4c7ef>:0)

The HorizontalCoordinateSystem only supports two axes so the error mostly makes sense.

Here is where things get a bit tricky...

I found this discussion about EPSG 4979 as well: https://github.com/OSGeo/PROJ/issues/394. There is mention of WKT v1 not supporting a third axis on a GEOGCS definition.

Both major versions of the WKT spec have a comment stating that 3D GEOGCS appears in OGC 06-103r4 as an extension where the third axis is optional. However, skimming through that document, I can't seem to find it (granted I am not very familiar with these specifications...). The WKT specifications themselves though seem to call out the vertical axis as seen below.

In WKT v1, under section 7.5, there is this statement:

iii. For geodetic CRSs having a three-dimensional ellipsoidal coordinate system, the name and abbreviation of the horizontal axes in a WKT string shall follow the requirements in (ii). The vertical axis name shall be ‘ellipsoidal height’; the vertical axis abbreviation shall be ‘h’ and should be included when abbreviations for the horizontal axes are included.

In WKT v2, under section 7.5.3, there is this statement:

d) For geographic CRSs having a three-dimensional ellipsoidal coordinate system, the name and abbreviation of the horizontal axes in a WKT string shall follow the requirements in c). The vertical axis name shall be 'ellipsoidal height'; the vertical axis abbreviation shall be 'h' and shall be included when abbreviations for the horizontal axes are included.

Also, in the same (or similar) sections for each version (v1 - 7.5, v2 - 7.5.2), there is a table that specifies the dimensions for the ellipsoidal type as 2 or 3.

It seems to me like it would be a good idea to support the optional third axis on GeographicCoordinateSystem though I am not sure what ramifications that would have on other parts of the library.

On an unrelated note, that WKT representation for EPSG 4979 confuses me a bit since I don't understand how it differentiates the units between the axes. The lat/lon (north/east) axes are in degrees which I see in the WKT but the height (up) should be in metres according to the EPSG entry but I don't see that defined anymore in the WKT representation. Maybe I am missing a default or assumption that gets made.

pchilds commented 3 years ago

Well Proj4 was updated to support this https://github.com/OSGeo/PROJ/issues/394 Would have to merge in these underlying changes, but if so may as well go the whole hog and target Proj 6+ (https://github.com/NetTopologySuite/ProjNet4GeoAPI/issues/89) rather than fighting all these little fires.

4979 seems to be a bit of an exception in how it handles 3d. For other cases, e.g. AHD_1971 geoid one normally sets up a COMPD_CS with PROJCS and VERT_CS members (COMPD_CS and VERT_CS are not supported by ProjNET though).

Incidentally there is ANGLEUNIT and LENGTHUNIT when one cannot implicitly derive which form UNIT takes. Incidoublydentally ProjNET does not support these.