如何使用GDAL python從LiDAR點（X，Y，Z）創建網格？

我對python編程真的很陌生，我只是想知道是否可以使用LiDAR點創建一個0.5米x 0.5米分辨率的規則網格。如何使用GDAL python從LiDAR點（X，Y，Z）創建網格？

我的數據採用LAS格式（以lasfile的形式從liblas導入文件中讀取），它們具有以下格式：X，Y，Z。 X和Y是座標。

這些點是隨機定位的，某些像素是空的（NAN值），在某些像素中有多個點。如果有更多的一點，我希望得到一個平均值。最後，我需要以TIF格式或Ascii格式保存數據。

我正在學習osgeo模塊和GDAL，但我誠實地說，我不知道如果osgeo模塊是最好的解決方案。

我與一些代碼，我可以學習和貫徹的幫助真的很高興，

在此先感謝您的幫助，我真的需要。

我不知道用這些參數獲取網格的最佳方法。

來源

2012-09-25 Gianni Spear

您可以使用直方圖功能與NumPy做分級，例如：

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

來源

2012-09-25 18:06:36 reptilicus

我上面添加一些代碼，我寫了 –

在這裏，你將不得不使用雖'np.histogram2d'。 – letmaik

嘗試LAStools，特別lasgrid或las2dem。

來源

2012-09-26 20:59:38

感謝邁克，但目標是使用外面的代碼。我需要在LAStools之外計算幾個新的idex。但謝謝你支持我！ :) –

這是一個有點晚，但也許這個答案將是你對別人有用，如果不是......

我這樣做與numpy的和熊貓，這是相當快。我使用的是TLS數據，可以在沒有任何麻煩的情況下，在數百萬個數據點上做到這一點。關鍵是將數據四捨五入，然後使用Pandas的GroupBy方法進行聚合並計算平均值。

如果你需要舍入到10的冪，你可以使用np.round，否則你可以通過修改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)

當我做這個，我已保存的輸出作爲.npy並轉換爲TIFF與圖片庫，然後在ArcMap地理參考。可能有一種方法可以用osgeo做到這一點，但我沒有使用它。

希望這可以幫助別人，至少...

來源

2013-03-18 22:56:44 cossatot

如何使用GDAL python從LiDAR點（X，Y，Z）創建網格？

回答

相關問題