Baixando dados LiDAR de MDT/MDS do GeoSampa em lote!

Jun 3, 2022 17 min read geo, LiDAR, MDT, MDS, gdal, GeoSampa, shell
Dados de MDT no visualizador LAZ integrado do GeoSampa

Documentações já disponíveis

Além da documentação oficial do GeoSampa existem tutoriais aqui na Internet ensinando como baixar os dados LiDAR de elevação do município de São Paulo. Aqui, o principal foco deste post é mostrar uma forma de baixar esses dados em lote.

Alguns links que abordam os dados sobre GeoSampa:

Passos para baixar uma quadrícula

Primeiro precisa ir ao site GeoSampa, e lá vai aparecer um mapa com a delimitação da cidade de São Paulo (se não aparecer clica na primeira opção no controle de camadas que fica ao lado direito da tela no Mapa Base e em Político-Administrativo) depois desce quase no final deste mesmo menu do lado direito, clica em Articulação de Imagens e depois seleciona MDT/MDS 2017 e a partir deste momento vai aparecer um grid em cima do mapa base. Dando ampliação no lugar de interesse, vai ver uns números nas quadrículas formados por 7 dígitos. Vamos entender este numero como o ID daquela quadrícula, ele é muito importante para poder baixar os dados pertencentes a esta quadricula.

Este grid cobrindo todo o município de São Paulo é formado por 5362 quadriculas no total. Vamos supor que a sua área de estudo está na quadrícula 2342-363. Para baixar os dados LiDAR disponíveis no site você precisa ir para a caixa de ferramenta do lado esquerdo, e clica em cima da ferramenta de Pesquisar (está com um ícone de lupa), seleciona a aba Download imagens/MDC e em Tipo seleciona o produto desejado MDS ou MDT 2017 e digita o ID da quadrícula anotado anteriormente (no caso, 2342-363) e clica ok. Vai acontecer uma ampliação automático naquela quadrícula e ao clicar em cima dela, aparecem opções de fazer download do arquivo informado mas também já tem um visualizador web integrado com Potree, muito interessante. Aliás, paralelamente, na caixa de ferramenta tem uma ferramenta LiDAR 3D que permite selecionar uma área maior para esta visualização.

Bom, para baixar um arquivo ou até 5, é tranquilo, é só repetir este processo, você anota o ID da quadrícula e vai entrando e baixando uma por uma. Se por acaso, você deu sorte e a sua área for bem maior do que isso (isso foi o meu caso, precisei da área toda) aí vem o problema, são 5362 quadriculas, fica complicado ir baixando as quadriculas uma por uma. Então, este post dá ênfase nesta parte.

Passos para baixar em lote

Primeiro, precisamos entender um pouco do link que baixa os arquivos da base do GeoSampa. Segue um exemplo para os dois tipos de dados, MDT e MDS na quadrícula 2342-363:

Não sei se notou que em algum lugar no meio, está o número da quadrícula desejada, conseguiu ver isso? Muito bem! Então, conhecendo o ID da quadrícula desejada, podemos alterar somente este valor e vamos conseguir baixar as outras quadrículas fazendo parte da nossa área de estudo. Isso já facilita um pouco se formos comparar com a opção de ficar fazendo o mesmo caminho clicando até baixar cada quadrícula.

Assim, podemos dividir o link em 3 partes: antes do ID(parte 1) + o ID + depois do ID(parte 3). A parte antes do ID e depois dele permanecem intactos. Então, encontrando uma forma de ir alterando o ID, podemos criar uma lista de links com todas as quadrículas. Mas o problema é que não temos a sequência dos 5362 IDs. Mas o bacana, é que nós conseguimos baixar esta camada de quadrícula no site GeoSampa mesmo. Para isso, vai em Download de arquivos >> Download de Camadas >> Articulação de Imagens >> Articulação MDT, aí dá a opção de baixar em KMZ ou em Shapefile, clica em Shapefile e vai dar a opção de baixar as quadriculas em SAD69 ou SIRGAS 2000, eu escolhi SIRGAS 2000 aí vai baixar um arquivo SIRGAS_SHP_quadriculamdt.zip contendo, obviamente, o shapefile. Baixei este arquivo para minha pasta temp.

Vendo as 3 partes dos links:

Produto Parte_1 ID Parte_3
MDT ht\tp://download.geosampa.prefeitura.sp.gov.br/PaginasPublicas/downloadArquivo.aspx?orig=DownloadMapaArticulacao&arq=MDT_2017%5C 2342-363 .zip&arqTipo=MAPA_ARTICULACAO
MDS ht\tp://download.geosampa.prefeitura.sp.gov.br/PaginasPublicas/downloadArquivo.aspx?orig=DownloadMDSLAZ&arq=MDS_2017%5C 2342-363 .laz&arqTipo=LAZ

