How to create a grid from LiDAR points (X,Y,Z) with GDAL python?

Question

I'm new really to python programming, and I was just wondering if you can create a regular grid of 0.5 by o.5 m of resolution using LiDAR points.

My data are in LAS format (reading with from liblas import file as lasfile) and they have the following format: X,Y,Z. Where X and Y are coordinates.

The points are randomly positioned and some pixel are empty (NAN value) and in some pixel there are more of one points. Where there are more of one point, I wish to obtain a mean value. In the end i need to save the data in a TIF format or Ascii format.

I am studying osgeo module and GDAL but I honest to say that i don't know if osgeo module is the best solution.

I am really glad for help with some code that i can study and implement,

Thanks in Advance for the help, I really need.

I don't know the best way to get a grid with these parameters.

Hanukkah

Answer 1

It's a bit late but maybe this answer will be useful for others, if not for you...

I have done this with Numpy and Pandas, and it's pretty fast. I was using TLS data and could do this with several million data points without any trouble on a decent 2009-vintage laptop. The key is 'binning' by rounding the data, and then using Pandas' GroupBy methods to do the aggregating and calculate the means.

If you need to round to a power of 10 you can use np.round, otherwise you can round to an arbitrary value by making a function to do so, which I have done by modifying this SO answer.

import numpy as np
import pandas as pd

# make rounding function:
def round_to_val(a, round_val):
    return np.round( np.array(a, dtype=float) / round_val) * round_val

# load data
data = np.load( 'shape of ndata, 3')
n_d = data.shape[0]

# round the data
d_round = np.empty( [n_d, 5] )
d_round[:,0] = data[:,0]
d_round[:,1] = data[:,1]
d_round[:,2] = data[:,2]

del data  # free up some RAM

d_round[:,3] = round_to_val( d_round[:,0], 0.5)
d_round[:,4] = round_to_val( d_round[:,1], 0.5)

# sorting data
ind = np.lexsort( (d_round[:,4], d_round[:,3]) )
d_sort = d_round[ind]

# making dataframes and grouping stuff
df_cols = ['x', 'y', 'z', 'x_round', 'y_round']
df = pd.DataFrame( d_sort)
df.columns = df_cols
df_round = df[['x_round', 'y_round', 'z']]
group_xy = df_round.groupby(['x_round', 'y_round'])

# calculating the mean, write to csv, which saves the file with:
# [x_round, y_round, z_mean] columns.  You can exit Python and then start up 
# later to clear memory if that's an issue.
group_mean = group_xy.mean()
group_mean.to_csv('your_binned_data.csv')

# Restarting...
import numpy as np
from scipy.interpolate import griddata

binned_data = np.loadtxt('your_binned_data.csv', skiprows=1, delimiter=',')
x_bins = binned_data[:,0]
y_bins = binned_data[:,1]
z_vals = binned_data[:,2]

pts = np.array( [x_bins, y_bins])
pts = pts.T

# make grid (with borders rounded to 0.5...)
xmax, xmin = 640000.5, 637000
ymax, ymin = 6070000.5, 6067000

grid_x, grid_y = np.mgrid[640000.5:637000:0.5, 6067000.5:6070000:0.5]

# interpolate onto grid
data_grid = griddata(pts, z_vals, (grid_x, grid_y), method='cubic')

# save to ascii
np.savetxt('data_grid.txt', data_grid)

When I've done this, I have saved the output as a .npy and converted to a tiff with the Image library, and then georeferenced in ArcMap. There is probably a way to do that with osgeo but I haven't used it.

Hope this helps someone at least...

Callant

Answer 2

You can use the histogram function in Numpy to do binning, for instance:

import numpy as np
points = np.random.random(1000)
#create 10 bins from 0 to 1
bins = np.linspace(0, 1, 10)
means = (numpy.histogram(points, bins, weights=data)[0] /
             numpy.histogram(points, bins)[0])

Amigo

Here you would have to use np.histogram2d though. – Brno

Answer 3

Try LAStools, particularly lasgrid or las2dem.

Departmentalize

Thanks Mike but the goal is using a code outside. I need to calculate several new idex outside LAStools. But Thanks for supporting me!!! :) – Hanukkah