[Liblas-devel] Precision of reprojection transforms

[Howard Butler]

On Apr 29, 2011, at 9:03 AM, Frank Warmerdam wrote:

> On 11-04-29 08:35 AM, Edwards, David wrote:
>> I am trying to use libLAS to read an LAS file and perform a coordinate
>> transformation.  Points in the file are stored in UTM coordinates.  I would
>> like to transform them into geocentric coordinates.  As best I can tell, the
>> version of libLAS I am using (1.6.1) does not support conversion to
>> geocentric coordinates, so my backup plan is to convert to WGS84 geodetic
>> coordinates (lat., long. , height above ellipsoid), and then perform the
>> geodetic-to-geocentric conversion myself.
>> 
>> The problem is that when I set up a ReprojectionTransform to convert the
>> points from UTM to geodetic coordinates, the resulting latitude and
>> longitude of each point has nonzero digits only to the hundredth of a
>> degree, even though the UTM coordinates of my input points have data to the
>> millimeter.  For example, here are the results of two such conversions:
>> 
>> UTM (Easting, Northing) --->  Geodetic (Longitude, Latitude) 831805.86
>> 4296378.26 --->  -77.18 38.75 831806.31 4296377.64 --->  -77.18 38.75
>> 
>> Even though the points have different UTM easting and northing coordinates,
>> they are transformed to the same latitude and longitude.  Is this limited
>> precision inherent in the way libLAS does its coordinate conversions, or am
>> I setting up the transform object(s) incorrectly?
>> 
>> I apologize if this is a duplicate post.  I tried to send this yesterday
>> before joining the mailing list.
> 
> Dave,
> 
> Reprojection is generally accurate to a small fraction of a millimeter.
> I believe the problem you are encountering is how you are formatting the
> coordinates for reporting - not the nature of the reprojection itself.

The problem is related to the formatting, but it is actually the data scaling, not the formatting.  X, Y, and Z coordinates in LAS are stored as scaled integers.  Your X, Y, Z values of your UTM data have two decimals of precision as defined by the scale and offset in the LAS file's header.  Reprojecting the data into DD causes the data to be descaled into doubles, reprojected, and then scaled back by the scale/offset values of the point's header (in this case 0.01 for X/Y).  This results in not enough precision, obviously.

You need to tell the ReprojectionTransform about what it should use for the scale/offset of the output data.  You can do this by creating a copy of your header, setting the scale/offset to what you need, and then giving a shared_ptr of that to the ReprojectionTransform constructor.  See line 1332 of apps/laskernel.cpp for an example.

Hope this helps,

Howard