Baixando com wget no bash

O wget é uma ferramenta principalmente de linha de comando que permite baixar dados na Internet passando como argumento o link do arquivo desejado. Por exemplo, se quisermos baixar os dados LiDAR da quadrícula 2342-363 com wget podemos simplesmente chamar wget link_do_arquivo e vai baixar o arquivo 2342-363.zip para o nosso disco local. A mesma ideia continua, se eu alterar os 7 dígitos do ID, e este ID for válido, consigo ir baixando outras quadrículas. Só que, estando na linha de comando podemos fazer muito mais do que isso! O wget consegue receber uma lista de links e vai iterando em cima da lista e vai baixando. Então, bastar criar uma lista de links, não é?!

Para isso, precisamos de uma coisa primeiro, é a lista dos IDs. Tendo a lista dos IDs é só juntar as três partes. Assumo que você já baixou aquele shapefile mencionado acima, ele vem compactado.

Descompactando o SIRGAS_SHP_quadriculamdt.zip

Você pode descompactar com o seu descompactador favorito, mas aqui estando no meu terminal com interpretador bash, vou utilizar a ferramenta unzip aqui mesmo.

Onde estou?
1

pwd

>>> /home/tredgi/temp

O que temos dentro desta pasta?

1

ls

>>> SIRGAS_SHP_quadriculamdt.zip
>>> xlimite

A pasta xlimite será utilizada logo

Extraindo com unzip
1

unzip SIRGAS_SHP_quadriculamdt

>>> Archive:  SIRGAS_SHP_quadriculamdt.zip
>>>    creating: SIRGAS_SHP_quadriculamdt/
>>>   inflating: SIRGAS_SHP_quadriculamdt/SIRGAS_SHP_quadriculamdt.shx  
>>>   inflating: SIRGAS_SHP_quadriculamdt/SIRGAS_SHP_quadriculamdt.shp  
>>>   inflating: SIRGAS_SHP_quadriculamdt/SIRGAS_SHP_quadriculamdt.prj  
>>>   inflating: SIRGAS_SHP_quadriculamdt/SIRGAS_SHP_quadriculamdt.dbf  
>>>  extracting: SIRGAS_SHP_quadriculamdt/SIRGAS_SHP_quadriculamdt.cpg
Depois da extração, temos esta pasta a mais
1

ls

>>> SIRGAS_SHP_quadriculamdt
>>> SIRGAS_SHP_quadriculamdt.zip
>>> xlimite
O que tem na pasta?
1
2

# entrando na pasta
ls ./SIRGAS_SHP_quadriculamdt

>>> SIRGAS_SHP_quadriculamdt.cpg
>>> SIRGAS_SHP_quadriculamdt.dbf
>>> SIRGAS_SHP_quadriculamdt.prj
>>> SIRGAS_SHP_quadriculamdt.shp
>>> SIRGAS_SHP_quadriculamdt.shx

Você pode abrir e verificar no seu programa de SIG favorito, inclusive fazer uma superposição de sua área de estudo para captar os IDs envolvidos. Como uma ajuda rápida, segue abaixo uma implementação do grid em cima de um mapa base de OpenStreetMap. Passando o mouse no centroide de uma quadrícula é revelado o seu ID.

Visualizador interativo com os IDs

GDAL no Bash

Para abrir e usar dados do shapefile vamos utilizar uma ferramenta chamada GDAL (Geospatial Data Abstraction Library), ela é uma biblioteca de código aberto para formatos de dados geoespaciais rasters e vetoriais. Se como eu, você estiver utilizando uma distribuição Linux baseada no Debian, basta rodar no seu terminal sudo apt-get install gdal-bin, caso não, vai encontrar no site da ferramenta a forma correta de fazer a instalação para o sistema que está utilizando.

Para listar informações sobre dados vetoriais se utiliza o comando ogrinfo acompanhado de uma ou mais opções/flags, como: -so (Summary Only) e -al (List all features of all layers). Vamos ver um resumo do shapefile com -so.

1

ogrinfo -so ./SIRGAS_SHP_quadriculamdt/SIRGAS_SHP_quadriculamdt.shp

>>> INFO: Open of `./SIRGAS_SHP_quadriculamdt/SIRGAS_SHP_quadriculamdt.shp'
>>>       using driver `ESRI Shapefile' successful.
>>> 1: SIRGAS_SHP_quadriculamdt (Polygon)

Podemos ver que o arquivo tem uma única camada chamada SIRGAS_SHP_quadriculamdt, vamos ver um resumo sobre ela

1

ogrinfo -so ./SIRGAS_SHP_quadriculamdt/SIRGAS_SHP_quadriculamdt.shp SIRGAS_SHP_quadriculamdt

