Vertical Coordinate Systems in .las, GeoTIFF and WKT

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

Vertical Coordinate Systems in .las, GeoTIFF and WKT

Frank Warmerdam
Folks,

I have recently done some work on vertical cs support in the .las (lidar)
format.  LAS uses geotiff tags for georeferencing, and liblas also returns
OGC WKT representations of coordinate system.

I have prepared a web page on my efforts, the first portion of which (down
about as far as "Sample Data") may be of broader interest.

The support for vertical coordinate systems is already covered in the
OGC CT specification (01-009) and the GeoTIFF specification, but in my
opinion little use has been made of this in the past.  So I'm interested
in feedback on my approach, and I'm hoping it can help as we move to
support vertical coordinate systems more widely in GeoTIFF, and other
coordinate system projects.

   http://liblas.org/wiki/VerticalCS

Best regards,
--
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, [hidden email]
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | Geospatial Programmer for Rent

_______________________________________________
Geotiff mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/geotiff
Reply | Threaded
Open this post in threaded view
|

Re: Vertical Coordinate Systems in .las, GeoTIFF and WKT

Max Martinez
Frank,

The description of your approach doesn't surprise me (sounds very
reasonable) but I wondered how you envisioned dealing with the
"ellipsoidal height" issue which is compounded (no pun intended) by the
use of GeoTIFF.

The issue is that I didn't think EPSG/ISO 19111/OGC support vertical
coordinate systems based on ellipsoidal height. Although the original
OGC specification may have allowed these vertical datum types, ISO
19111/OGC AT 2 does not and I don't believe it is possible to have a
standard vertical coordinate system of this type in the EPSG database.
This not only makes it impossible to describe a projected coordinate
system with ellipsoidal height, but also causes problems in the purely
geographic case (because you still can't do 2d geographic + vertical
based on ellipsoid). Are you envisioning people using the 3d geographic
crs codes and omitting the vertical information?

This has been a topic before on the GeoAPI mailing lists (which
continues to support a geoidal vertical datum type even in their
proposed 3.0 version which is intended for OGC standardization) but
remains somewhat unresolved. Roel Nicolai was involved in that
discussion but I didn't quite get a satisfactory answer for why it has
been eliminated from the model. I've attached that conversation (in case
I'm being dense).

I don't know if it is an issue for las data, but in the general
promotion of the vcs information in GeoTIFF, it is likely to come up and
cause interoperability issues unless an explicit statement is made about
how to handle lat, lon, ellipsoidal height and whether to and how to
handle projected + ellipsoidal height in geotiff. (And, what happens
when translating from WKT [where it is expressible] to GeoTIFF in these
cases)?

Max

-----Original Message-----
From: [hidden email]
[mailto:[hidden email]] On Behalf Of Frank Warmerdam
Sent: Monday, January 11, 2010 1:51 PM
To: [hidden email]; GeoTIFF
Subject: [Geotiff] Vertical Coordinate Systems in .las, GeoTIFF and WKT

Folks,

I have recently done some work on vertical cs support in the .las
(lidar)
format.  LAS uses geotiff tags for georeferencing, and liblas also
returns
OGC WKT representations of coordinate system.

I have prepared a web page on my efforts, the first portion of which
(down
about as far as "Sample Data") may be of broader interest.

The support for vertical coordinate systems is already covered in the
OGC CT specification (01-009) and the GeoTIFF specification, but in my
opinion little use has been made of this in the past.  So I'm interested
in feedback on my approach, and I'm hoping it can help as we move to
support vertical coordinate systems more widely in GeoTIFF, and other
coordinate system projects.

   http://liblas.org/wiki/VerticalCS

Best regards,
--
---------------------------------------+--------------------------------
------
I set the clouds in motion - turn up   | Frank Warmerdam,
[hidden email]
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | Geospatial Programmer for Rent

_______________________________________________
Geotiff mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/geotiff

Roel,

> As soon as you want to use elliposidal heights in a different CRS you need to transform them using the orginal latitude and longitude belonging to each height, which is why ISO 19111 won't let you separate them from the heights.

My basis for calling the argument weak is my (perhaps mis-) impression that a similar (though not identical) potentially error-prone procedure is required for transforming geoidal heights that are part of a compound (projected + vertical) CRS...which is, in fact, allowed by the specification. What is it about the geoidal case that allows it to be viewed as less problematic in this regard?

