Subscribe to this thread
Home - General / All posts - Linked LiDAR point clouds - accessing true X, Y, Z values
LandSystems73 post(s)
#31-Mar-20 05:41

With the current situation, I am finally finding time to start looking a bit further into some of the LAS tools that have been appearing over recent builds. My aim is to be able to use M9 for performing a series of QA checks on a whole of region LiDAR supply forming part of New Zealand’s National Elevation Model (https://www.linz.govt.nz/data/linz-data/elevation-data).

Thus far, I have found the librarylas dataport used in-conjunction with the ‘Index LiDAR files’ tool to be great additions that work really well. I have found however, that some of my LAS files do not contain projection information so by way of an improvement, it would be nice to be able to set a projection for all files using the index tool as in my case a I know that all the files have the same known projection.

Anyway, Using these two tools I have my indexes built and my LAS points linked into the project with the data using the specialised point cloud index. Now I want to start working with it. The first thing I need to do is derive the real LAS X, Y and Z values using VectorValue function:

--SQL9

VectorValue(GeomCoordXYZ([Geom], 0), 0) AS [X]

VectorValue(GeomCoordXYZ([Geom], 0), 1) AS [Y]

VectorValue(GeomCoordXYZ([Geom], 0), 2) AS [Z]

Here I am coming unstuck as the linked data schema is read only, so I can’t do what I originally planned which was to create computed columns for these values. My next thought that I might be able to use a query in the project along the lines of:

--SQL9

TABLE CALL TableCacheIndexGeoms(

(

SELECT VectorValue(GeomCoordXYZ([Geom], 0), 0) AS [X]

VectorValue(GeomCoordXYZ([Geom], 0), 1) AS [Y]

VectorValue(GeomCoordXYZ([Geom], 0), 2) AS [Z]

[CLASSIFICATION][Geom]

FROM [Data Source]::[LAS]

)

, TRUE

);

To derive the correct X, Y and Z values within the without any data replication. While this appears to work, takes considerable time and I think reverts back to the standard RTree index, thus loses the benefits of the specialised index.

I am therefore wondering if I am approaching this all the wrong way and ask if anyone knows of a way to make use of the specialised point cloud index while retaining access to the correct X, Y and Z values. Hopefully this is possible and if so I would be very happy to be pointed in the right direction.

Thanks in advance

adamw


10,447 post(s)
#02-Apr-20 06:50

The call to TableCacheIndexGeoms does create the standard RTREE index on the Geom field, indeed, which does lose the benefits of the specialized spatial index for point clouds in the LAS / LAS library data source. We will add support for the specialized spatial index to MAP files in the future, too.

Why do you need to create a new table though? If you want to have direct access to X / Y / Z of existing data and you do not need to alter them, you can just separate the SELECT query you have inside the call to TableCacheIndexGeoms into a query component and link that as a drawing, taking geometry from the Geom field. This will use the specialized spatial index in the LAS dataport yet you will see X / Y / Z of every point.

Is it that you want to alter the X / Y / Z values? If so, given that we also want to continue using the specialized spatial index, I guess the process might be: create a query that would select all fields from LAS, but replace the contents of the Geom field, then export the result table of that query back into a LAS and link that. If all of the geoms should be altered using the same rules, that will work. If some of the geoms should be altered and others not, that should also work. If we are talking about making various manual adjustments interactively, I guess the best would be to import LAS - edit data in Manifold (yes, using the slow index) - export LAS.

LandSystems73 post(s)
#02-Apr-20 07:55

Thanks Adam, I don't need to edit the LAS data, just extract the correct X,Y,Z values without any replication and retaining the point cloud index. I think the 'TRUE' in my TableCacheIndexGeoms was misleading so apologies for the red herring.

If you want to have direct access to X / Y / Z of existing data and you do not need to alter them, you can just separate the SELECT query you have inside the call to TableCacheIndexGeoms into a query component and link that as a drawing, taking geometry from the Geom field. This will use the specialized spatial index in the LAS dataport yet you will see X / Y / Z of every point.

Is exactly what I want to do. I will give this a whirl in the morning. Many thanks.

LandSystems73 post(s)
#05-Apr-20 00:13

Hi Adam,

I had thought I had understood what you were meaning, but perhaps not. I did a little test with 3 or 4 adjacent las files which I indexed and linked using the librarylas dataport ...

>> Clipboard-1.png

I then re-wrote my query to remove the TableCacheIndexGeoms part as shown in the screen grab below. The top pane shows the query text, the middle, the query result with the correct X, Y and Z values and the bottom shows the linked table [LAS] with the raw X, Y, Z values.

>> Clipboard-2.png

I then link my query [LAS XYZ] as a new drawing.

>> Clipboard-3.png

But the drawing remains steadfastly empty of points. I ran the query and left it for some time in case there was some form of caching going on in the background but to no avail. Hopefully I have missed something simple or misinterpreted your meaning.

My aim from this point is to be able to query the corrected X, Y, Z values further (for example to look for anomalies or to filter out ground strikes) and then to interpolate these to a grid using one of the transform templates.

Manifold 9.0.171.0 65GB RAM

Attachments:
Clipboard-1.png
Clipboard-2.png
Clipboard-3.png

adamw


10,447 post(s)
#05-Apr-20 08:25

The query doesn't need significant time to compute - the result table of the SELECT, as written, returns data on demand, so normally you should see data streaming in nearly immediately.

It seems to me you did everything correct, but there's one thing I forgot to mention - sorry for that. Unlike the table in the LAS library dataport, the table returned by the SELECT does not know the bounds of its geometry data. For a general case of a SELECT with WHERE / other clauses, computing these bounds could take significant time, so the query engine does not do that, because frequently these bounds just aren't needed. So, since the table returned by the SELECT does not know the bounds of its geometry data, when you open a drawing based on that table in a map window, the map window does not know where to zoom to see the data. The data is there and the table will return them, it is just that the window observes a region where there is no data.

To fix this, drop a different component covering the extents of your data into the window and zoom to where the data should be using that component as a guide. Eg, you can add to the map the original drawing based on the LAS library table, zoom to it, then turn it off in the Layers pane so that it doesn't interfere with the display of the drawing based on the query. Or, if the component is correctly georegistered, you can add a base layer from Bing / Google and zoom to the place on the Earth where the LAS data is. Either of the above should let you see the data in the drawing based on the query.

LandSystems73 post(s)
#06-Apr-20 00:14

Thanks for that Adam. I hadn’t appreciated this so good to know. I followed your instructions and all good as I can now see that the data really was there all the time.

In My previous post, I indicated that what I might then like to do is …

to be able to query the corrected X, Y, Z values further (for example to look for anomalies or to filter out ground strikes) and then to interpolate these to a grid using one of the transform templates.

With my sample, I decided to filter out groundstrikes and interpolate to a new grid

>> Clipboard-1

Pressing ‘Add Component’ or ‘Edit Query’ however does nothing. I am guessing that this is due to the query not knowing the bounds, therefore making it unable to build the interpolation query? If this is the case, is there any easy way round this as it is a bit of a roadblock to gridding point clouds without first having to import them?

Any pointers again much appreciated.

Attachments:
Clipboard-1.png

adamw


10,447 post(s)
#06-Apr-20 09:06

OK, I see what is going on.

You are exactly right - the culprit is that the query does not know the bounds of its geometry data. These bounds are needed to generate the transform, so if they are not known the transform cannot be generated.

There is a simple trick which you could use, however. The map window can automatically compute the bounds for layers which do not know them already. It just needs a little kick - you have to ask it to zoom to the extents of the layer explicitly. So, in your map window, right-click the layer tab for LAS XYZ Drawing and invoke Zoom. This will compute the bounds of the layer by reading all the geometry data inside of it - running a progress dialog if the amount of geometry data is big. After that, you should be able to use the Interpolate transforms as well as other transforms that use layer bounds. We will likely alter the Transform pane to ask the map window to compute the layer bounds if they are needed for the transform, so that it all just works automatically without any tricks required.

One more thing right away - after you create an interpolated image, it will likely look either all-black or all-white, due to the default style of 0 = black and 255 = white. To fix this, go to the Style pane, select R / G / B channels in the channel list, right-click one of them and invoke one of the Autocontrast commands, or Full Range if you want. This will set the values for black / white to the ones that make more sense for the image. We will perhaps do this as part of the transform in the future as well.

adamw


10,447 post(s)
#06-Apr-20 10:06

One more thing regarding using a drawing based on a query - you have to set the coordinate system of that drawing to that of the original data. There are multiple ways to do this. If you already created the drawing, the easiest is to copy / paste the coordinate system from one drawing to another - open the original drawing (from the LAS library dataport), open the Contents pane, go to the Component pane inside the Contents pane, click the button next to the coordinate system readout, invoke Copy, open the drawing based on the query, go to the Component pane again, click the button next to the coordinate system readout, invoke Assign Initial Coordinate System - Paste (or Repair Initial Coordinate System - Paste, if you already set the coordinate system to something but aren't sure you set it correctly).

I guess you may have already done something like this, because it seems that your drawing based on a query lines up with the original drawing, but I thought I'd spell it out explicitly just to be sure.

adamw


10,447 post(s)
#06-Apr-20 10:14

And one more thing.

Since you are interested in the scaled X / Y / Z values, it might be worth mentioning that the raw X / Y / Z values are already exposed in the original table, but they need scaling. Figuring out the scale factor isn't difficult - you just open the drawing and click an arbitrary point, then look at the X / Y / Z values reported in the Values tab (the values stored in the file) and at the same values reported in the Coordinates tab, and divide them in pairs. In many files, the scale factor is simply 100. I don't remember off-hand if there's any offset - if there is, we need two points to figure out what it is. What this saves you is using a SELECT query and living with the fact that it doesn't know its bounds, etc - you can just use the original table and use the values from the X / Y / Z fields exposed by the dataport, applying necessary scaling during computations.

This may, however, break with multiple files inside a library - if they are scaled differently, then scales for different points in the common table will be different. But usually, scale factors between different files in the same data set are also the same.

LandSystems73 post(s)
#07-Apr-20 21:09

Thank you so much for taking the time to lay this all out so clearly for me. It is very much appreciated. Your last paragraph ...

This may, however, break with multiple files inside a library - if they are scaled differently, then scales for different points in the common table will be different.

Might be one reason why this wouldn't be a good idea, but wouldn't adding an option to the dataport to automatically scale X, Y and Z make much of this go away?

adamw


10,447 post(s)
#08-Apr-20 08:26

Yes, we'll consider it (perhaps an option to add fields for scaled X / Y / Z so that you still have access to the raw X / Y / Z values stored within the file - or just no option and always report both scaled and raw values).

LandSystems73 post(s)
#08-Apr-20 10:10

Any or either would be great. Thanks for considering it as an option.

radiowebst47 post(s)
#27-Apr-20 22:35

is there some way to just hover and see the Z value with the coordinates like M8 does?

LandSystems73 post(s)
#27-Apr-20 23:54

If you ALT-Click on the image, the Record pane will open. Click the Pixels Tab and and value of the pixel clicked will be the one which is marked.

Manifold User Community Use Agreement Copyright (C) 2007-2021 Manifold Software Limited. All rights reserved.