Neste post descrevo, e submeto a opiniões, os passos adotados para um
processamento semi-automático de
vetorização a partir de imagem de satélite Sentinel-2
tendo em vista a vetorização de áreas de matas
com suficiente distinção entre natural=wood e landuse=forest (mata natural e mata cultivada).
Importante:
esta metodologia é aqui proposta como um auxilío destinado somente ao mapeamento de coberturas de terreno (landcover).
O processo todo resulta em simplificação, como curvas com espaçamento entre nós não menor que 10m (a resolução da imagem é 10m/px).
Não serve para objetos pequenos como lotes, praças, etc., pois removeria detalhes que podem ser vistos nas imagens padrão do OSM.
O propósito é:
-poder gerar desenho de grandes áreas de mata densa, em interiores do território, não densamente urbanizados;
-com diferenciação de vegetação usando índices apropriados, para forest e wood, o que não é facilmente, e/ou comumente, distinguido nos desenhos sobre as imagens padrão.
Mesmo assim, deve ser verificado visualmente o resultado ao final do processamento, junto com as imagens padrão do OSM.
DESCRIÇÃO INICIAL DOS DADOS:
IMAGENS: SENTINEL-2
Cobertura: ~100x100Km
Documentação original:
https://www.sentinel-hub.com/develop/documentation/eo_products/Sentinel2EOproducts
Licença compatível com OSM:
conforme https://wiki.openstreetmap.org/wiki/Sentinel-2
Fonte das imagens:
https://earthexplorer.usgs.gov/ camadas Sentinel-2
Bandas utilizadas: B11
(testadas todas as bandas, B11 se mostrou a melhor para o objetivo; condição:nuvem=0)
Índices utilizados: NDVI; CRE
NDVI : Normalized Difference Vegetation Index
Fórmula: NDVI = (B08 - B04) / (B08 + B04)
CRE : Chlorophyll Red-Edge (abbrv. Chlred-edge)
Fórmula: CRE = (B07 / B05) ^ -1
PROCESSAMENTO - PASSOS:
1)ESCOLHA DA IMAGEM DE SATÉLITE - SENTINEL:
Escolhida imagem em data que não gere discrepâncias (para mais ou para menos) de atividade vegetal:
preferencialmente entre equinócio e solstício; evitar alto inverno e alto verão.
Caso Bom Jesus, RS:
https://earthexplorer.usgs.gov/ - Sentinel - "T22JEP_20180420T132231_B11"
Data: 2018/04/20
Hora: 13:22 (09:22 Brasil; gera alguma sombra)
(A imagem abaixo pode ser aberta em nova guia para visualizar detalhes em zoom maior.)
2)Examinar no QGIS valores do raster; procurar valores limites de distinções de classes:
-Usar a imagem Bing para localizar e marcar exemplos claramente identificáveis dos tipos de vegetação a serem distinguidos em classes.
Exige conhecimento da vegetação típica básica, e sobretudo poder identificá-la em imagem de satélite em cores naturais, como Bing.
Além de evidentemente poder distinguir os demais elementos, como área urbanizada, campo ralo, trilhas e corpos de água.
A identificação dos tipos wood e forest é possível de fazer somente sobre imagem Bing. Basta saber distingui-los.
O que o presente método propõe é usar uma seleção amostral, feita sobre a imagem Bing, de objetos dos quais se tenha seguro conhecimento, que sejam representantes destes 2 tipos, wood e forest, bem como dos demais objetos que não são estas 2 classe, para, nas imagens Sentinel, avaliar os valores típicos de pixels dos objetos selecionados, testando se os mesmos limites de valores são encontrados em matas típicas em algumas diversas porções da imagem.
Os valores típicos em uma imagem podem apresentar alguma variação em outras imagens e outros locais. Por isso todos os passos do método deve ser repetidos em trabalhos com outras imagens ou territórios.
Uma vez que se encontre os valores limites típicos na imagem, usá-los como parâmetros para deduzir os objetos dos tipos wood e forest na imagem toda.
-Áreas de mata cultivada isoladas são mais fáceis de distinguir no Bing. Não tanto quando plantadas em meio a mata natural.
Assim é preferível procurar primeiramente áreas destacadas de plantação, em campo aberto. Como lotes mais ou menos regulares de mata plantada.
As espécies plantadas mais comuns nesta região de enfoque são sobretudo o Pinus, para indústria de celulose, e menos comumente Eucalipto, mais para construção e lenga.
As florestas de Pinus, e florestas plantadas em geral, costumam apresentar espécimes regularmente espaçados, formando uma malha regular densa e identificável. Procurar por estas feições.
São plantadas em campos planos e também em encostas.
Nas encostas podem ser detectadas por várias trilhas de serviço para extração, escalonadas.
Eucalipto costuma ser plantado ao redor de sedes de fazendas em campo aberto.
A superfície apresenta maior homogeneidade de tamanho de indivíduos. Forma uma superfície mais regular.
-As matas nativas ou mais velhas, com menor atividade biológica, e por isso índices diversos nas imagens, são o restante de matas a identificar. Em geral apresentam a superfície da cobertura mais irregular. E em locais onde não haja estradas ou trilhas de serviço indicando acesso constante.
Estas áreas exclusivas com respectivos objetos-tipo servirão de amostragem de elementos das classes para avaliação dos valores de pixels.
Podem ser feitos polígonos para marca-los, em layer à parte.
Com estes objetos-tipo, serão verificados e anotados os valores típicos dos pixels nas respectivas imagens de referência, examinadas alternando os layers sobre o mesmo objeto amostral: B11, NDVI, CRE. E anotados os limites de valores em relação aos objetos das demais classes.
Estes valores limites típicos serão usados para distinção automática na imagem toda.
CASO: Matas nos municípios de Bom Jesus e Monte Alegre dos Campos, RS:
SENSOR / ÍNDICE :: B11 :: NDVI :: CRE
Min. :: 511 :: -0.87 :: 0.2
Max. :: 2709 :: +0.87 :: 0.75
MATA NATIVA (n.=wood) :: 1200-1800::0.30 /0.80 ::
MATA PLANTADA (l.=forest) ::150-1200 :: 0.70/0.80 ::
Mata na sombra(encosta.Sul) ::150-400 ::0.30/0.65 :: 0.4 - 0.5
Mata no sol (encosta.N.) :: 1200-1800 :: 0.65/0.80:: 0.25-0.3
campo ralo (null) :: 1500 -3000 :: :
urbanização :: 1700-3000 ::-0.10/0.20 :
água parada (açude) :: 50-150 ::-0.10/0.10 :: 0.8-1.0
água corrente (rio) :: 150-300 :: -0.30/0.60 :: 0.5-0.8
LIMITE P/ ÁGUA = 0 :: 0.5
Resultado observado:
B11 distingue bem entre mata wood e forest.
Em diferentes áreas do terreno abrangido na imagem, as mesmas classes de matas apresentam alguma variação, o que pode gerar ainda algum grau de imprecisão. No entanto o limite de matas x campos ou urbanização é bem distinto. Não separa rios de forest; para isto necessita contrastar B11 com outros índices mais adequados, a seguir.
-NDVI e CRE destacam o que não é vegetação, como corpos de água.
[IMAGEM-2]
3)CORREÇÕES:
Corrigir o que não é bem detectado ou misturado quando somente usando na B11:
mata na sombra (encostas Sul); água corrente
Método: usados os índices NDVI e CRE para contrastar B11
OPERAÇÃO: B11+(1-NDVI)+CRE
Raster calculator - SINTAXE: "T22JEP_20180420T132231_B11@1" + ((1-"NDVI@1")*1800) + ("CRE@1" *2600)
Agregados empiricamente os fatores ponderadores NDVI*1800 e CRE*2600, para que fiquem acima de B11 (~2000) e equilibrá-los.
RESULTADO: "B11+1-NDVI+CRE.tif" : (10m/px)
[IMAGEM-3]
4)TESTAR CLASSES VISUALMENTE:
-Igual à etapa (2), para B11+1-NDVI+CRE:
Resultado:
5 CLASSES :
(1)forest <= 2150 <
(2)wood <= 2600 <
(3)campo alto ou ciliar <= 3500 <
(4)campo ralo ou trilha <= 3950 <
(5)água, estrada, cidade, lavourado
5)PROCESSAR A IMAGEM, REDUZINDO ÀS CLASSES TESTADAS:
CONDICIONAIS "if/then/else":
Raster calculator:
SINTAXE: (("@1"<x)* operação ) => se sim=1(opera); se não=0(não opera)
SINTAXE :
(1* ( "B11+1-NDVI+CRE@1" <=2150 ))
+ (2* (( 2150 < "B11+1-NDVI+CRE@1" ) AND ( "B11+1-NDVI+CRE@1" <= 2600 ) ))
+ (3* (( 2600 < "B11+1-NDVI+CRE@1" ) AND ( "B11+1-NDVI+CRE@1" <= 3500 ) ))
+ (4* (( 3500 < "B11+1-NDVI+CRE@1" ) AND ( "B11+1-NDVI+CRE@1" <= 3950 ) ))
+ (5* ( "B11+1-NDVI+CRE@1" > 3950 ))
Resultado:
Imagem com valores limitados às classes: "B11-5classes.tif"
[IMAGEM-4]
6)SIMPLIFICAR ÁREAS, REMOVER PIXELS ISOLADOS:
OPERAÇÕES: ~1min cada
QGIS: Processing / Toolbox /
* GDAl / Analysis / Sieve
* GRASS / r.neighbor
1 Sieve T3 C8 - Remove até 3px(10x30m)
2 Neig Max 3 Ci - Aumenta borda do maior, 1px cada lado (10m+10m)
3 Sieve T3 C4
4 Neig Mode 3 Sq
5 Sieve T6 C8 - Remove até 6px(20x30m)
Resultado:
5 classes filtradas: "B11-5classes-Filtrado.tif" (10m/px)
[IMAGEM-5]
7)REDUZIR A IMAGEM SOMENTE ÀS CLASSES DE INTERESSE (forest e wood; 1 e 2) E IGUALAR DEMAIS CLASSES (3,4,5) A "ZERO":
QGIS: Raster calculator: (( "@1" > 2 ) * 0) + (( "@1" <=2 ) * "@1" )
8)CRUZAR O RASTER PROCESSADO COM OS VETORES EXISTENTES NO OSM, DE RIOS E ESTRADAS, PARA MELHORAR RECORTES NAS MATAS:
-Baixar dados do OSM/overpass: highway=* or waterway=* and type:way
-Criar um novo layer vetor SHP contendo as linhas de highway e waterway que servirão para recortar.
-Adicionar um campo de número real, com valor "0"(zero) para todas as linhas/polígonos a serem filtrados
-Gerar um buffer sobre todas as linhas:
o buffer deve ser maior que a resolução da imagem destino e menor que os menores grupos de pixels após as filtragens.
Para não gerar pixels ligados só nos cantos nem elminar a mais do que as filtragens anteriores.
Utilizado buffer de 20m; suficiente para imagem de 10m/px. Contendo valor classe=0.
-Salvar o buffer como SHP
RASTERIZAR O VETOR SOBRE A IMAGEM:
-Fazer uma cópia simples do raster anteriormente processado, a ser alterada
-Adicionar o layer vetor do buffer sobre a imagem destino:
QGIS: Raster / Conversion / Rasterize: SHP com o buffer; campo p/ valor "0"(zero); imagem destino
ÚLTIMA FILTRAGEM:
simplificar; remover últimos pequenos clumps de até 20 pixels isolados (4x5 pixels) que sobram em alguns recortes :
Sieve T20 C4
Resultado:
"B11-0-1-2-Filtrado+OSM.tif"
Não é objetivo do processo captar áreas com menos de 40x40m nas distinções, o que equivale a pouco, o que gera muitas imprecisões.
Somente áreas de matas com maior homogeneidade, que possuem dimensões muito maiores.
[IMAGEM-6]
9) ELIMINAR DEMAIS CLASSES NÃO UTILIZADAS - Igualar "zero" a "null":
(ajuda a reduzir a quantidade de polígnos no poligonize):
QGIS: gdal_translate / -a_nodata <X>
(X = valor da classe p/null)
Resultado:
Raster Classes 1+2+null : "B11-1-2-null+OSM.tif"
10)VETORIZAR - Gerar polígonos do raster:
QGIS: Raster / conversion / poligonize
MANTER A MESMA CRS NO PROJETO
11)SUAVIZAR POLÍGONOS - despixelizar:
QGIS: GRASS commands -> v.generalize.smooth:
(https://grass.osgeo.org/grass64/manuals/v.generalize.html)
* Sliding Averaging:
look_ahead = 3 (mín=3, sempre ímpar)
slide = 0.7
Resultado: ~5min
"Sliding-Average-L3-SL07-OK.shp"
LEGENDA:
Verde claro: mata nativa (natural=wood)
Verde escuro: mata cultivada (landuse=forest)
Branco: null
12)ADICIONAR TAGS OSM E LIMPAR:
ELIMINAR DO SHP POLÍGONOS TRUNCADOS DE BORDA DE IMAGEM ou
RESTRINGIR POLÍGONOS A MUNICÍPIO
Salvar como novo layer: "B11+OSM-limpo.shp"
-Selecionar todo polígono que encosta nas bordas do layer e remover manualmente (pode fazer um buffer menor que os limites do layer);
-Adicionar tags (Table Manager + Field Calculator): landuse=forest; natural=wood;
-Remover tags sem uso (area; classe).
Salvar como WGS84-EPSG:4326
13) Opcional: Selecionar por máscara de município
(pois o JOSM fica muito pesado se abrir mais de 10.000 polígonos)
QGIS: Vector / Spatial Query:
"SHP-origem" + Intersects + "SHP-máscara" - selection / invert / remove
PROCESSOS NO JOSM:
1)ABRIR CADA SHP SEPARADAMENTE :
-Sem baixar dados do OSM
-salvar como .osm
Já entra como polígono ou multipol(outer/inner); curvas com ~1nó/10m (mesma resolução das imagens raster)
Monte Alegre dos Campos:
QGIS: SHP 5,2MB; features/polígonos= 3.519
JOSM: "B11-MA-4326.osm" : 23,8MB; nodes= 243.070
2)SIMPLIFICAR:
Select: type:way / Symplify way
Resultado: nodes= 136.853
3) VALIDAR SEPARADAMENTE
-Sem baixar dados do OSM
Resultado:
Errors: (0)
Warnings: (28)
Self-intersecting polygon ring (3)
Self-intersecting ways (25)
Todos foram casos de auto-intersecção em cruz (quina) - imagem abaixo.
Talvez possa evitar previamente no QGIS com
QGIS: GRASS commands -> v.generalize.smooth / Displacement
(https://grass.osgeo.org/grass64/manuals/v.generalize.html)
Solução:
Alteração manual: "unglue(G)" nos nós de cada intersecção; afastar nós e ways
Cuidar: testar todos os nós da intersecção após o unglue; manter polígonos e multipolígonos fechados ao final
VALIDAR NOVAMENTE:
Resultado: VALIDAÇÃO COMPLETA (0)
4) TESTES DE CONFLITO COM O EXISTENTE:
Baixar dados do OSM;
VALIDAR:
Resultado: VALIDAÇÃO COMPLETA (0)
Pronto para aprovação e upload
5) PARA UPLOAD DOS NOVOS DADOS AO OSM:
Proposto salvar em changeset exclusivo, com a nota:
"adição de polígonos de matas (wood;forest) semi-automatizados com validação manual"
O arquivo .osm resultante deste teste para verificação pode ser baixado neste repositório
(não feito upload ao OSM).
Comentários:
Onde acho que ainda pode melhorar:
Na filtragem, manter o mesmo buffer de aumento e retorno de limites de áreas vizinhas.
Ficou um pouco a mais, 10m maior, pra dentro da mata. Por um lado evita invadir outras coisas.
Mas usando os vetores do OSM para buffer de estradas e waterway, isso já resolve.
A imagem classificada tem 10m/px.
A área cobre 100.000mx100.000m.
Na minha opinião tem que simplificar e eliminar clusters pequenos, como <40x40m (~16px), pois:
-o OSM usa simplificação de áreas de terreno, exemplo maior as coast_lines;
-manter áreas pequenas de <40x40m, tumultua o resultado com objetos pequenos sem alta precisão para a imagem, objetos que não são objetivo do processo.
O objetivo são as áreas de mata densa e mais homogêneas.
Por favor, fiquem à vontade para tecer comentários se desejarem,
como se vale ou não vale a pena investir neste tipo de processo para desenhar áreas de mata densa,
em regiões de interior, para o OSM.
OpenStreetMap Blogs
.
