Análise e visualização de dados LiDAR – Parte 5

Para o nosso exemplo de hoje vamos dar um zoom e encontrar um edifício para analisar (se você aumenta o zoom, os edifícios se tornam clicáveis).

Podemos ver que o edifício #81719 (identificador da chave primária do banco de dados) tem um número limitado de atributos, incluindo uma elevação de 1.446,43! Nossos dados LIDAR variam cerca de 450 metros, de modo que a elevação nos edifícios é provavelmente medida em pés.

Como podemos calcular a elevação para o edifício #81719 utilizando a tabela LIDAR? Veja como a lógica funciona visualmente:

1. Comece com o edifício.

2. Localize todos os patches que se cruzam naquele edifício.

3. Transforme essas manchas em pontos individuais.

4. Filtre esses pontos usando o contorno do edifício.

5. Finalmente a média dos pontos de elevações dará um valor para o edifício.

Veja o procedimento acima em SQL:

-- We run a set of subqueries sequentially using
-- the "with" keyword
WITH
-- Get the one building we are interested in
building AS (
  SELECT geom FROM buildings
  WHERE buildings.gid = 81719
),
-- All the patches that intersect that building
patches AS (
  SELECT pa FROM medford
  JOIN building ON PC_Intersects(pa, geom)
),
-- All the points in that patch
pa_pts AS (
  SELECT PC_Explode(pa) AS pts FROM patches
),
-- All the points in our one building
building_pts AS (
  SELECT pts FROM pa_pts JOIN building
  ON ST_Intersects(geom, pts::geometry)
)
-- Summarize those points by elevation
SELECT
  Avg(PC_Get(pts, 'z')) AS lidar_meters
FROM building_pts;

O resultado do SQL é 441.075 metros, que em pés é 1447.1, quase exatamente o mesmo valor do arquivo de edifícios, 1.446,43!

Ok! Porém podemos acrescentar a elevação a partir da nossa imagem LIDAR para todos os nossos edifícios? Sim, mas vai demorar algum tempo de processamento. Primeiro vamos adicionar uma coluna para aceitar o valor, e então executar uma atualização na nossa tabela do banco de dados.

-- Add column for our calculate Z value
ALTER TABLE buildings ADD COLUMN z real;
-- Update into the column
UPDATE buildings SET z = elevs.z
FROM (
  -- For every building, all intersecting patches
  WITH patches AS (
    SELECT
      buildings.gid AS buildings_gid,
      medford.id AS medford_id,
      medford.pa AS pa
    FROM medford
    JOIN buildings
    ON PC_Intersects(pa, geom)
  ),
  -- Explode those patches into points, remembering
  -- which building they were associated with
  pa_pts AS (
    SELECT buildings_gid, PC_Explode(pa) AS pts FROM patches
  )
  -- Use the building associations to efficiently
  -- spatially test the points against the building footprints
  -- Summarize per building
  SELECT
    buildings_gid,
    Avg(PC_Get(pts, 'z')) AS z
  FROM pa_pts
  JOIN buildings
  ON buildings.gid = buildings_gid
  WHERE ST_Intersects(buildings.geom, pts::geometry)
  GROUP BY buildings_gid
) AS elevs
-- Join calculated elevations to original buildings table
WHERE elevs.buildings_gid = gid;

Nossa tabela está atualizada! Confira as elevações originais e calculadas (temos que converter a coluna z de metros para pés para compará-la com a coluna elevação):

SELECT z AS z_m, z*3.28084 AS z_ft, elevation AS elevation_ft
FROM buildings;

Os valores estão bem próximos, veja:

   z_m | z_ft | elevation_ft 
--------- + ------------------ + --------------- 
 438,128 | 1.437,42663560303 | 1434.43000000 
 433,556 | 1.422,4291678418 | 1413.63200000 
 435,489 | 1.428,76987573853 | 1406.25400000 
 439,244 | 1.441,09014682129 | 1422.58200000 
 433,738 | 1.423,02460105347 | 1416.93600000 
 429,648 | 1.409,60687857788 | 1403.92400000 
 437,264 | 1.434,59264585083 | 1425.84000000 
 430,607 | 1.412,75115040894 | 1404.03675300

Vamos colocar a diferença de elevação (em metros) na tabela, veja como:

-- Add column for our calculated height
ALTER TABLE buildings ADD COLUMN height real;
-- Update the building heights by subtracting elevation of
-- the nearest street from the elevation of the building
UPDATE buildings SET height = heights.height
FROM (
  WITH candidates AS (
    SELECT
      b.gid AS building_gid,
      s.gid AS street_gid,
      s.z AS streets_z,
      b.z as buildings_z
    FROM buildings b, streets s
    WHERE ST_DWithin(b.geom, s.geom, 0.001)
    ORDER BY
      building_gid,
      ST_Distance(b.geom, s.geom)
  )
  SELECT DISTINCT ON (building_gid)
    building_gid, street_gid,
    buildings_z - streets_z AS height
  FROM candidates
) AS heights
WHERE heights.building_gid = buildings.gid;

No próximo post, finalizaremos nossa série de posts sobre LiDAR, não perca!