Exercício 04 - FIltragem DTM
Neste exercício, utilizaremos a grade da aula anterior para remover as árvores (filtragem de nuvens de pontos) utilizando o método do Block Minimum.
Já sabemos:
- Instalar bibliotecas e carregar modulos python
- Abrir e ler arquivo LAS (lake, da aula 01)
- Preencher as falhas
!pip install laspy
importar bibliotecas e abrir drive
import numpy as np
import laspy
from google.colab.patches import cv2_imshow
from google.colab import drive
drive.mount('/content/drive')
import numpy as np
import laspy
from google.colab import drive
drive.mount('/content/drive')
Abrir arquivo Lake.laz da pasta data do seu drive. Crir um Raster e preencehr as falhas.
res = 1.0 # resolução do raster em metros
# criar grade
xmin, xmax = x.min(), x.max()
ymin, ymax = y.min(), y.max()
ncols = int(np.ceil((xmax - xmin) / res))
nrows = int(np.ceil((ymax - ymin) / res))
#cria matriz de nrows × ncols e coloca NaN em todas as células.
raster = np.full((nrows, ncols), np.nan)
col = ((x - xmin) / res).astype(int) # transforma em coordenada linha coluna
row = ((ymax - y) / res).astype(int)
# Preencher cada célula com a menor elevação encontrada
for r, c, valor in zip(row, col, z):
if np.isnan(raster[r, c]) or valor < raster[r, c]:
raster[r, c] = valor
# Guardar esta grade para fazer outros exercícios
raster_original = raster.copy()
Preencher as falhas usando dilatação.
for i in range(1, nrows - 1):
for j in range(1, ncols - 1):
if np.isnan(raster[i, j]):
vizinhanca = raster[i-1:i+2, j-1:j+2]
valores_validos = vizinhanca[~np.isnan(vizinhanca)] #pega somente os valores que existem na janela 3×3.
if len(valores_validos) > 0:
raster[i, j] = valores_validos.min() # aqui escolhe o menor da vizinhanca
Filtragem
Nesta etapa, vamos barrer toda a grade.
- Para cada ponto, buscamos na vizinhança (definida por raio) a cota mínima.
- Calculamos a diferença entre o pixel central e este mínimo
- Se a diferença for maior que um limiar de cota (limiar), a célula recebe o valor NaN. Ou seja, o ponto é removido
- Caso contrário, o valor fica inalterado (o ponto é terreno)
print("Z mínimo:", np.min(z))
print("Z máximo:", np.max(z))
janela = 5
raio = janela // 2
limiar = 2.0 # diferença de altura em metros
iteracoes = 5
for k in range(iteracoes):
c=0
raster_novo = raster_terreno.copy()
for i in range(raio, nrows - raio):
for j in range(raio, ncols - raio):
if np.isnan(raster_terreno[i, j]):
continue vizinhanca = raster_terreno[ i-raio:i+raio+1, j-raio:j+raio+1 ]
if np.all(np.isnan(vizinhanca)):
continue
z_min = np.nanmin(vizinhanca) # minimo sem NaN
delta_z = raster_terreno[i, j] - z_min
if delta_z > limiar:
raster_novo[i, j] = np.nan
c+=1
print('iteracao', k, 'removidos', c)
raster_terreno = raster_novo
#################################
plt.figure(figsize=(10, 8))
plt.imshow(
raster_terreno,
cmap='terrain',
extent=[xmin, xmax, ymin, ymax],
origin='upper'
)
plt.colorbar(label='Elevação (m)')
plt.show()
Interpole nos lugares vazios
# Coordenadas dos pixels que ainda possuem dados
linhas, colunas = np.where(~np.isnan(raster_terreno))
x_validos = xmin + colunas * res
y_validos = ymax - linhas * res
z_validos = raster_terreno[linhas, colunas]
pontos_validos = np.column_stack((x_validos, y_validos))
# Interpolação
interpolado = griddata(
pontos_validos,
z_validos,
(xi, yi),
method='linear'
)
# Mantém os dados existentes e preenche somente os NaN
#raster_dtm = raster_terreno.copy()
#mask = np.isnan(raster_dtm)
#raster_dtm[mask] = interpolado[mask]
raster_dtm = interpolado.copy()
Ver em 3d
import matplotlib.pyplot as plt
fig = plt.figure(figsize=(12, 9))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(
xi,
yi,
raster_dtm,
cmap='terrain',
linewidth=0,
antialiased=True
)
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Elevação (m)')
plt.show()
ou em 2d
filtragem em Lastools
Vamos fazer o mesmo em Lastools?
Jorge Centeno: centeno@ufpr.br