Subscribe to this thread
Home - General / All posts - Grid maths in M9 - Transfering lidar texture to a coarser DEM
danb


1,642 post(s)
#30-Jul-18 03:23

I have a task which I really want to undertake in Manifold 9 due to the volume of data involved. Unfortunately however it will involve using grid maths and conditional logic to manipulate and combine grids (as would be achieved via Surface > Transform in Mfd8).

I am attempting to apply the surface ‘texture’ of more detailed, but geographically restricted lidar data to a coarser but geographically extensive DEM. This might sound like an odd thing to wish to do, but the additional information encoded in the lidar ‘texture’ is very useful for fine guidance of covariates used in digital soil mapping of low lying areas.

To this end, I have created a simple proof of concept in Mfd8 which I have outlined below, but ideally would like to complete the process in Mfd9 due to the volume of lidar data rather than processing tile by tile in Mfd8.

Is this sort of processing currently available in Mfd9 given that while the two grids are expressed in the same projection, the pixel dimensions and offsets are different?

BTW Merge in Mfd9 won’t work for me in this case due to height differences between the DEM and Lidar data creating a large height step along the join line.

Mfd8 Processing steps (please feel free to make any suggestions for improvements as this is really just a proof of concept):

  1. Harvest the lidar texture (from [LIDAR_DEM])

This I have done as simply as possible, by generating a highly smoothed 'average' lidar surface using the median pixel value over an arbitrary wide window here 5 pixels deep (11x11 window). I subtract this 'average' surface from the unmodified lidar to retain the texture without the heights.

[LIDAR_DEM] - MedValue([LIDAR_DEM], 5)

The resulting grid is called [LIDAR_TEXTURE].

  1. Combine the lidar texture and the coarse DEM data ([BASE_DEM])

Combination is achieved using conditional statements. If the lidar texture is null (here -9999), use the coarser base dataelse use the coarse DEM plus the lidar texture. Because the range of values in the texture is small, this doesn't really change the overall height of the base DEM, rather only imposes the local texture which is more important than the height.Combination by this method will also not create a significant step in heights across the lidar/coarse DEM boundary.

IIf([LIDAR_TEXTURE] = -9999, [BASE_DEM], [BASE_DEM] + [LIDAR_TEXTURE])

The resulting grid is called [BASE RESULT].

  1. The result of this conditional statement does leave some pixel artefacts along the join line between the datasets (not sure why). I have smoothed these using something like:

[BASE RESULT FINAL]:

IIf([BASE RESULT]>5000, MedValue([BASE RESULT], 4), [BASE RESULT])

Attachments:
Clipboard-1.jpg


Landsystems Ltd ... Know your land | www.landsystems.co.nz

adamw

8,037 post(s)
#01-Aug-18 13:46

It all looks fairly simple except for the first step with the median filter which is doable, but more complex than it needs to be because we don't have the median filter as a built-in function (coming).

Step 1: Do median filter on the more detailed but covering a smaller area LIDAR DEM. This currently requires a script. Tell me if you want it. Write the result into the second tile field in the same table as the original image.

Step 2: Subtract the result of the median filter from the LIDAR DEM, leaving only the texture. This is just a Subtract transform between the tile fields. If we were writing the result of the median filter into a separate image, this would have required a join between the two images on X and Y.

Step 3: Fill invisible pixels in the texture with 0. This is me getting rid of IIf - let's just convert it into an addition. The filling of invisible pixels is done using TileCombine(tile, TileMakeNew(cx, cy, 0)), invisible pixels in the first tile get filled from the second tile.

Step 4: Project the texture to the exact coordinate system of the less detailed but covering a bigger area base DEM, using the exact same tile size as the base DEM. The coordinate system, including offsets and scales, and the tile size have to be exactly the same, this will make sure the X and Y values in the base DEM and the projected texture correspond to the same tiles.

Step 5: Add two images together. This requires a join on X and Y. Adding tiles is done using just tile1+tile2.

Step 6: Smooth the result using the median filter again.

We are planning to make it easier to combine images dynamically automatically projecting them during transforms, this should make a fair bit of the above automatic.

danb


1,642 post(s)
#04-Aug-18 00:34

Apologies for the delay in responding. Many thanks for this Adam, I suspected that the infrastructure would be largely there but not the GUI implementation. Very glad also to hear that it is coming as I believe M9 will ultimately transform much of my work with lidar and point clouds.

With regards your kind offer for a median filter script, if a builtin filter is coming, I don't wish to distract you, though I am interested if it will demonstrate a way to generate custom filters further down the track.


Landsystems Ltd ... Know your land | www.landsystems.co.nz

adamw

8,037 post(s)
#04-Aug-18 08:57

Let's wait until the built-in filter then (yes, it is coming). We will provide examples of custom filters regardless.

tjhb

8,168 post(s)
#13-Aug-18 02:54

I'm going to have a decent go at writing a median filter in the next week or so. I need to learn to write Manifold 9 filters in a standardised way for my own purposes.

I would be really grateful if Dan could test the particular filter, perhaps against Manifold 8; and if Adam could cast an eye over the eventual code and suggest areas for improvement (including as to speed and efficient data access).

danb


1,642 post(s)
#13-Aug-18 05:05

That would be great Tim and a really useful learning exercise for me. I look forward to testing your filter in due course.


Landsystems Ltd ... Know your land | www.landsystems.co.nz

adamw

8,037 post(s)
#13-Aug-18 07:29

Of course I'll cast an eye, I love things like this.

tjhb

8,168 post(s)
#13-Aug-18 08:45

Thank you both. This will be fun.

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