Agora que temos dados de elevação e altura em nossos edifícios podemos fazer um indicador ainda mais atraente visualmente.

1. Na guia “Dados”, clique no botão “Atualizar tipo de recurso”, para que GeoServer esteja ciente dos novos atributos z e z_diff.

gs_reloadft

2. Dentro do diretório de dados, localize a pasta workspaces/opengeo/lidar/buildings

3. Criar um arquivo de texto com o nome height.ftl no diretório mencionado acima com o conteúdo: ${height.value}

4. É isso aí! Agora abra a camada KML no Google Earth e dê zoom até se aproximar dos edifícios para ver o resultado.

5. http://localhost:8080/geoserver/wms/kml?layers=opengeo:buildings&mode=refresh&kmscore=50&format_options=lookatbbox:bbox=-122.8808,42.3311,-122.8806,42.3313

ge_buildings

Se você explorar a visualização em 3D dos edifícios por um tempo, você pode se deparar com alguns bairros ímpares, como este.

ge_buildings_trees1

Se você perceber na imagem acima 3 edifícios parecem um pouco fora do lugar nesta área residencial, e a torre no quintal? Algo está errado.

O que estamos vendo é o efeito das árvores nos dados LIDAR. Os pulsos de laser estão refletindo em volta das árvores, bem como nas casas fazendo com que a altitude média de pontos dentro dos polígonos das edificações sejam artificialmente infladas.

Como podemos corrigir isso?

Os dados LIDAR refletidos das árvores irá gerar mais de uma reflexão de retorno para cada pulso de saída, e quanto mais longe a reflexão, maior o tempo de retorno da luz. Ao invés de gerar nossas alturas dos edifícios como a média de todos os pontos LIDAR, queremos a média apenas os últimos retornos para cada pulso; os retornos que são os mais profundos.

Uma solução simples é dar peso aos retornos, de modo que em áreas de tipos de retorno mistos os retornos mais profundas têm mais peso. Podemos fazer isso ser substituindo a nossa função AVG no cálculo LIDAR com uma média ponderada:

-- numerator
Sum(CASE WHEN PC_Get(pts,'ReturnNumber') = 1
     THEN PC_Get(pts, 'z')
     WHEN PC_Get(pts,'ReturnNumber') = 2
     THEN 10*PC_Get(pts, 'z')
     ELSE 100*PC_Get(pts, 'z') END) /
-- denominator
Sum(CASE WHEN PC_Get(pts,'ReturnNumber') = 1
    THEN 1
    WHEN PC_Get(pts,'ReturnNumber') = 2
    THEN 10
    ELSE 100 END) AS z

Agora estamos prontos para voltar a executar o nosso cálculo de elevação e nosso cálculo de altura para ver o efeito sobre a saída do final KML.

-- Update into the buildings Z 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,
    -- Use a weighted average that heavily favors
    -- multi-return pulses
    Sum(CASE WHEN PC_Get(pts,'ReturnNumber') = 1
         THEN PC_Get(pts, 'z')
         WHEN PC_Get(pts,'ReturnNumber') = 2
         THEN 10*PC_Get(pts, 'z')
         ELSE 100*PC_Get(pts, 'z') END) /
    Sum(CASE WHEN PC_Get(pts,'ReturnNumber') = 1
        THEN 1
        WHEN PC_Get(pts,'ReturnNumber') = 2
        THEN 10
        ELSE 100 END) AS z
  FROM pa_pts
  JOIN buildings
  ON buildings.gid = buildings_gid
  WHERE ST_Intersects(buildings.geom, pts::geometry)
  -- Only use the last returns in this calculation
  AND PC_Get(pts,'ReturnNumber') = PC_Get(pts,'NumberOfReturns')
  GROUP BY buildings_gid
) AS elevs
-- Join calculated elevations to original buildings table
WHERE elevs.buildings_gid = gid;

-- 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;

Agora temos uma nova visualização para os edifícios em 3D. Você pode verificar que a torre está a metade do tamanho, e as edificações foram reduzidas em muitos andares.

ge_buildings_trees2

Este tutorial é uma tradução e adaptação livre do artigo “Analyzing and Visualizing LIDAR” publicado no site da Boundless.

Fonte: Boundless