[solved] Georeferencing LAS file for ION

Cesium ion
1

Hi!

I’m trying to get custom point data using LAS into Cesium, however, they’re ending up being slightly off…

I have a bunch (hundreds) of GPS tracks in kml that I process, project to WGS84, rotate and move to a local coordinate system. I write my point data into an LAS file using laspy. However, when I load the data into ION and set the origin to what I used for my local coordinate system, it does not respect the origin of the dataset.

Any suggestions how to get custom point data into Cesium, or how tiles are generated from the LAS files?

Cheers

2

Hi @juraj, thanks for your question and welcome to the forum! Would you be able to share some of your data with us and more details about your process (e.g., steps for projecting, rotating, and moving the data), so we can investigate internally? We can try to reproduce the problem on our end.

Also, how far off are the point clouds from their correct location? If you have a screenshot, that would be helpful.

1 Like
3

Hi @Matt_Boyd-Surka !

Thanks a lot for a quick reply! I’ve set up a small example here: GitHub - jurajpalenik/las-to-cesium: Getting custom point data into Cesium

For the wgs projection I re-implemented the Cesium Cartesian3.fromRadians function into python. [1]

For the rotation I use the spherical normal, tangent and bi-tangent. This I’m not 100% sure about, however, the fromRadians function is using the spherical normal to offset points on the ellipsoid, so it should correspond.

For the centering I simply project a reference location and subtract this from the data points.

In the end I shuffle the axis such that Z is altitude, X is longitude, Y is latitude.

Here is a screenshot of the misplaced data.

result.PNG
result.PNG707×1464 377 KB

I did a bunch of experiments with re-centering the data and I suspect that ION might be computing a bounding box of the data. Centre of the bounding box would be then used as the new origin. A clue towards this is the position of the 3D cursor when placing the dataset. I’ve centred the data towards the centre of mass (using mean function), however, the cursor appears in the middle of the BB. Centering the data around the middle of the BB got me somewhat closer, however, still not the same position.

Let me know if you have more questions. :slight_smile:

Cheers
Juraj

[1] cesium/Cartesian3.js at 1.83 · CesiumGS/cesium · GitHub

4

@juraj Thank you for putting together the example and providing the screenshot and other info—that’s all very helpful. Just wanted to let you know that we’re investigating this, and I’ll update you when we have a result or if we have any more questions.

5

Hi @juraj, one related question I should’ve asked earlier: What is your goal with converting the KML files to LAS? CesiumJS can bring in KML files as data sources without needing to convert, if that is helpful.

6

Hi @Matt_Boyd-Surka !

Thank you very much for looking into the issue! I’m sorry I didn’t see your message earlier. Yeah, I know that Cesium can show KML data directly. That’s actually how I check whether the points are displayed correctly.

My workflow is that I process multiple KML flights and generate point cloud data that I would like to geo-locate and display using different visual encodings. At the moment I am using the Cesium3DTilesetStyle property to set size and color of the points. I believe I should also be able to animate the dataset by providing data for multiple timesteps.

Using the LAS file was a convenience choice as I was working with the laspy library before and noticed that Cesium ION supported loading of point cloud datasets in the LAS format. In theory, I should be able to display the point data using some other format like .pnts, .b3dm, or .gltf. However, I am not sure how I would tell Cesium to billboard the point data loaded from .gltf, as there are no faces to be rendered in the point cloud.

Cheers
Juraj

7

Hi @juraj, thanks for clarifying, that’s helpful. The workflow I would recommend is that you embed the CRS of your data in the LAS file before tiling it in ion. When I ran the Python script you provided, I saw that it output 6.40429397482618 lon, 60.65092204339487 lat, 0 alt as the geodetic origin of the data. When I passed this origin into the tiler manually via command line, the KML and LAS data were aligned after tiling. (Almost. When I zoom way in, I still see a slight difference, and I’m wondering if that is the result of error being introduced during the initial processing.)

Some information about the tiler and LAS files:

The variable length records may contain information about the Coordinate Reference System (CRS); in LAS 1.3 and earlier the CRS is stored as GeoTIFF keys and in LAS 1.4 the CRS may be stored as Well Known Text. The tiler supports all versions of the spec, 1.0 through 1.4, though most LAS files are either LAS 1.2 or LAS 1.4.