>>> INFO: Open of `./SIRGAS_SHP_quadriculamdt/SIRGAS_SHP_quadriculamdt.shp'
>>>       using driver `ESRI Shapefile' successful.
>>> 
>>> Layer name: SIRGAS_SHP_quadriculamdt
>>> Metadata:
>>>   DBF_DATE_LAST_UPDATE=2020-02-18
>>> Geometry: Polygon
>>> Feature Count: 5362
>>> Extent: (313179.846000, 7343480.630333) - (360683.191000, 7416432.613000)
>>> Layer SRS WKT:
>>> PROJCRS["SIRGAS 2000 / UTM zone 23S",
>>>     BASEGEOGCRS["SIRGAS 2000",
>>>         DATUM["Sistema de Referencia Geocentrico para las AmericaS 2000",
>>>             ELLIPSOID["GRS 1980",6378137,298.257222101,
>>>                 LENGTHUNIT["metre",1]]],
>>>         PRIMEM["Greenwich",0,
>>>             ANGLEUNIT["degree",0.0174532925199433]],
>>>         ID["EPSG",4674]],
>>>     CONVERSION["UTM zone 23S",
>>>         METHOD["Transverse Mercator",
>>>             ID["EPSG",9807]],
>>>         PARAMETER["Latitude of natural origin",0,
>>>             ANGLEUNIT["degree",0.0174532925199433],
>>>             ID["EPSG",8801]],
>>>         PARAMETER["Longitude of natural origin",-45,
>>>             ANGLEUNIT["degree",0.0174532925199433],
>>>             ID["EPSG",8802]],
>>>         PARAMETER["Scale factor at natural origin",0.9996,
>>>             SCALEUNIT["unity",1],
>>>             ID["EPSG",8805]],
>>>         PARAMETER["False easting",500000,
>>>             LENGTHUNIT["metre",1],
>>>             ID["EPSG",8806]],
>>>         PARAMETER["False northing",10000000,
>>>             LENGTHUNIT["metre",1],
>>>             ID["EPSG",8807]]],
>>>     CS[Cartesian,2],
>>>         AXIS["(E)",east,
>>>             ORDER[1],
>>>             LENGTHUNIT["metre",1]],
>>>         AXIS["(N)",north,
>>>             ORDER[2],
>>>             LENGTHUNIT["metre",1]],
>>>     USAGE[
>>>         SCOPE["Engineering survey, topographic mapping."],
>>>         AREA["Brazil - between 48°W and 42°W, northern and southern hemispheres, onshore and offshore."],
>>>         BBOX[-33.5,-48,5.13,-42]],
>>>     ID["EPSG",31983]]
>>> Data axis to CRS axis mapping: 1,2
>>> qmdt_cod: String (20.0)

Neste resumo, podemos ver que a camada tem 5362 feições (features) de tipo polígono, e o mais importante, o nome da coluna contendo essas feições é qmdt_cod. Com a mesma ferramenta ogrinfo podemos listar estas feições com a flag -al. Vamos tentar, mas lembra que são 5362? Pois é, não dá para exibir tudo isso aqui, vou pegar somente o finalzinho com o comando tail.

1

ogrinfo -al ./SIRGAS_SHP_quadriculamdt/SIRGAS_SHP_quadriculamdt.shp | tail -n 12

>>> OGRFeature(SIRGAS_SHP_quadriculamdt):5359
>>>   qmdt_cod (String) = 3336-444
>>>   POLYGON ((331725.429 7373051.662,332256.363 7373057.815,332263.037 7372481.037,331732.124 7372474.883,331725.429 7373051.662))
>>> 
>>> OGRFeature(SIRGAS_SHP_quadriculamdt):5360
>>>   qmdt_cod (String) = 3336-442
>>>   POLYGON ((331718.736 7373628.441,332249.6915 7373634.592,332256.363 7373057.815,331725.429 7373051.662,331718.736 7373628.441))
>>> 
>>> OGRFeature(SIRGAS_SHP_quadriculamdt):5361
>>>   qmdt_cod (String) = 3336-414
>>>   POLYGON ((331712.043 7374205.219,332243.02 7374211.369,332249.6915 7373634.592,331718.736 7373628.441,331712.043 7374205.219))

São as três ultimas feições, note que a última é apresentada como a feição 5361, é porque começou no zero. Também, precisei colocar para apresentar 12 linhas porque as informações de uma feição fica em várias linhas, daí a opção -n 12.

Filtrando com grep

