Nós vamos carregar os dados LIDAR em uma tabela no PostgreSQL, e depois fazer a análise espacial usando o PostGIS, por isso, vamos precisar de um banco de dados com algumas extensões habilitadas.

  • Crie um novo banco de dados chamado LIDAR
  • Habilite as extensões pointcloud, postgis, e pointcloud_postgis.

No PostgreSQL, execute os seguintes comandos:

CREATE EXTENSION postgis;
CREATE EXTENSION pointcloud;
CREATE EXTENSION pointcloud_postgis;

Feito isso, agora nós podemos carregar os dados.

PDAL , ou “biblioteca de abstração de dados de nuvem de pontos”, serve para manipular dados de nuvens de pontos espaciais. Ele inclui uma biblioteca de código (para integração em aplicações) e uma ferramenta de linha de comando. Nuvens de pontos são frequentemente capturados por sensores LIDAR mas também podem ser capturados por muitos outros tipos de sensores.

Podemos usar o PDAL através da linha de comando para ler alguns metadados do nosso arquivo LAZ:

pdal info --input 20090429_42122c8225_ld_p23.laz --schema

O parâmetro –schema nos mostra informações de todas as dimensões contidas no arquivo, que neste caso são muitas, veja:

“X”, “Y”, “Z”, “Intensity”, “ReturnNumber”, “NumberOfReturns”, “ScanDirectionFlag”, “EdgeOfFlightLine”, “Classification”, “ScanAngleRank”, “UserData”, “PointSourceId”, and “Time”.

Nós também pode ler os metadados do arquivo:

pdal info --input 20090429_42122c8225_ld_p23.laz --metadata --xml

O retorno é um pouco confuso pois o XML não vem formatado, infelizmente, mas se você ler o arquivo formatado pode encontrar muita informação interessante.

  • A referência espacial dos dados é WGS84 Geográfico
  • Os limites dos dados são (-122.8874999,42.3125), (- 122.8749998,42.325)
  • Os dados foram criados utilizando o software “TerraScan”

Nós vamos construir um “pipeline” para ler nosso arquivo LAZ, e em seguida, escrever os dados em nosso banco de dados.

pdal_flow

O “arquivo de pipeline” é um XML que descreve o processamento. Cada processo envolve o processo precedente, resultando numa estrutura “nesting dolls“, no qual o primeiro processo (o leitor) está no centro e o último (o escritor) está no lado de fora.

Logo abaixo é código do nosso arquivo de pipeline. Note que estamos usando “EPSG:4326” como sistema de referência espacial, já que é o que aprendemos com os metadados:

  • Nosso leitor é um drivers.las.reader,
  • Nosso escritor é um drivers.pgpointcloud.writer e
  • No meio, estamos aplicando um filters.chipper
  • Para mais informações sobre filtros de pipeline e leitor/escritor de PDAL veja a documentação de referência.

Copie o código abaixo para um arquivo de pipeline com o nome de laz2pg.xml:

<?xml version="1.0" encoding="utf-8"?>
<Pipeline version="1.0">
  <Writer type="drivers.pgpointcloud.writer">
    <Option name="connection">dbname=lidar user=postgres</Option>
    <Option name="table">medford</Option>
    <Option name="srid">4326</Option>
    <Filter type="filters.chipper">
      <Option name="capacity">400</Option>
      <Filter type="filters.cache">
        <Reader type="drivers.las.reader">
          <Option name="filename">20090429_42122c8225_ld_p23.laz</Option>
          <Option name="spatialreference">EPSG:4326</Option>
        </Reader>
      </Filter>
    </Filter>
  </Writer>
</Pipeline>

Agora estamos prontos para executar o processo de carregamento de dados:

pdal pipeline laz2pg.xml

Quando o processo estiver completo, haverá uma nova tabela na base de dados:

                           Table "public.medford"
 Column |    Type    |                      Modifiers
--------+------------+------------------------------------------------------
 id     | integer    | not null default nextval('medford_id_seq'::regclass)
 pa     | pcpatch(1) |
Indexes:
    "medford_pkey" PRIMARY KEY, btree (id)

Note que o tipo da coluna “pa” na base de dados, é pcpatch (1). O pcpatch refere-se ao tipo de dados, que é um conjunto de nuvem de pontos, agrupados em uma área quadrada, um “patch” de dados. O (1) refere-se ao “formato” dos pontos dentro do patch: quantas dimensões cada ponto tem, e o que essas dimensões são.

SELECT * FROM pointcloud_formats WHERE pcid = 1;

A informação é difícil de ler no formato de saída do PostgreSQL, mas com uma versão formatada pode ser mais fácil de entender.

Podemos usar nosso conhecimento do banco e as funções da extensão pointcloud para saber mais sobre os nossos dados, veja:

-- How many points are in our cloud? (11418635)
SELECT Sum(PC_NumPoints(pa))
FROM medford;
-- What is the average elevation of the first patch? (439.384)
WITH pts AS (
  SELECT PC_Explode(pa) AS pt
  FROM medford LIMIT 1
)
SELECT Avg(PC_Get(pt,'Z')) FROM pts;
-- What does the first point look like?
WITH pts AS (
  SELECT PC_Explode(pa) AS pt
  FROM medford LIMIT 1
)
SELECT PC_AsText(pt) FROM pts LIMIT 1;
-- {
--  "pcid":1,
--  "pt":[-122.887,42.3125,439.384,42,1,1,1,0,1,6,181,343,419629,1.14073e+07,0]
-- }
-- How many patches do we have? (28547)
SELECT Count(*)
FROM medford;
-- What is the min/max elevation in our cloud? (421.20/467.413)
SELECT
  Min(PC_PatchMin(pa, 'z')) AS min,
  Max(PC_PatchMax(pa, 'z')) AS max
FROM medford;
-- What does a patch look like as a geometry?
SELECT st_asewkt(pa::geometry) FROM medford LIMIT 1;
-- SRID=4326;POLYGON((-122.8874998 42.3125002,-122.8874998 42.312613,
-- -122.8874414 42.312613,-122.8874414 42.3125002,-122.8874998 42.3125002))
-- What does a point look like as a geometry?
WITH pts AS (
  SELECT PC_Explode(pa) AS pt
  FROM medford LIMIT 1
)
SELECT ST_AsEWKT(pt::geometry) FROM pts LIMIT 1;
-- SRID=4326;POINT(-122.8874601 42.3125002 439.384)

Há mais informações sobre a extensão pointcloud disponíveis na página de documentação da extensão.

No próximo post veremos como disponibilizar os dados LiDAR em um mapa.