Max

-----Original Message-----
From: [hidden email] [mailto:[hidden email]]
Sent: Wednesday, June 03, 2009 4:41 AM
To: Max Martinez; [hidden email]; [hidden email]
Cc: [hidden email]; [hidden email]
Subject: RE: [GeoAPI.wg] Compound CRS Question: Does 2D Projected + Ellipsoidal Height not make sense?

Max,

Geodetically the argument is not weak at all.
Software developers have a different approach than geodesists. Technical integrity of software is a different thing than geodetic technical integrity.
If your software doesn't fall over, then you have a software-technically good product. However, it may have given you the wrong answer geodetically and you may not even notice.

The revision of ISO 19111 that was initiated from OGC, and resulted in the current version, was primarily aimed at enabling interoperabilty of geospatial software, but ISO 19111 also contains geodetic constraints, some of which boil down to strongly recommended best practices, others prevent downright geodetic errors. One example of the latter is that you are not free to choose the ellipsoid with the Geodetic Datum, a fairly straightforward case.
Much debate has arisen over the axis order and although I do not want stir up the debate again, I will give an example: we had a contractor, working in the Arabian Gulf, who fed geographic coordinates into a coordinate transformation algorithm in the order (longitude, latitude), whilst the software expected (latitude, longitude). The software didn't keel over, but its results were in error by a few hundred metres.

Using ellipsoidal heights on their own is not a geodetic error to begin with, but it very quickly leads to errors. Geodesists are aware of this but non-geodesists don't see the danger and cannot be blamed for that. This discussion is a case in point.
As soon as you want to use elliposidal heights in a different CRS you need to transform them using the orginal latitude and longitude belonging to each height, which is why ISO 19111 won't let you separate them from the heights. To just warn a user who has little or no geodetic knowledge that stripping off latitude and longitude and only retain ellipsoidal height is 'dangerous' is not going to help:  the user will believe that if the software hasn't keeled over after he did that, the 'danger' didn't materialise.
I have plenty of examples in Shell where that is the case: we use ArcGIS and that used to warn the user upon import of a dataset with a different CRS that he/she will need to transform the dataset. Most of our users would click "OK' after which the software did not transform, but happily merged the dataset with the others in the project, expecting that the user had been warned enough about it. Some users would even tick the box "Don't warn me about this again", as they found it a nuisance warning message ("This happens every time I import a dataset, but nothing ever seems to actually go wrong!"). This is not a opportunity to gratuitously bitch at one piece of software, because there is hardly any geospatial software around that maintains geodetic integrity for 100%. Apparently geodesy is a very hard-to-understand subject.
For that reason ISO 19111 precludes some options that software developers would happily allow and non-geodetic users would happily use.


Best regards,
Roel

Roel Nicolai
Principal Geodesist

Shell International Exploration and Production B.V.
The Hague, The Netherlands - Trade Register no. 27002688
Data Management and Geomatics
Kessler Park 1, 2288 GS Rijswijk
(PO Box 60, 2280 AB Rijswijk)
The Netherlands
Tel: +31 70 447 3652
Tel (mob): +31 6 1097 3652
Fax: +31 70 447 3695
Email: [hidden email]
Internet: http://www.shell.com

This message, any attachment and response string are confidential and may be legally privileged. It is intended only for the use of the parties to whom it is addressed. If you are not the addressee indicated in this message please notify the sender immediately by reply email and destroy this message. All information and attachments remain the property of Shell.


-----Original Message-----
From: Max Martinez [mailto:[hidden email]]
Sent: 02 June 2009 16:52
To: Nicolai, Roel R SIEP-EPT-SCD; [hidden email];
[hidden email]
Cc: [hidden email];
[hidden email]
Subject: RE: [GeoAPI.wg] Compound CRS Question: Does 2D Projected +
Ellipsoidal Height not make sense?


Thanks, Roel. That has definitely clarified the intention of ISO 19111.