Como você pode ver, o campo qmdt_cod está entregando para nós uma coisa muito importante, os 7 dígitos do id da quadrícula, então se pudéssemos filtrar somente os dígitos, seria ótimo! Podemos chamar outra utilidade dentro do bash que é o grep. Vamos fazer um pipe seguido pelo comando grep filtrando somente este padrão com as tags -Eo (extended-regexp + only-matching) e para não mostrar tudo aqui, outro pipe pegando somente os 10 primeiros id.

A ajuda/documentação dos comandos pode ser acessada pelo terminal mesmo, digitando o_comando + --help (ex: ogrinfo --help) ou caso queira mais detalhes pode a acessar o manual deles sem sair do shell com man (ex: ogrinfo echo) 
1

ogrinfo -al ./SIRGAS_SHP_quadriculamdt/SIRGAS_SHP_quadriculamdt.shp | grep -Eo "[0-9]{4}-[0-9]{3}" | head -n 10

>>> 3212-353
>>> 3212-344
>>> 3212-343
>>> 3211-464
>>> 3211-463
>>> 3211-454
>>> 3211-453
>>> 3211-444
>>> 2242-323
>>> 3335-313

Aqui com [0-9]{4}-[0-9]{3} estou informando que quero somente o padrão dos 7 dígitos, 4 dígitos seguido de um traço e depois 3 dígitos. Para mostrar todos os IDs, não precisa colocar o pipe e o comando head.

Se está acompanhando comigo, deve perceber que acabamos de encontrar a forma de gerar os 5362 IDs. Podemos contar eles aqui mesmo com wc -l, antes de ir mais longe:

1

ogrinfo -al ./SIRGAS_SHP_quadriculamdt/SIRGAS_SHP_quadriculamdt.shp | grep -Eo "[0-9]{4}-[0-9]{3}" | wc -l

>>> 5362

Certo?! É sempre bom verificar. Tendo os IDs podemos gerar todos os 5362 links. Vamos ver como fazer isso:

Juntando as partes

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

# guardando as duas partes dentro de variáveis
parte1="http://download.geosampa.prefeitura.sp.gov.br/PaginasPublicas/downloadArquivo.aspx?orig=DownloadMapaArticulacao&arq=MDT_2017%5C"
parte3=".zip&arqTipo=MAPA_ARTICULACAO"
# guardando a a lista dos id, dentro de uma estrutura de lista, no bash
quadriculas=($(ogrinfo -al ./SIRGAS_SHP_quadriculamdt/SIRGAS_SHP_quadriculamdt.shp | grep -Eo "[0-9]{4}-[0-9]{3}" | head -n 10))
# iterando com um for-loop em cima da lista dos id, e ir juntando com as duas partes
for i in ${quadriculas[@]}
do
echo $parte1$i$parte3
done

>>> http://download.geosampa.prefeitura.sp.gov.br/PaginasPublicas/downloadArquivo.aspx?orig=DownloadMapaArticulacao&arq=MDT_2017%5C3212-353.zip&arqTipo=MAPA_ARTICULACAO
>>> http://download.geosampa.prefeitura.sp.gov.br/PaginasPublicas/downloadArquivo.aspx?orig=DownloadMapaArticulacao&arq=MDT_2017%5C3212-344.zip&arqTipo=MAPA_ARTICULACAO
>>> http://download.geosampa.prefeitura.sp.gov.br/PaginasPublicas/downloadArquivo.aspx?orig=DownloadMapaArticulacao&arq=MDT_2017%5C3212-343.zip&arqTipo=MAPA_ARTICULACAO
>>> http://download.geosampa.prefeitura.sp.gov.br/PaginasPublicas/downloadArquivo.aspx?orig=DownloadMapaArticulacao&arq=MDT_2017%5C3211-464.zip&arqTipo=MAPA_ARTICULACAO
>>> http://download.geosampa.prefeitura.sp.gov.br/PaginasPublicas/downloadArquivo.aspx?orig=DownloadMapaArticulacao&arq=MDT_2017%5C3211-463.zip&arqTipo=MAPA_ARTICULACAO
>>> http://download.geosampa.prefeitura.sp.gov.br/PaginasPublicas/downloadArquivo.aspx?orig=DownloadMapaArticulacao&arq=MDT_2017%5C3211-454.zip&arqTipo=MAPA_ARTICULACAO
>>> http://download.geosampa.prefeitura.sp.gov.br/PaginasPublicas/downloadArquivo.aspx?orig=DownloadMapaArticulacao&arq=MDT_2017%5C3211-453.zip&arqTipo=MAPA_ARTICULACAO
>>> http://download.geosampa.prefeitura.sp.gov.br/PaginasPublicas/downloadArquivo.aspx?orig=DownloadMapaArticulacao&arq=MDT_2017%5C3211-444.zip&arqTipo=MAPA_ARTICULACAO
>>> http://download.geosampa.prefeitura.sp.gov.br/PaginasPublicas/downloadArquivo.aspx?orig=DownloadMapaArticulacao&arq=MDT_2017%5C2242-323.zip&arqTipo=MAPA_ARTICULACAO
>>> http://download.geosampa.prefeitura.sp.gov.br/PaginasPublicas/downloadArquivo.aspx?orig=DownloadMapaArticulacao&arq=MDT_2017%5C3335-313.zip&arqTipo=MAPA_ARTICULACAO

