How to create a grid of LiDAR points (X, Y, Z) using Python GDAL?

Question

I'm new to python programming, and I'm just wondering if it is possible to create a regular grid with a resolution of 0.5 by 0.5 m using LiDAR points.

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

The points are arranged randomly, and some pixel is empty (NAN value), and in some pixel more than one point. Where there is more than one point, I want to get the average value. In the end, I need to save the data in TIF format or Ascii format.

I am studying the osgeo module and GDAL, but I honestly say that I do not know if the osgeo module is the best solution.

I am very happy for helping with some code that I can learn and implement,

Thanks to Advance for the help I really need .

I do not know how best to get a grid with these parameters.

Answer 1

A little late, but perhaps this answer will be useful to others, if not for you ...

I did this with Numpy and Pandas, and it is pretty fast. I used TLS data and could do it with several million data points without any problems on a decent 2009 laptop. The key is "binning" by rounding data, and then using Pandas' GroupBy methods to aggregate and calculate funds.

If you need to round to 10, you can use np.round, otherwise you can round to an arbitrary value by making a function for this, which I did by changing 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 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 did this, I saved the output as .npy and converted to tiff with an image library, and then attached to the geometry in ArcMap. There is probably a way to do this with osgeo, but I have not used it.

Hope this helps someone at least ...

Answer 2

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

 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]) 

Answer 3

Try LAStools , especially lasgrid or las2dem .