Reproyección de datos LiDAR con libLAS, error: ¿"latitud o longitud exceden los límites"?

Question

Tengo un archivo lidar (archivo.las) con la siguiente referencia espacial:

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"]]

Luego lo reproyecto a WGS84 con las2las:

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

Y lo consigo:

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"]]

Pero cuando intento revertir la transformación: de WGS84 a NAD83/UTM 11N, obtengo un 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.

¿Por qué tengo este error, si sólo estoy invirtiendo la transformación (originalmente de NAD83 a WGS84, y luego de WGS84 a NAD83)?

Answer 1

Establezca la escala y el desplazamiento al reproyectar a WGS84, por ejemplo:

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

Explicación

Te ha pillado una limitación/característica del formato de archivo las. Internamente, un archivo las almacena las dimensiones espaciales de cada punto (x, y, z) como una combinación de valores largos por punto (valor entero con signo de 4 bytes) y valores de escala y desplazamiento globales del archivo. Al leer un archivo las, los valores x, y, z de cada punto se calculan multiplicando el valor entero con signo de ese punto por el factor de escala global del archivo, y añadiendo después el desplazamiento global del archivo:

x = x_in_file * x_scaling + x_offset

Al reproyectar los datos de EPSG:26911 a EPSG:4326, como ha hecho, no está indicando las2las para utilizar cualquier nuevo valor de escala o desplazamiento; por lo tanto, las2las mantiene los mismos valores de escala y desplazamiento que leyó del archivo fuente. Cuando se pasa de dos sistemas de coordenadas proyectados (por ejemplo, de UTM a plano de estado) este puede pero al pasar de un sistema de coordenadas proyectado a otro geodésico (por ejemplo, de UTM a WGS84) es casi seguro que se obtendrán resultados erróneos.

Ejemplos

Para ilustrar esto, ejecuté un par de pruebas sencillas en un archivo de prueba EPSG:26915 contenido en el repositorio de fuentes de PDAL ( https://github.com/PDAL/PDAL/blob/master/test/data/las/utm15.las ). En primer lugar, una simple reproyección sin escala ni desplazamiento:

$ 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

Puede ver que los datos contenidos en mi archivo wgs84.las han sido recortados a dos decimales de precisión. Eso significa que mi archivo las tiene una resolución de datos de unos 0,8 km - ¡no es bueno!

Escala

Podemos aumentar la precisión con el botón --scale a 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

Eso está mejor. Elegí 0,0000001 (1e-7) para mis factores de escala x e y porque esto me da un poco menos de 1cm de precisión en esta latitud y longitud. Tal vez tengas que ajustar tus factores de escala si estás en un lugar donde las líneas de longitud están muy juntas (por ejemplo, el norte de Alaska), pero en general puedes usar 1e-7 o 1e-8 y las cosas deberían funcionar. Yo utilizo la ecuación del haversino para calcular la precisión de la distancia, concretamente https://github.com/mapado/haversine .

Offset

Esto funcionará perfectamente si sus datos están limitados a un área espacial pequeña. Sin embargo, a partir de cierto punto (que depende de la escala), un entero con signo de 4 bytes no podrá representar todos los valores de los datos. En este punto, tendrá que desplazar su archivo de datos láser - piense que es como "volver a poner a cero" el origen del archivo de datos láser más cerca del centro de los datos.

El cálculo de un buen desplazamiento suele requerir cierto conocimiento a priori: yo suelo tomar los valores mínimos de latitud y longitud de mis datos, truncar esos valores a la unidad y utilizar los enteros resultantes como desplazamiento. Para mi archivo de datos de prueba, esto sería así:

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

Nótese la diferente sintaxis para el desplazamiento y para el escalado - esto es una extraña peculiaridad de libLAS.

Viaje de ida y vuelta

Llegados a este punto, debería ser capaz de crear un lasfile WGS84 correctamente escalado y desplazado a partir de sus datos UTM. Asegurémonos de que nuestro archivo WGS84 vuelve a su sistema de coordenadas proyectado originalmente (sin incluir ninguna información de desplazamiento esta vez):

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

Desechar el desplazamiento es otra extraña peculiaridad de libLAS que podría ser un error (aunque tengo que investigar más). De todas formas, no debería perjudicar a tus datos tener ese pequeño desplazamiento en el archivo UTM de ida y vuelta.