I’m not familiar with laspy, but hopefully it allows you to put CRS information in the LAS header. An alternative solution would be to store the points in ECEF.

The tiler may recompute the bounding box, so the original bounding box can’t necessarily be relied upon to georeference the data after the fact. Positions are stored relative-to-center.

Hopefully this information is helpful, but please feel free to ask any follow-up questions. I’m happy to clarify or investigate further.

For your reference, below is the tileset.json file outputted by the point-cloud tiler when no CRS or origin is supplied:

{
   "asset":{
      "extras":{
         "ion":{
            "georeferenced":false,
            "movable":true
         }
      },
      "version":"1.0"
   },
   "geometricError":4958.606108514771,
   "root":{
      "boundingVolume":{
         "box":[
            0.0,
            0.0,
            613.4449999999999,
            1223.7150000000001,
            0.0,
            0.0,
            0.0,
            2067.16,
            0.0,
            0.0,
            0.0,
            613.4449999999999
         ]
      },
      "children":[
         {
            "boundingVolume":{
               "box":[
                  -263.78499999999985,
                  755.9499999999998,
                  976.365,
                  1223.7150000000001,
                  0.0,
                  0.0,
                  0.0,
                  2067.16,
                  0.0,
                  0.0,
                  0.0,
                  613.4449999999999
               ]
            },
            "content":{
               "uri":"0/0.pnts"
            },
            "geometricError":0.0,
            "transform":[
               1.0,
               0.0,
               0.0,
               0.0,
               0.0,
               1.0,
               0.0,
               0.0,
               0.0,
               0.0,
               1.0,
               0.0,
               263.78499999999997,
               -755.95,
               -362.92,
               1.0
            ]
         }
      ],
      "geometricError":4958.606108514771,
      "refine":"ADD"
   }
}

And here it is after I specified the geodetic origin (note the “transform” block toward the end):

{
   "asset":{
      "extras":{
         "ion":{
            "georeferenced":true,
            "movable":true
         }
      },
      "version":"1.0"
   },
   "geometricError":4958.606108514771,
   "root":{
      "boundingVolume":{
         "box":[
            0.0,
            0.0,
            613.4449999999999,
            1223.7150000000001,
            0.0,
            0.0,
            0.0,
            2067.16,
            0.0,
            0.0,
            0.0,
            613.4449999999999
         ]
      },
      "children":[
         {
            "boundingVolume":{
               "box":[
                  -263.78499999999985,
                  755.9499999999998,
                  976.365,
                  1223.7150000000001,
                  0.0,
                  0.0,
                  0.0,
                  2067.16,
                  0.0,
                  0.0,
                  0.0,
                  613.4449999999999
               ]
            },
            "content":{
               "uri":"0/0.pnts"
            },
            "geometricError":0.0,
            "transform":[
               1.0,
               0.0,
               0.0,
               0.0,
               0.0,
               1.0,
               0.0,
               0.0,
               0.0,
               0.0,
               1.0,
               0.0,
               263.78499999999997,
               -755.95,
               -362.92,
               1.0
            ]
         }
      ],
      "geometricError":4958.606108514771,
      "refine":"ADD",
      "transform":[
         -0.11154340883487353,
         0.9937595624422921,
         0.0,
         0.0,
         -0.866210285377921,
         -0.09722678568387881,
         0.4901292621860905,
         0.0,
         0.4870706411302129,
         0.05467068867395802,
         0.8716497612854136,
         0.0,
         3114085.343231913,
         349271.51667731546,
         5537068.5816861205,
         1.0
      ]
   }
}
8

Hi @Matt_Boyd-Surka!

Sorry for a late reply.

Thank you very much, this was very helpful!

I ended up using ECEF in the WGS-84 projection and store those in the LAS, specifying the coordinate reference system in the variable length record using the WKT and the LAS format version 1.4.

If anyone finds this thread, the code in the github project is updated and correctly georeferences the data.

Cheers
Juraj

1 Like
9

@juraj Thank you very much for reporting back with your solution and for sharing your project! The community really appreciates it!