1

There are a number of tools out there that can create a raster from a LiDAR point cloud by interpolating elevation, intensity or class values, but these tools generally limit you to interpolating those particular fields.

In particular I use the Python bindings for Whitebox. Unfortunately if you want to interpolate e.g., the number of returns per pulse* you have to trick the program by somehow swapping the field name for number of returns and a compatible field (like elevation).

I've tried to do this in-memory with laspy, but the memory limit is rather low. In these cases I've had to use a tool like las2txt and subsequently txt2las to swap the headers which is quite slow.

Another option is to convert the LiDAR cloud to vector points which have more general support for interpolation, but I'd guess that's even slower.

Is there a straightforward way in Python to swap LAS field headings in memory, allowing me to create as raster that is the result of interpolating any arbitrary LiDAR field? Pure Python solutions are preferred.

*Interpolating the elevation would result in a raster that represents the elevation at a given cell (i.e., a digital elevation model). Interpolating the number of returns per pulse results in a raster that represents the average number of returns per pulse in a given cell.

Improve this question
9
  • Welcome to GIS SE! Could you please clarify the following: "if you want to interpolate e.g., the number of returns"? Could you provide an example or description of your desired end product? – Aaron♦ Aug 4 '19 at 0:45
  • Yes, I just updated the OP with an example. "Interpolating" a field should result in a raster that represents the average value of that field at each cell. – Sky Jones Aug 4 '19 at 3:41
  • Be careful when trying to interpolate non-continuous phenomena like “# of returns” because in dense vegetation areas for example, you may get no points, but interpolating your values wouldn’t show that; the result would be an interpolation if neighbouring values. – Jae Aug 5 '19 at 10:57
  • Something like a “point density” solution would likely be more relevant to a “# of returns” question, but the example you’ve given of “returns per pulse” is difficult to map since all returns of any given pulse may very well span multiple output raster cell locations. – Jae Aug 5 '19 at 11:14
  • My apologies for not being more explicit, the field I intended to give as the example is 'number of returns (given pulse)' as specified by the LAS 1.2 file format. This is of interest to me since vegetation splits returns, so it can be identified by areas where that field is elevated. This is distinct from the number of pulses in a given area. I agree that interpolating can be problematic in areas with a low pulse density, but I'm limited to doing my analysis with rasterized LiDAR data – Sky Jones Aug 5 '19 at 15:29
1

You are over complicating this a bit. Certain types of data are suitable for interpolation whereas other data is more well-suited for simple binning. In the lidar realm, elevation and intensity are commonly interpolated but, other lidar attributes (class, return number, time-stamp, strip-id) are simply binned to a pre-existing grid. From an analytic standpoint this is a straight-forward intersection operation between the point cloud, its attributes and a raster or vector grid. You may have to come up with a rule set to deal with multiple points intersecting a given grid cell.

If you stop and think about some of the other attributes that one may want to grid from a las file, they are not really suitable for interpolation as they represent descriptive, nominal data characteristics. Since it does not represent continuous data, this includes the classified point descriptor. You can also use a binning approach to derive things like returns per meter^2 at any given resolution that the grid represents. I have performed something akin to a sensitivity test, across different raster resolutions, in this way.

Improve this answer
1
  • Your comments regarding interpolating discrete data make sense; interpolating class was a bad example. The main field I want to interpolate is number of returns per pulse which, as you noted, is technically discrete as well. However I'm using number of returns per pulse as a proxy for vegetation density so I'm okay with pretending it's continuous. I think the binning approach is a good idea, and it'll give a similar result. My only concern is that I'll have areas of NoData, so I'll have to interpolate in the end anyway. But it gives me a new way to look at it, thanks! – Sky Jones Aug 6 '19 at 17:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.