I would comment, however, that the argument for not allowing this appears to me (I'll emphasize that I am a non-geodesist) to be weak. You don't seem to be saying "you shouldn't do it because you can't make it work" but rather "it's disallowed because you have to be careful to make it work". And "make it work" means a "correct" CRS transformation which will always introduce error anyway even under ideal circumstances. Am I misinterpreting the answer below?

I think software developers typically have a different prioritization of "practical" matters (as pointed out by Martin's messages on this subject) which puts the consideration "is there already data like this out there" at the top of the list. Secondarily, you warn the user when they might be doing something non-kosher (something that "may be incorrect"), but you let them do it anyway unless you KNOW it is incorrect.

Max

-----Original Message-----
From: [hidden email] [mailto:[hidden email]]
Sent: Tuesday, June 02, 2009 8:11 AM
To: Max Martinez; [hidden email]; [hidden email]
Cc: [hidden email]; [hidden email]
Subject: RE: [GeoAPI.wg] Compound CRS Question: Does 2D Projected + Ellipsoidal Height not make sense?

Hello Martin, Max,

I have been away on leave for a week and a bit, so am picking up the conversation only now.

ISO 19111 explicitly 'forbids' the use of ellipsoidal heights on their own, i.e. as a separate Vertical CRS.  The geodetic reason for that is, that ellipsoidal height is an inseparable part of a 3D coordinate tuple. Positions determined by GPS give you Geocentric (X, Y, Z) to start with. This coordinate tuple can be converted to (lat, lon, ell height).
By inseparable I mean geodetically inseparable, because of course anyone manipulating a data file can split off the ellipsoidal heights. Similarly you cannot split off the latitude; latitudes do not mean much on their own. A analagous but less strong and less obvious case exists for ellipsoidal heights.

The practical reason for ISO 19111 not permitting ellipsoidal heights on their own is that it becomes too easy to make errors, not just as a secondary effect in the horizontal position as described by Martin, but primarily in the interpretation of the height numbers. If the Projected CRS you wish to use is associated with a different Geodetic Datum that the ellipoidal heights are and you do not realise that, you make interpretation errors. The ellipsoidal heights would have to be transformed to the new Geodetic Datum first and you can only do that by rejoining them first with their original latitudes and longitudes. So why split them off to begin with?
Another potential problem is that your new Geographic CRS is by definition a 2D system, such as e.g. European Datum 1950 (ED50). By definition ellipsoidal heights on ED50 do not exist.
However, the biggest risk is that people will NOT transform the elliposidal heights, because they believe they have an intrinsic spatial meaning on their own and then you have a dataset with errors.

A Projected CRS + ellipsoidal height cannot be modelled in an ISO 19111-compliant way, so I can hardly advise you how to do it. In that context Martin's statement that:
"...ISO 19111 wants to force us to express "2D Projected CRS + ellipsoidal height" as a ProjectedCRS backed by a 3D cartesian CS, in the same way that a "2D Geographic CRS + ellipsoidal height" is backed by a GeographicCRS backed by a 3D ellipsoidal CS."
... is incorrect. ISO 19111 tells you not to do that at all.

Martin mentions that this practice is happening anyway. I would be surprised if WKT allows the expression of a Compound CRS consisting of a Projected CRS and an elliposidal height (the term Projected 3D CRS does not exist in ISO 19111 !!). If this appears in some OGC document, please let me know, because it will have to be changed.

The trouble with all these things is that you may be able to do something mathematically, or as part of data manipulation, i.e. you will not an error message like "division by zero". However, the operation may be incorrect in the geodetic sense, which means that after your calculations, when you make the association back to the real world, you may end up with position errors. The idea of the modeling constraints that exist in ISO 19111 is to prevent the conditions that may lead to such errors.  

Perhaps superfluously I add a drawing from one of my presentations on the consequences of not transforming ellipsoidal heights.


Best regards,
Roel

Roel Nicolai
Principal Geodesist

Shell International Exploration and Production B.V.
The Hague, The Netherlands - Trade Register no. 27002688
Data Management and Geomatics
Kessler Park 1, 2288 GS Rijswijk
(PO Box 60, 2280 AB Rijswijk)
The Netherlands
Tel: +31 70 447 3652
Tel (mob): +31 6 1097 3652
Fax: +31 70 447 3695
Email: [hidden email]
Internet: http://www.shell.com

This message, any attachment and response string are confidential and may be legally privileged. It is intended only for the use of the parties to whom it is addressed. If you are not the addressee indicated in this message please notify the sender immediately by reply email and destroy this message. All information and attachments remain the property of Shell.




-----Original Message-----
From: Max Martinez [mailto:[hidden email]]
Sent: 02 June 2009 03:15
To: Martin Desruisseaux; GeoAPI SWG
Cc: Nicolai, Roel R SIEP-EPT-SCD; [hidden email];
[hidden email]
Subject: RE: [GeoAPI.wg] Compound CRS Question: Does 2D Projected +
Ellipsoidal Height not make sense?


Thanks, Martin,

I'm not 100% sure I know exactly what you mean by "backed by a 3D Cartesian CS". I guess I could envision, instead of the usual 2D conversion operation from the base 2D geographic there being some sort of 3D passthrough conversion operation from a base 3D Geographic where the ellipsoidal height passes through unchanged. Is that what you mean? (thus the CS of the target of the passthrough operation would be a 3D Cartesian). There doesn't seem to be a lot of direct support in ISO 19111 for that. Table 16 on page 25 of ISO 19111 2007 doesn't even allow for a third axis when the CRS type is projected. Perhaps a simple omission but in any case its confusing.

Max

-----Original Message-----
From: Martin Desruisseaux [mailto:[hidden email]]
Sent: Wednesday, May 27, 2009 4:53 PM
To: GeoAPI SWG
Cc: Max Martinez; [hidden email]; [hidden email]; [hidden email]
Subject: Re: [GeoAPI.wg] Compound CRS Question: Does 2D Projected + Ellipsoidal Height not make sense?

Hello Max

Max Martinez a écrit :
> On a separate email thread, the issue of the inability to easily
> formulate a compound CRS using a 2D Projected CRS for the horizontal and
> a 1D ellipsoidal height system for the vertical using the ISO 19111
> model caused me to wonder "does this not make sense"? If not, is there a
> simple explanation that the geodesist can give to the layman? If, on the
> other hand, it does make sense, what would be the recommended way to
> construct such a beast in the ISO 19111 model?

If I were to guess (please Roel correct me if I'm wrong) I would said that ISO
19111 wants to force us to express "2D Projected CRS + ellipsoidal height" as a
ProjectedCRS backed by a 3D cartesian CS, in the same way that a "2D Geographic
CRS + ellipsoidal height" is backed by a GeographicCRS backed by a 3D
ellipsoidal CS.

I would guess that they want to "force" as a matter of principle, since most
mathematical operation like conversions from Geographic to Geocentric CRS
require a 3D CRS. For example when doing a transformation from a 3D Geographic
CRS to an other 3D Geographic CRS using the Molodenski transformation, a change
in the ellipsoidal height produces slight changes in the resulting latitudes and
longitudes. If the 3D Geographic CRS are actually implemented as a compound of
2D Geographic CRS + 1D ellipsoidal height, it would be very easy for an
implementor to just forget the height, thus introducing a slight error in the
result. With ISO 19111 forcing us to a "true" 3D CRS, this risk or error may be
slightly reduced.

The above was for Geographic CRS, but maybe by extension a similar rational can
be imagined for Projected CRS.

However in practice a few difficulties arise. There is this "Well Known Text"
(WKT) format in the wild, which is widely used. That format express 3D projected
CRS as a compound of 2D CRS + 1D height (ellipsoidal or not). When parsing such
WKT, an implementation may need to *temporarily* be able to express that
construct, even if it is going to be merged in a "real" 3D CRS later. A possible
compromise could be that an implementation may very well create such constructs
internally, but they should be as short lived as possible and replaced by 3D
single CRS before being visible to the user.

I'm not sure if I'm answering the question, it is just some thoughs that could
be developped more if I'm not going in the wrong direction...

        Regards,

                Martin

P.S.: We have been silent on the OGC mailing list (there is a little bit more
trafic on the SourceForge list), but we are in the process of creating a
specification draft to be submitted to the next OGC meeting at Boston. A time
slot for GeoAPI has been reserved at Wednesday 1500-1600 hrs.




_______________________________________________
GeoAPI.wg mailing list
[hidden email]
https://lists.opengeospatial.org/mailman/listinfo/geoapi.wg

_______________________________________________
Geotiff mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/geotiff