Baixando de fato

Acabamos de criar 10 links para baixar 10 quadriculas de MDT! Note que eu limitei a 10 para não exibir 5 mil linhas aqui. Tendo uma lista basta integrá-los no wget e pronto. Bastaria trocar o comando echo no for-loop com wget e ao invés de exibir os links no terminal, o wget iria baixá-los para a nossa maquina. Ainda, mais do que isso, podemos baixar vários arquivos em simultâneo. É bem provável que você tenha vários núcleos no seu processador, então podemos rodar o comando em paralelo com a ferramenta parallel (os nomes são bem óbvios, não é? rs) assim distribuir as tarefas entre os núcleos. Vamos ver como fica.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

# guardando as duas partes dentro de variáveis
parte1="http://download.geosampa.prefeitura.sp.gov.br/PaginasPublicas/downloadArquivo.aspx?orig=DownloadMapaArticulacao&arq=MDT_2017%5C"
parte3=".zip&arqTipo=MAPA_ARTICULACAO"
# guardando a a lista dos id, dentro de uma estrutura de lista, no bash
quadriculas=($(ogrinfo -al ./SIRGAS_SHP_quadriculamdt/SIRGAS_SHP_quadriculamdt.shp | grep -Eo "[0-9]{4}-[0-9]{3}" | head -n 2))
# iterando com um for-loop em cima da lista dos id, e ir juntando com as duas partes
for i in ${quadriculas[@]}
do
echo $parte1$i$parte3
done | parallel "wget --content-disposition {}"

O parallel por padrão, vai rodar 10 jobs simultaneamente, mas você pode especificar com a opção -j n onde n é a quantidade de jobs a serem executados simultaneamente. Note que neste caso aqui, baixando dados, a velocidade da sua Internet é fator limitante (não adianta colocar para baixar 100 arquivos por vez, não tendo uma Internet que aguenta isso, pois simplesmente não vai mais rápido rsrs).

1

ls

>>> 3212-344.zip
>>> 3212-353.zip
>>> SIRGAS_SHP_quadriculamdt
>>> SIRGAS_SHP_quadriculamdt.zip
>>> xlimite

Você pode descompactar os arquivos e terá os arquivos com extensão LAZ para trabalhar. Aqui podemos fazer isso tranquilamente com

1

unzip "*[0-9].zip"

>>> Archive:  3212-353.zip
>>>  extracting: MDT_3212-353_1000_v0.laz  
>>> 
>>> Archive:  3212-344.zip
>>>  extracting: MDT_3212-344_1000_v0.laz  
>>> 
>>> 2 archives were successfully processed.

Por ter outros arquivos .zip na mesma pasta, precisei restringir para somente aqueles cujos nomes contem um número antes da extensão .zip. Senão, eu poderia chamar diretamente unzip "*.zip".

1

ls

>>> 3212-344.zip
>>> 3212-353.zip
>>> MDT_3212-344_1000_v0.laz
>>> MDT_3212-353_1000_v0.laz
>>> SIRGAS_SHP_quadriculamdt
>>> SIRGAS_SHP_quadriculamdt.zip
>>> xlimite

Filtros mais específicos

Por ser uma demonstração, eu estou limitando o download a somente a dois arquivos, mas isso fica a seu critério baixar os arquivos desejados, inclusive todos eles de uma vez. Aqui tem a vantagem, de poder selecionar conforme um padrão qualquer filtrando com o grep. Vamos supor que você queira os arquivos cujo primeiro digito seja 2 ou 4, o segundo digito 3, o quarto 2, o quinto 2 e o sétimo 1, é só informar isso no padrão do grep e pronto, fica assim.

1

ogrinfo -al ./SIRGAS_SHP_quadriculamdt/SIRGAS_SHP_quadriculamdt.shp | grep -Eo "[2|4]3[0-9]2-2[0-9]1"

>>> 4312-231
>>> 4312-221
>>> 4312-211
>>> 2322-231
>>> 4312-261
>>> 4312-251
>>> 4312-241
>>> 2322-261
>>> 2342-231
>>> 2342-221
>>> 2342-261
>>> 2342-251
>>> 2342-241

Tendo isso, basta colocá-lo na repetição do for-loop!

