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