GDAL LiDAR – Reprojecting LiDAR Data with libLAS and Troubleshooting Latitude or Longitude Limit Errors

coordinate systemerrorgdalliblaslidar

I have a lidar file (file.las) with the following spatial reference:

Spatial Reference:           
PROJCS["NAD83 / UTM zone 11N",
    GEOGCS["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.2572221010002,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6269"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4269"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-117],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","26911"]]

Then I reproject it to WGS84 with las2las:

las2las --a_srs EPSG:26911 --t_srs EPSG:4326 -i file1.las -o output.las

And I get: 

Spatial Reference:           
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]]

But when I try to revert the transformation: from WGS84 to NAD83/UTM 11N, I get an error:

las2las --a_srs EPSG:4326 --t_srs EPSG:26911 -i output.las -o wgs2utm.las
ERROR 1: latitude or longitude exceeded limits
error: Could not project point for ReprojectionTransform::latitude or longitude exceeded limits0.

Why I got this error, if I am just reversing the transformation (originally from NAD83 to WGS84, and then from WGS84 to NAD83)?.

Best Answer

Set the scaling and offset when reprojecting to WGS84, e.g.:

las2las --a_srs EPSG:26911 --t_srs EPSG:4326 -i file1.las -o output.las --scaling 1e-7 1e-7 0.01 --offset <something close to your data's longitudes>,<something close to your data's latitudes>,0

Explanation

You've been caught by a limitation/feature of the las file format. Internally, a las file stores each point's spatial dimensions (x, y, z) as a combination of per-point long values (4 byte signed integer value) and file-global scale and offset values. When reading a las file, each point's x, y, and z values are calculated by multiplying that point's signed integer value by the file-global scaling factor, then adding the file-global offset:

x = x_in_file * x_scaling + x_offset

When you reproject the data from EPSG:26911 to EPSG:4326, as you did, you are not instructing las2las to use any new scaling or offset values; therefore, las2las keeps the same scaling and offset values that it read from the source file. When going between two projected coordinate systems (e.g. UTM to state plane) this might be ok, but when going between projected and geodetic coordinate systems (e.g. UTM to WGS84) this is almost assuredly going to give you wonky results.

Examples

To illustrate this, I ran a couple of simple tests on a EPSG:26915 test file contained in the PDAL source repository (https://github.com/PDAL/PDAL/blob/master/test/data/las/utm15.las). First, a simple reprojection with no scaling and offset:

$ las2las -i utm15.las -o wgs84.las --t_srs EPSG:4326 && lasinfo wgs84.las | grep "Bounding Box"
  Bounding Box:     -93.35, 41.58, -93.35, 41.58

You can see that the data contained in my wgs84.las file has been clipped to two decimals of precision. That means my lasfile has a data resolution of about 0.8km — not good!

Scaling

We can dial up the precision with the --scale option to las2las:

$ las2las -i utm15.las -o wgs84.las --t_srs EPSG:4326 --scale 0.0000001 0.0000001 0.01 && lasinfo wgs84.las | grep "Bounding Box"
  Bounding Box:     -93.3515626, 41.5771484, -93.3515626, 41.5771484

Now that's better. I chose 0.0000001 (1e-7) for my x and y scaling factors because this gives me a bit less than 1cm precision at this latitude and longitude. You maybe have to adjust your scaling factors if you're in a place where the longitude lines get pretty close together (e.g. northern Alaska) but in general you can use 1e-7 or 1e-8 and things should work. I use the haversine equation to calculate the distance precision — specifically https://github.com/mapado/haversine.

Offset

This will work fine and dandy if your data are constrained to a small spatial area. However, at a certain point (scaling-dependent) a 4-byte signed integer will not be able to represent all of the values in your data. At this point you will need to offset your lasfile — think of it like "rezeroing" the origin of the lasfile closer to the center of your data.

Calculating a good offset usually requires some a-priori knowledge — I usually take the minimum latitude and longitude values for my data, truncate those values to the ones place, and use the resulting integers as my offset. For my test data file, this would look like:

$ las2las -i utm15.las -o wgs84.las --t_srs EPSG:4326 --scale 1e-7 1e-7 0.01 --offset -93,41,0

Note the different syntax for offset and for scaling — this is a strange quirk of libLAS.

Round trip

At this point, you should be able to create a properly scaled and offset WGS84 lasfile from your UTM data. Let's make sure by roundtripping our wgs84 file back to its original projected coordinate system (not including any offset information this time):

$ las2las -i wgs84.las -o utm15-roundtrip.las --t_srs EPSG:26915 --scale 0.01 0.01 0.01

Throwing out the offset is another weird quirk of libLAS that might actually be a bug (though I need to investigate further). Anyways, it shouldn't hurt your data to have that little offset in the roundtrip UTM file.

Related Question