Eu sei que aqui, seria muito melhor rodar filtros espaciais, pois nada garante que essas quadriculas sejam vizinhos ou tenham qualquer outra relação espacial. A boa notícia é que, dá para fazer isto, com o GDAL, de uma olhadinha na documentação.

Seleção baseada em relação espacial

Até o momento, estamos filtrando as quadrículas pelo padrão dos números mesmo. Esta forma não garanta nenhuma relação espacial entre as mesma (pelo menos, eu não encontrei como escolher duas quadriculas que sejam vizinhas uma da outra de uma forma segura). Você deve estar pensando, mas se são dados com referenciais espaciais, porque não faz uma seleção pelo atributo espacial? E sim, você está certx! Os dados espaciais são especiais, e devemos tratá-los como tal.

Na realidade, é muito provável que você tenha um polígono delimitando a sua área e queira baixar somente as quadrículas cuja área interceptam a área do seu limite. Podemos fazer isso aqui também com o GDAL.

O GDAL aceita comandos SQL, que é uma forma muito poderosa de fazer manipulações e consultas em dados estruturados como tabelas. O GDAL aceita dialetos de SQL como por exemplo, o SQLite que por sua vez estendido pelo SpatiaLite dando a ele a capacidade de trabalhar com dados contendo atributos espaciais.

Eu tenho outro arquivo shapefile com nome de limite contendo um polígono delimitando uma área que estou interessado, gostaria de selecionar somente as quadrículas que interceptam a área deste meu polígono. Usarei o dialeto SQLite para fazer isso no GDAL, enquanto o SpatiaLite é chamado em segundo plano.

1
2
3
4

ogrinfo -dialect sqlite -sql \
"SELECT * FROM SIRGAS_SHP_quadriculamdt AS quadriculas, './xlimite/limite.shp'.limite AS limite WHERE ST_Intersects(quadriculas.geometry, limite.geometry)" \
./SIRGAS_SHP_quadriculamdt/SIRGAS_SHP_quadriculamdt.shp
 # o backslash serve para quebrar linhas compridas!

