Neste post descrevo, e submeto a opiniões, os passos adotados para um teste de processamento semi-automático de vetorização a partir de imagem de satélite Sentinel-2, com controle de parâmetros e validação manuais feitos pelo usuário, 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 apresentada como uma preparação para uma possibilidade de proposta de mapeamento para a comunidade OSM, como um metodo auxiliar ao mapeamento, e que tem como foco somente o mapeamento voltado a coberturas de terreno (landcover). Por enquanto trata-se de testes, offline. Não feito upload.
-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). Por isso não serve para objetos pequenos como lotes, praças, etc., pois não mapearia 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)
[IMAGEM-1: abaixo. 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 as imagens em cores naturais, Sentinel-TCI (True Color - ver imagem abaixo ) e 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. 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. Neste método, deve ser feita ainda preferencialmente sobre a imagem TCI do conjunto Sentinel-2, uma vez que esta apresenta o conjunto dos mesmos objetos no mesmo instante de tempo (áreas de florestas podem ter sido cortadas em diferentes períodos). Basta saber distinguir os tipos, o que é viável com algum pouco treinamento do olho.

O que o presente método propõe é usar uma seleção amostral, feita sobre as imagens em cores naturais, 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 dos demais sensores e ínidices 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 na imagem em cor natural. Não tanto quando plantadas em meio a mata natural, o que exige um olhar mais treinado aos detalhes das feições típicas (ver na imagem abaixo).
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 lenha.

As florestas de Pinus, e florestas plantadas em geral, costumam apresentar espécimes regularmente espaçados, formando uma malha regular densa, com uma superfície mais regular, como um tapete plano. 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 de florestas plantadas de modo geral costuma apresentar maior homogeneidade de tamanho de indivíduos.

-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 quanto ao tamanho dos indivíduos. Mais facilmente encontradas em locais onde não haja estradas ou trilhas de serviço indicando acesso constante, e em faixas variáveis ao longo das margens de cursos d’água, devido ao estatuto legal de proteção das matas marginais.

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 imagens de referência a serem examinadas alternando, sobre o mesmo objeto amostral, os respectivos layers: 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 de tipos de vegetais na imagem toda.

---

CASO: Matas nos municípios de Bom Jesus e Monte Alegre dos Campos, RS:
QGIS: View menu : Identify feature : testar valores de pixels individuais

---

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
valor limite para água :: < 150 :: < 0.3 :: >0.5

---

Resultado observado:
* B11 distingue bem entre os tipos de matas “wood” e “forest”.
Para isto ver também a documentação Sentinel sobre as aplicações dos diversos sensores e índices:
https://www.sentinel-hub.com/develop/documentation/eo_products/Sentinel2EOproducts
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 para uma distinção absolutamente exata. No entanto o limite de matas x campos ou urbanização é bem distinto. Por outro lado, não separa rios de forest; para isto necessita contrastar B11 com outros índices mais adequados, a seguir.
* NDVI e CRE destacam mais nitidamente o que é e o que não é vegetação, como corpos de água, trilhas, urbanização, etc.

[IMAGEM-2]

---

3)CORREÇÕES:

A imagem B11 distingue bem entre os tipos de matas wood e forest, mas não distingue bem matas na sombra (encostas Sul) e água corrente. Os índices NDVI e CRE fazem esta distinção. É necessário então corrigir a imagem B11 acrescentando as distinções de NDVI e CRE.

OPERAÇÃO: B11+(1-NDVI)+CRE
QGIS: Raster calculator: SINTAXE: "T22JEP_20180420T132231_B11@1" + ((1-"NDVI@1")*1800) + ("CRE@1" *2600)
Agregados empiricamente os fatores ponderadores NDVI1800 e CRE2600, para que fiquem acima dos valore limites dos tipos em B11 (~2000) , e equilibrar estes índices entre si. Isto é, depende de testes de tentativa-erro-acerto até chegar aos valores de equilíbrio. Precisa ser verificado em cada trabalho com outros conjuntos de imagens diferentes, em outros locais.

RESULTADO: “B11+1-NDVI+CRE.tif” : (10m/px)

[IMAGEM-3]

---

4)TESTAR CLASSES VISUALMENTE:
-Igual à etapa (2), para B11+1-NDVI+CRE:

