PNP (Prof. Centeno)

Exercício 03 - Modelos 3D

Neste exercício, utilizaremos a linguagem Python para visualizar os dados e filtrar um DTM utilizando o método do Block Minimum.

Instalar bibliotecas e carregar modulos python
Abrir e ler arquivo LAS (lake, da aula 01)
Imprimir numero de pontos; minimos e maximos
Visualizar o arquivo em 3D (aula 01)

!pip install laspy lazrs
!pip install laspy

A seguir:importar bibliotecas e abrir drive import matplotlib.pylab as plt
import numpy as np
import laspy
from google.colab.patches import cv2_imshow
from google.colab import drive
drive.mount('/content/drive')
import matplotlib.pylab as plt
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. Imprima numero de pontos; mínimos e máximos de cada coordenada. nome = '/content/drive/My Drive/data/lake.laz' las = laspy.read(nome)
x = np.asarray(las.x)
y = np.asarray(las.y)
z = np.asarray(las.z)
print("Pontos:", len(x))
print("X:", x.min(), x.max())
print("Y:", y.min(), y.max())
print("Z:", z.min(), z.max())

Agora, visualize a distribuição dos dados em XYZ (3D)
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x,y,z, c=las.z, cmap='viridis', s=1 )
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')'

criar um RASTER com estes pontos

criar um raster

definir o tamanho da célula do raster/grid. Depende da desidade dos pontos.
calcular o tamanho da grade em função do mínimo/máximo das coordendas X e Y
preencer as células com os valores do LiDAR. usar sempre a menor cota.

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

#ver o raster
plt.figure(figsize=(10, 8))
plt.imshow(
   raster,
   cmap='terrain',
   extent=[xmin, xmax, ymin, ymax],
   origin='upper'
)
plt.colorbar(label='Elevação (m)')
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Rasterização da nuvem de pontos')
plt.show()

Preencher as falhas usando dilatação.

varrer toda a grade
buscar um pixel NAN (sem valor)
varrer a vizinhança desse pxiel e buscar o menor valor válido de um vizinho, apenas (o menor)
atribuir esse valor ao central

raster= raster_original.copy()
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
plt.figure(figsize=(10, 8))
plt.imshow( raster, cmap='terrain', extent=[xmin, xmax, ymin, ymax], origin='upper' )
plt.colorbar(label='Elevação (m)')
plt.show()

Outra opção é interpolar uma superfície a patir dos pontos disponíveis. Para isto,

os pontos LiDAR são usados para calcular uma triangulação de Delaunay.
A seguir, cada célula da grade é preenchida. Para isto, é verificado em qual triângulo o centro da célula cai e
É interpolado um valor usando estes três vértices.

from scipy.interpolate import griddata
pontos = np.column_stack((x, y))
#cria grade com ncols, nlins, variando os valores entre xmin - xmax, e o mesmo para y
xi, yi = np.meshgrid(
np.linspace(xmin, xmax, ncols),
np.linspace(ymax, ymin, nrows)
)
# Calcula a interpolação para toda a grade
interpolado = griddata( #interpola en cada celula
   pontos,
   z,
   (xi, yi),
   method='linear'
)
# Recuoera os dados originais para preencher os locais onde ja existiam valores para preservar os valores originais
raster_interp = raster_original.copy()
# Preenche somente onde havia NaN
mask = np.isnan(raster_interp)
raster_interp[mask] = interpolado[mask]
#############################
plt.figure(figsize=(10, 8))
plt.imshow( raster_interp, cmap='terrain', extent=[xmin, xmax, ymin, ymax], origin='upper' )
plt.colorbar(label='Elevação (m)')
plt.show()