>>> INFO: Open of `./SIRGAS_SHP_quadriculamdt/SIRGAS_SHP_quadriculamdt.shp'
>>>       using driver `ESRI Shapefile' successful.
>>> 
>>> Layer name: SELECT
>>> Geometry (GEOMETRY): Polygon
>>> Geometry (GEOMETRY): Polygon
>>> Feature Count: 8
>>> Extent (GEOMETRY): (322446.413000, 7392554.977000) - (324587.080000, 7393734.207000)
>>> Extent (GEOMETRY): (322677.886173, 7392807.102998) - (324387.224103, 7393552.486424)
>>> SRS WKT (GEOMETRY):
>>> PROJCRS["SIRGAS 2000 / UTM zone 23S",
>>>     BASEGEOGCRS["SIRGAS 2000",
>>>         DATUM["Sistema de Referencia Geocentrico para las AmericaS 2000",
>>>             ELLIPSOID["GRS 1980",6378137,298.257222101,
>>>                 LENGTHUNIT["metre",1]]],
>>>         PRIMEM["Greenwich",0,
>>>             ANGLEUNIT["degree",0.0174532925199433]],
>>>         ID["EPSG",4674]],
>>>     CONVERSION["UTM zone 23S",
>>>         METHOD["Transverse Mercator",
>>>             ID["EPSG",9807]],
>>>         PARAMETER["Latitude of natural origin",0,
>>>             ANGLEUNIT["degree",0.0174532925199433],
>>>             ID["EPSG",8801]],
>>>         PARAMETER["Longitude of natural origin",-45,
>>>             ANGLEUNIT["degree",0.0174532925199433],
>>>             ID["EPSG",8802]],
>>>         PARAMETER["Scale factor at natural origin",0.9996,
>>>             SCALEUNIT["unity",1],
>>>             ID["EPSG",8805]],
>>>         PARAMETER["False easting",500000,
>>>             LENGTHUNIT["metre",1],
>>>             ID["EPSG",8806]],
>>>         PARAMETER["False northing",10000000,
>>>             LENGTHUNIT["metre",1],
>>>             ID["EPSG",8807]]],
>>>     CS[Cartesian,2],
>>>         AXIS["(E)",east,
>>>             ORDER[1],
>>>             LENGTHUNIT["metre",1]],
>>>         AXIS["(N)",north,
>>>             ORDER[2],
>>>             LENGTHUNIT["metre",1]],
>>>     USAGE[
>>>         SCOPE["Engineering survey, topographic mapping."],
>>>         AREA["Brazil - between 48°W and 42°W, northern and southern hemispheres, onshore and offshore."],
>>>         BBOX[-33.5,-48,5.13,-42]],
>>>     ID["EPSG",31983]]
>>> Data axis to CRS axis mapping: 1,2
>>> SRS WKT (GEOMETRY):
>>> PROJCRS["SIRGAS 2000 / UTM zone 23S",
>>>     BASEGEOGCRS["SIRGAS 2000",
>>>         DATUM["Sistema de Referencia Geocentrico para las AmericaS 2000",
>>>             ELLIPSOID["GRS 1980",6378137,298.257222101,
>>>                 LENGTHUNIT["metre",1]]],
>>>         PRIMEM["Greenwich",0,
>>>             ANGLEUNIT["degree",0.0174532925199433]],
>>>         ID["EPSG",4674]],
>>>     CONVERSION["UTM zone 23S",
>>>         METHOD["Transverse Mercator",
>>>             ID["EPSG",9807]],
>>>         PARAMETER["Latitude of natural origin",0,
>>>             ANGLEUNIT["degree",0.0174532925199433],
>>>             ID["EPSG",8801]],
>>>         PARAMETER["Longitude of natural origin",-45,
>>>             ANGLEUNIT["degree",0.0174532925199433],
>>>             ID["EPSG",8802]],
>>>         PARAMETER["Scale factor at natural origin",0.9996,
>>>             SCALEUNIT["unity",1],
>>>             ID["EPSG",8805]],
>>>         PARAMETER["False easting",500000,
>>>             LENGTHUNIT["metre",1],
>>>             ID["EPSG",8806]],
>>>         PARAMETER["False northing",10000000,
>>>             LENGTHUNIT["metre",1],
>>>             ID["EPSG",8807]]],
>>>     CS[Cartesian,2],
>>>         AXIS["(E)",east,
>>>             ORDER[1],
>>>             LENGTHUNIT["metre",1]],
>>>         AXIS["(N)",north,
>>>             ORDER[2],
>>>             LENGTHUNIT["metre",1]],
>>>     USAGE[
>>>         SCOPE["Engineering survey, topographic mapping."],
>>>         AREA["Brazil - between 48°W and 42°W, northern and southern hemispheres, onshore and offshore."],
>>>         BBOX[-33.5,-48,5.13,-42]],
>>>     ID["EPSG",31983]]
>>> Data axis to CRS axis mapping: 1,2
>>> Geometry Column 1 = GEOMETRY
>>> Geometry Column 2 = GEOMETRY
>>> qmdt_cod: String (0.0)
>>> id: Integer64 (0.0)
>>> OGRFeature(SELECT):0
>>>   qmdt_cod (String) = 3313-164
>>>   id (Integer64) = 1
>>>   GEOMETRY = POLYGON ((324048.476 7393151.04,324041.537 7393727.821,324573.243 7393734.207,324580.161 7393157.427,324048.476 7393151.04))
>>>   GEOMETRY = POLYGON ((324269.531983566 7392936.00389127,322694.699332767 7392807.10299815,322677.886172794 7393429.18991715,323159.863425352 7393552.48642362,323653.049451225 7393518.86010367,324387.224103377 7393535.67326364,324269.531983566 7392936.00389127))
>>> 
>>> OGRFeature(SELECT):1
>>>   qmdt_cod (String) = 3313-163
>>>   id (Integer64) = 1
>>>   GEOMETRY = POLYGON ((323516.79 7393144.633,323509.83 7393721.415,324041.537 7393727.821,324048.476 7393151.04,323516.79 7393144.633))
>>>   GEOMETRY = POLYGON ((324269.531983566 7392936.00389127,322694.699332767 7392807.10299815,322677.886172794 7393429.18991715,323159.863425352 7393552.48642362,323653.049451225 7393518.86010367,324387.224103377 7393535.67326364,324269.531983566 7392936.00389127))
>>> 
>>> OGRFeature(SELECT):2
>>>   qmdt_cod (String) = 3313-154
>>>   id (Integer64) = 1
>>>   GEOMETRY = POLYGON ((322985.103 7393138.207,322978.122 7393714.99,323509.83 7393721.415,323516.79 7393144.633,322985.103 7393138.207))
>>>   GEOMETRY = POLYGON ((324269.531983566 7392936.00389127,322694.699332767 7392807.10299815,322677.886172794 7393429.18991715,323159.863425352 7393552.48642362,323653.049451225 7393518.86010367,324387.224103377 7393535.67326364,324269.531983566 7392936.00389127))
>>> 
>>> OGRFeature(SELECT):3
>>>   qmdt_cod (String) = 3313-153
>>>   id (Integer64) = 1
>>>   GEOMETRY = POLYGON ((322453.415 7393131.761,322446.413 7393708.546,322978.122 7393714.99,322985.103 7393138.207,322453.415 7393131.761))
>>>   GEOMETRY = POLYGON ((324269.531983566 7392936.00389127,322694.699332767 7392807.10299815,322677.886172794 7393429.18991715,323159.863425352 7393552.48642362,323653.049451225 7393518.86010367,324387.224103377 7393535.67326364,324269.531983566 7392936.00389127))
>>> 
>>> OGRFeature(SELECT):4
>>>   qmdt_cod (String) = 3313-332
>>>   id (Integer64) = 1
>>>   GEOMETRY = POLYGON ((324055.416 7392574.258,324048.476 7393151.04,324580.161 7393157.427,324587.08 7392580.647,324055.416 7392574.258))
>>>   GEOMETRY = POLYGON ((324269.531983566 7392936.00389127,322694.699332767 7392807.10299815,322677.886172794 7393429.18991715,323159.863425352 7393552.48642362,323653.049451225 7393518.86010367,324387.224103377 7393535.67326364,324269.531983566 7392936.00389127))
>>> 
>>> OGRFeature(SELECT):5
>>>   qmdt_cod (String) = 3313-331
>>>   id (Integer64) = 1
>>>   GEOMETRY = POLYGON ((323523.751 7392567.85,323516.79 7393144.633,324048.476 7393151.04,324055.416 7392574.258,323523.751 7392567.85))
>>>   GEOMETRY = POLYGON ((324269.531983566 7392936.00389127,322694.699332767 7392807.10299815,322677.886172794 7393429.18991715,323159.863425352 7393552.48642362,323653.049451225 7393518.86010367,324387.224103377 7393535.67326364,324269.531983566 7392936.00389127))
>>> 
>>> OGRFeature(SELECT):6
>>>   qmdt_cod (String) = 3313-322
>>>   id (Integer64) = 1
>>>   GEOMETRY = POLYGON ((322992.085 7392561.423,322985.103 7393138.207,323516.79 7393144.633,323523.751 7392567.85,322992.085 7392561.423))
>>>   GEOMETRY = POLYGON ((324269.531983566 7392936.00389127,322694.699332767 7392807.10299815,322677.886172794 7393429.18991715,323159.863425352 7393552.48642362,323653.049451225 7393518.86010367,324387.224103377 7393535.67326364,324269.531983566 7392936.00389127))
>>> 
>>> OGRFeature(SELECT):7
>>>   qmdt_cod (String) = 3313-321
>>>   id (Integer64) = 1
>>>   GEOMETRY = POLYGON ((322460.418 7392554.977,322453.415 7393131.761,322985.103 7393138.207,322992.085 7392561.423,322460.418 7392554.977))
>>>   GEOMETRY = POLYGON ((324269.531983566 7392936.00389127,322694.699332767 7392807.10299815,322677.886172794 7393429.18991715,323159.863425352 7393552.48642362,323653.049451225 7393518.86010367,324387.224103377 7393535.67326364,324269.531983566 7392936.00389127))