Depende igualmente de fazer testes de tentativa-erro-acerto até chegar aos limites significativos e adequados à realidade que possa ser observada nas imagens naturais.
QGIS: Raster Properties :
Style / Categorized Renderer: Equal Interval
Color interpolation: Discrete : testar manualmente classes e valores limites

Resultado para esta imagem do caso-teste:
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”:
QGIS: 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:
QGIS: Vector / Geoprocessing / Buffer : fixed distance
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, como destino; -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 mesma CRS no projeto)
Gera SHP com polígonos para as 2 classes somente.

---

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”

[IMAGEM 7 e abaixo]
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:
QGIS: Table Manager + Field Calculator : landuse=forest; natural=wood;
-Remover tags sem uso (area; classe).
Salvar com CRS para OSM: 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;
-verificar existência de objetos de landcover previamente mapeados no OSM (como landuse=grass;orchard;…, natural=grassland;scrub;wetland;…, etc): se houver, não substituí-los ou alterá-los; ao contrário, remover no novo conjunto os polígonos que os intersectarem.O método não objetiva alterar objetos das mesmas classes de foco (wood e forest) previamente existentes, apenas acrescentar novos que não estejam em conflito com aqueles.

VALIDAR:
Resultado: VALIDAÇÃO COMPLETA (0)

Pronto para aprovação e upload

---

5) PARA UPLOAD DOS NOVOS DADOS AO OSM:

Se aceito o método como válido para o OSM , proposto salvar em changesets exclusivos, com a nota:
“adição de polígonos de matas (wood;forest) semi-automatizados com controle de parâmetros e validação manuais”

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:

-A área da imagem cobre cerca de 100km x 100km. A imagem classificada tem 10m/px. O método permite obter cerca de 250 novos nós de polígonos por km2 (exemplo Monte Alegre dos Campos: 548km2, 136.853 nós), o que para a área inteira significa cerca de 2.500.000 nós. É algo que levaria muito tempo para ser mapeado no OSM somente manualmente sobre imagem Bing (ou muito possivelmente nunca viesse a ser mapeado). É possível ainda pensar em reduzir esta quantidade de nós, simplificando sobretudo linhas mais retas.

Na minha opinião de todo modo é muito recomendável simplificar e eliminar pequenos clusters de pixels, como <40x40m (~16px), que estejam “internos às matas”, ou inclusive descartar todos destas dimensões “fora das matas”, pois:
-o OSM usa simplificação de áreas de terreno, exemplo maior são 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 é mapear as áreas grandes de mata densa e mais homogêneas, sem invadir áreas de campos ou urbanização, o que já é suficientemente bem obtido pela limitação dos índices, não sendo necessários portanto este tipo de pequenas “ilhas”.

-Na filtragem, seria melhor 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 distinção “matas” x “o que não é mata”, é bem precisa.
A distinção matas “wood” x matas “forest”, por outro lado, é muito boa na maior parte, mas nem sempre permanece exatamente a mesma, sob os mesmos valores limites, em todas as porções da imagem. Há poucas situações em que uma é classifica no lugar da outra, por variações dos índices. Uma distinção absoluta poderia ser comprovada somente com uma análise completa de campo. O que provavelmente nunca será viável para o OSM, e pouco provável de ser feito mesmo em qualquer órgão ambiental governamental.

Em resumo, somente por mapeamento totalmente manual e sobre imagem natural como Bing, sem recorrer aos recursos de maior precisão de distinção de tipos vegetais que sensores e índices como do Sentinel proporcionam para uma classificação mais precisa, e sem contar com ferramentas de desenho automatizado (vetorização) sobre classes previamente determinadas, com controle do usuário, para desenhar grandes áreas, como o QGIS oferece:
-dificilmente se mapeará semelhantes grandes áreas de matas;
-se ainda assim forem manualmente mapeadas, dificilmente terão a precisão de distinções que sensores específicos de satélite e indices derivados oferecem como possibilidade.

Resta a comunidade avaliar:
-se é desejável e/ou aceitável esta forma de mapeamento semi-automatizado de matas;
-se é aceitável o grau de imprecisão inerente por generalização e simplificação.

Mesmo com o limite de precisão inerente ao método, pode-se observar aí melhores resultados de mapeamento, tanto em qualidade da geometria como em distinção dos tipos, do que comumente em mapeamentos manuais existentes no OSM, em geral bastante simplificados, pela própria limitação dos recursos de imagem natural.

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.