私はPythonプログラミングが初めてで、LiDARポイントを使用して0.5 x o.5 mの解像度の通常のグリッドを作成できるかどうか疑問に思っていました.
私のデータは LAS 形式 (liblas インポート ファイルから lasfile として読み取る) であり、次の形式を持っています: X、Y、Z。X と Y は座標です。
ポイントはランダムに配置され、一部のピクセルは空 (NAN 値) であり、一部のピクセルには複数の 1 つのポイントがあります。1 点が多い場合は、平均値を取得したいと考えています。最後に、データを TIF 形式または Ascii 形式で保存する必要があります。
私は osgeo モジュールと GDAL を研究していますが、osgeo モジュールが最適なソリューションであるかどうかはわかりません。
勉強して実装できるコードを手伝ってくれて本当にうれしいです。
助けてくれてありがとう、私は本当に必要です。
これらのパラメーターを使用してグリッドを取得する最良の方法がわかりません。
少し遅れていますが、この回答は、あなたにとってではないにしても、他の人にとって役立つかもしれません...
私は Numpy と Pandas でこれを行いましたが、かなり高速です。私は TLS データを使用していましたが、まともな 2009 年製のラップトップで何の問題もなく、数百万のデータ ポイントでこれを行うことができました。重要なのは、データを丸めて「ビニング」し、Pandas の GroupBy メソッドを使用して集計を行い、平均を計算することです。
10 の累乗に丸める必要がある場合は、np.round を使用できます。それ以外の場合は、関数を作成して任意の値に丸めることができます。これは、この 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でそれを行う方法がありますが、私はそれを使用していません。
これが少なくとも誰かに役立つことを願っています...
たとえば、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])