Pode ver que retornou somente 8 quadriculas, que são aquelas que atendam as condições espaciais informadas no comando acima. Agora, é só usar o grep para pegar somente os 7 dígitos para construir o seu link.

1
2
3

ogrinfo -dialect sqlite -sql \
"SELECT * FROM SIRGAS_SHP_quadriculamdt AS quadriculas, './xlimite/limite.shp'.limite AS limite WHERE ST_Intersects(quadriculas.geometry, limite.geometry)" \
./SIRGAS_SHP_quadriculamdt/SIRGAS_SHP_quadriculamdt.shp | grep -Eo "[0-9]{4}-[0-9]{3}"

>>> 3313-164
>>> 3313-163
>>> 3313-154
>>> 3313-153
>>> 3313-332
>>> 3313-331
>>> 3313-322
>>> 3313-321

Agora sei que essas 8 quadrículas são vizinhas e cobrem a minha área de interesse!

Conclusão

Seguindo a mesma lógica, você pode fazer isso na linguagem/ferramenta que preferir, pois, os passos consistem em criar uma lista dos IDs e depois juntá-los com as duas outras partes do link. Depois integrar isso numa iteração com alguma ferramenta que consegue baixar eles de forma automática para a sua máquina local.

No próximo post, mostrarei uma rotina de como sair dos arquivos LAZ para GeoTIFF. Não é sempre o caso, mas muitas vezes precisamos os dados no formato tif para integrar com outros dados para rodar alguma análise específica. Ah! Caso não se sentiu muito confortável com a linha de comando e o bash em específico, dá uma olhada neste outro post , onde eu entrei em alguns detalhes do fundamento desta maravilhosa ferramenta.

Até a próxima!

LiDAR parallel wget gdal MDT MDS bash comand-line free data open data
Réginal Exavier
Réginal Exavier

Related