Vetorização de matas com Sentinel-2
Introdução
Proposta de método de mapeamento para vetorização semi-automática de matas (natural=wood e landuse=forest) com imagens Sentinel-2
Publicado na talk-br em:
https://lists.openstreetmap.org/pipermail/talk-br/2018-September/012423.html
Em desenvolvimento.
***IMPORTANTE***:
LER COMENTÁRIOS NA PÁGINA DE DISCUSSÃO DESTE ARTIGO
https://wiki.openstreetmap.org/wiki/Talk:Vetorização_de_matas_com_Sentinel-2
OBJETIVO:
Contribuir para o mapeamento no OSM de coberturas de matas,
gerando geometrias em modo semi-automatizado para mapeamento exclusivamente voltado a áreas de matas densas,
encontradas principalmente em interiores de território, isto é, áreas não densamente urbanizadas,
com diferenciação de vegetação para as classes natural=wood e landuse=forest,
utilizando imagens e índices apropriados fornecidos pelos Satélites Sentinel-2,
sendo que feitos manualmente pelo usuário o controle de parâmetros e a validação completa.
(Caso não desejar ler as explicações sobre a proposta, pode pular às etapas no passo-a-passo, ou ao resultado final, onde se encontram os links de acesso aos arquivos .osm resultantes)
JUSTIFICATIVA:
Áreas de matas (wood e forest) de modo geral são pouco mapeadas no OSM, pelo próprio fato de cobrir grandes áreas,
e o mapeamento em geral concentrar-se ao redor de áreas mais urbanizadas.
Permanecem assim grandes vazios, sobretudo em um país de dimensões continentais como o Brasil,
onde uma enorme maioria do território não é urbanizado, sendo dominantemente rural, ou mesmo natural.
Quando são mapeadas tais áreas, pode-se observar frequentemente as seguintes principais limitações:
- geometrias bastante simplificadas;
- não muitos casos de adequada distinção entre wood (mata natural) e forest (mata cultivada), com base unicamente em imagens em cores naturais (sobretudo em regiões onde ocorram os 2 tipos de matas em áreas densas adjacentes ou interpenetradas, necessitando maior conhecimento para fazer a distinção).
Para distinguir vegetação adequadamente, como matas do tipo natural=wood (matas naturais, mais velhas, com menor atividade biológica) e landuse=forest(matas cultivadas, mais jovens, com maior atividade biológica), existem recursos próprios de sensoriamento remoto voltados à detecção de vegetação, como os do satélite Sentinel-2.
As imagens Sentinel são atualizadas quase diariamente, mais recentes, portanto, que as encontradas no OSM nas camadas Bing e similares.
Apresentam licença compatível com OSM, conforme indicado na wiki OSM:
https://wiki.openstreetmap.org/wiki/Sentinel-2
https://wiki.openstreetmap.org/wiki/User:Ff5722/Using_Sentinel-2_imagery
Existem ainda ferramentas para processamento de imagens e vetores, como fornecidas pelo QGIS,
e validação específica para o OSM, como fornecidas pelo JOSM.
A resolução das imagens geradas no processo segue a originalmente fornecida pelo Sentinel, de 10m/px, e gera detalhes geométricos com cerca de 1nó/10m, resolução compatível com mapeamento de landcover.
O método não se destina a mapear objetos com detalhes menores que a resolução das imagens, 10m/px, os quais devem ser mapeados nas imagens em cores naturais disponíveis ao OSM.
O OSM segue a prática de simplificação para grandes áreas, como coast-lines, boundary, tipos de landcover, etc.
Neste método é proposto ainda a simplificação das geometrias usando as ferramentas do QGIS e do JOSM.
Ainda assim, o resultado possibilita melhor precisão de geometrias do que muitas vezes encontrado atualmente no OSM para este tipo de áreas,
ou como gerados com imagens Landsat.
O processo todo envolve controle manual de parâmetros e de validação. Fornece uma distinção segura entre matas e o que não é mata (como urbanização, campos, lavouras, corpos d'água), não gerando conflitos com o mapeamento de demais áreas.
O processo pode gerar cerca de 150 nós por km2 (isto em municípios rurais com razoável variação de cobertura vegetal, em geometrias e tipos, como nos testados;
em áreas rurais com vegetação mais homogênea, ou mais extensa porém uniforme, ou em municípios mais densamente urbanizados, provavelmente bem menos).
O Brasil tem 8.516.000km2 e 5570 municípios. A média de tamanho dos municípios é de 1529km2.
Isto significa que o método pode aportar cerca de 230.000 nós por município em média, para áreas de matas.
Não é de esperar encontrar 5570 mapeadores no Brasil que aportem 230.000 novos nós em cada município.
Para cobertura vegetal, ou para qualquer outra classe de elementos geográficos.
E dificilmente haverá no mundo, para cada município ou região, um mapeador que dê conta de mapear toda sua cobertura vegetal.
Assim, é de se esperar que não fossem manualmente mapeadas todas áreas de matas; exclusivavemente de modo manual possivelmente nunca o sejam.
CONSIDERAÇÕES GERAIS:
Os valores limites para distinção de classes de matas (wood e forest)
são obtidos empiricamente, por amostragem e por verificação visual, em cada conjunto de imagens (data set).
O método exige, portanto, controle ativo do usuário, e que sejam refeitos os passos para cada conjunto diferente de imagens (data set).
Não é, pois, um método totalmente automatizado. A qualidade final dos resultados depende diretamente da avaliação do usuário.
Este método de edição aqui apresentado propõe também, necessariamente, que ao final de cada processo,
seja disponibilizado à comunidade link temporário para o arquivo final "*.osm", devidamente taggeado e 100% validado, zipado, sem subir ao OSM,
publicando-o na lista oficial talk-br em tópico identificado com nome "Submissão de resultado de vetorização de matas: região de <nome>",
para que a comunidade, se assim o desejar, possa verificar diretamente os resultados na íntegra.
Após este período de verificação, não havendo objeções sobre o que é ordinariamente exigido pelo OSM,
pede-se considerar o resultado aprovado para upload dos novos dados.
O método possibilita que mapeadores interessados possam mapear várias áreas de mata de municípios ou regiões de seu interesse,
adequadamente e mais rapidamente, e com melhor qualidade do que muitas vezes encontrado em mapeamento exclusivamente manual de vegetação no OSM.
Reiterando, este método destina-se somente a áreas de matas tipo natural=wood e landuse=forest.
Para tanto, submeto esta proposta à apreciação da comunidade, esperando obter aprovação para uso deste método para mapeamento no OSM.
Etapas do processo: passo-a-passo
Testes feitos por Biomas do Brasil, conforme http://snif.florestal.gov.br/pt-br/os-biomas-e-suas-florestas.
Bioma / Local do teste:
1) Amazônia : Municípios de Porto Velho/RO e Humaitá/AM : processada vetorização; 100% validado para OSM.
2) Cerrado : Município de Ribeirão Cascalheira/MT : processada vetorização; 100% validado para OSM.
3) Mata Atlântica : Municípios de Bom Jesus e Monte Alegre dos Campos/RS : processada vetorização; 100% validado para OSM.
4) Caatinga : Município de São Raimundo Nonato/PI : Não apresentou resultado eficiente; matas dominantemente arbustivas, copas afastadas, indivíduos esparsos.
5) Pampa : Município de Santana do Livramento/RS : (a fazer vetorização e validação OSM) .
6) Pantanal : Município de Poconé/MT : (a fazer vetorização e validação OSM) .
O processos é por tentativa, erro-acerto. Pode, portanto, necessitar variar um pouco em alguns aspectos de etapas utilizadas, para cada região, bioma, ou imagem, onde será aplicado. O objetivo é chegar a uma adequada distinção de classes em valores de pixels, que permita gerar adequada vetorização de polígonos para os respectivos elementos de vegetação.
Segue abaixo descrição geral e resumida das etapas do processo, para o
caso exemplo: Bioma Mata Atlântica.
(Para descrição completa dos passos efetuados em cada bioma, acessar página de registros:
https://wiki.openstreetmap.org/wiki/User:SergioAJV/Sentinel-2_vectorizing_tests)
No QGIS
Caso exemplo utilizado nesta página descritiva geral: Bioma Mata Atlântica.
1)ESCOLHA DA IMAGEM - SENTINEL:
Fonte e descrição dos dados:
Fonte: https://earthexplorer.usgs.gov/ (necessita cadastro)
Data Set Sentinel: 14 imagens (B01 a B12; B08A; TCI):
Cloud Cover: 0-10%
Imagens e índices utilizados diretamente no processo:
Bandas: B04; B08; B11 (Identificador genérico: T2..._B11.jp2
Índices:
NDVI: Normalized Difference Vegetation Index:
Fórmula (Bandas envolvidas): NDVI = (B08-B04)/(B08+B04)
EVI2: Enhanced vegetation Index:
Fórmula (Bandas envolvidas): EVI2 = 2.5*((B08-B04)/(B08+(2.4*B04)+1))
Ver também: https://www.sentinel-hub.com/develop/documentation/eo_products/Sentinel2EOproducts
2)MARCAR ÁREAS COM NUVENS:
A avaliação com a Banda B11 fica seriamente prejudicada onde há qualquer nebulosidade, causando "buracos" no resultado.
Usar a imagem Banda B10, especialmente indicada para distinguir nuvens. Se houver nuvens, marcar polígnos em SHP próprio, para posterior mapeamento manual no JOSM.
Resultado: Sem nuvens nos municípios em foco.
3)EXAMINAR VALORES DO RASTER POR AMOSTRAGEM:
Criar um layer .shp de polígonos para marcar áreas amostrais para os objetos típicos, conforme listados na tabela abaixo,
sobre os quais serão verificados os valores de pixels. Procurar áreas de matas típicas mais nitidamente identificáveis no Bing;
áreas menos imediatamente identificáveis, como onde os tipos de matas se interpenetram; e áreas de objetos que "não são matas",
como campos, açudes, urbanização, etc.
Anotar valores limites pré-aproximados de distinções de classes dos objetos amostrais, em uma tabela como abaixo:
Caso exemplo: Mata Atlântica Municípios de Bom Jesus e Monte Alegre dos Campos/RS SENSOR ÍNDICES PONDERAÇÃO | B11 | NDVI | EVI2 | B11+((1-NDVI)*4000)) min.QGIS|512 |+0.177 | +0.281 | 1636 max.QGIS| 2709| +0.785| +1.707| 5710 Matas: wood |1200/1800| 0.30/0.80| 1.20/1.60 | 2055/2575 corte|~ > 1200 | | forest | 150/1200| 0.70/0.80| 1.30/1.78 | <=2055 M.Sombra| 150-400 | 0.30/0.65| M.Sol |1200-1800| 0.65/0.80| corte|~ < 1850 | | Não é mata (null): campos |1500/3000| 0.45/0.63| 0.80/ 1.20| 3100/4200 urbaniz.|1700/3000|-0.10/0.20| corte|~ > 600 | | açude | 50/150 |-0.10/0.10|-0.40/-0.30| 5000/5400 rio | 150/300 |-0.30/0.60|-0.10/ 0.00| 3800/4300
IMPORTANTE: Obtenção dos valores de classes a partir de observação empírica.
A cada diferente conjunto de imagens a utilizar, o processo deve ser refeito e os valores avaliados novamente.
4)CORREÇÃO DE DISTINÇÕES NA IMAGEM:
Para boa separação de classes de vegetação e demais elementos geográficos,
pode ser necessário juntar as qualidades próprias de distinções de cada imagem e/ou índice, gerando nova imagem.
Foi observado que há casos, conforme os biomas, em que basta imagem como B11;
em outros biomas, necessita usar índices como NDVI ou EVI2, ou agregar índices e imagem, para melhor distinção.
Por exemplo, no caso teste Mata Atlântica, acima:
-B11 sozinho apresentou boa distinção wood/forest, porém não diferencia bem corpos d'água,
ficando estes valores de pixels próximos da faixa de valores de matas, o que causa misturas pontuais, no processo de geração de polígonos.
Além disso, B11 tem resolução um pouco mais limitada, de 20m/px.
-Os índices NDVI/EVI agregam maior distinção de "vegetação" versus "o que não é vegetação",
como corpos d'água, e maior resolução, 10m/px;
Porém em muitos casos não destacam tanto a distinção wood/forest (florestas velhas e novas).
Assim, pode ser conveniente agregar imagem e índice.
Neste caso, foi utilizado um fator ponderador de x4000 para NDVI, testado empiricamente.
O suficiente para filtrar e separar o que não é vegetação.
Operação: Imagens B11 + ((1-NDVI)*4000)
No QGIS - Raster calculator:
("T22JEP_20180420T132231_B11@1" + ((1-"NDVI-20180420T132231@1")*4000))
Resultado: B11NDVI4000.tif c/ resolução da imagem resultante de 10m/px
Testar e anotar os valores amostrais limites nesta imagem e completar na tabela conforme etapa anteriores.
IMPORTANTE: Esta etapa é um pré-processamento, é feita com observação e testes empíricos de valores, ajustando-os e preparando-os para filtragem de classes, conforme as amostras previamente conhecidas.
Baseia-se em testagem visual, empírica, por tentativa-erro-acerto.
Neste caso, tanto para o valor para o fator ponderador de NDVI (aqui x4000), quanto para os valores limites das classes.
A cada diferente conjunto de imagens a utilizar, o processo deve ser refeito e os valores avaliados novamente.
5)TESTAR VISUALMENTE DISTINÇÃO EM 3 CLASSES:
Testar visualmente os valores limites aproximados na imagem obtidos nas etapas anteriores para destacar 3 classes:
No QGIS - Classes: Layer Properties / style / pseudocolor / discrete: 3 classes :
1-forest: testado: < = 2055
2-wood : testado: < = 2575
0-null: tudo o que não é mata.
Alterar empiricamente, por testes, para valores próximos, para mais e para menos,
de modo a chegar a uma distinção visual suficientemente representativa,
tanto aplicada às áreas amostrais como à imagem toda.
6)REDUZIR ÀS CLASSES TESTADAS:
Caso exemplo: Mata Atlântica / Municípios de Bom Jesus e Monte Alegre dos Campos/RS
Reprocessar a imagem anterior, gerando uma nova imagem com os valores limites das 3 classes encontrados na etapa acima.
Operação:
No QGIS - Raster calculator (modelo do código completo abaixo):
(2* ( "B11NDVI4000@1" <= 2055)) +(1* (( 2055 < "B11NDVI4000@1" ) AND ( "B11NDVI4000@1" <= 2575)) ) +(0* ( "B11NDVI4000@1" > 2575))
Resultado: B11-3class.tif c/3 classes:0,1,2
IMPORTANTE: Obtenção dos valores para as 2 classes a partir de observação empírica.
A cada diferente conjunto de imagens a utilizar, o processo deve ser refeito e os valores avaliados novamente.
7)CRUZAR O RASTER PROCESSADO COM OS VETORES EXISTENTES NO OSM, DE RIOS E ESTRADAS, PARA MELHORAR RECORTES NAS MATAS:
Operações:
1- OSM overpass: "(highway=* or waterway=*) and (type:way or type:relation)"
2- QGIS: Vector / Geoprocessing / Buffer : fixed distance : 0,0001 (10m)
3- QGIS: Table Manager + Field Calculator : Adicionar campo "classe", com valor ="0", para todos os polígonos.
Sobre cópia da imagem anterior processada:
4- QGIS: Raster / Conversion / Rasterize : SHP-buffer; campo "0"(zero); tiff destino
Resultado: B11-3cl+overpass.tif
8)FILTRAGENS - SIMPLIFICAR ÁREAS, REMOVER PIXELS ISOLADOS:
Operações: ~1min cada
QGIS: Processing / Toolbox :
GDAL / Analysis / Sieve
GRASS / r.neighbor
SAGA / Raster filter / Majority
Caso exemplo: Mata Atlântica / Municípios de Bom Jesus e Monte Alegre dos Campos/RS.
Sobre imagem anterior processada:
Etapas de filtragem:
1- Neighbors: Operation: minimum / Size:3 / Square 2- Sieve : Threshold: 20 / Pixel connection: 4 3- Neighbors: Operation: maximum / Size:3 / Circular 4- Sieve : Threshold: 20 / Pixel connection: 8 5- Sieve : Threshold: 8 / Pixel connection: 4
Resultado: 3 classes filtradas: B11-Sieve.tif (c/10m/px)
NOTA: Não foi possível eliminar casos de pixels pertencentes a área mesmo valor conectados por uma única conexão em quina, formando um anel.
Os fitros, ao corrigir quina em uma área, criam quinas em outras.
Isto gera erros de self-intesecting ao vetorizar para polígonos/multipolígonos.
Estes erros precisam então ser corrigidos nas etapas seguintes, de processamento do layer vetor gerado a partir do raster.
9) ELIMINAR CLASSES NÃO UTILIZADAS:
Igualar "zero" a "null"(ajuda a reduzir a quantidade de polígnos no poligonize):
OPERAÇÃO: Sobre imagem anterior processada:
-QGIS: gdal_translate : -a_nodata < x > : (x = valor da classe p/null)
Resultado: Raster Classes 1+2+null : B11-12null.tif
10)VETORIZAR RASTER:
OPERAÇÃO: Sobre imagem anterior processada:
-QGIS:: Raster / conversion / poligonize : / campo "classe" (manter mesma CRS no projeto) ~2min
Resultado: B11polig.shp
ADICIONAR TAGS OSM:
-Adicionar campos "landuse" e "natural" respectivamente para "forest" e "wood".
-Remover campos sem utilidade para OSM.
-QGIS:: Table Manager + Field Calculator : landuse=forest; natural=wood;
11)GENERALIZAR POLÍGONOS:
Despixelizar e simplificar.
-QGIS:: GRASS commands -> v.generalize - (https://grass.osgeo.org/grass64/manuals/v.generalize.html):
A) Suavizar / despixelizar polígnos:
v.generalize.smooth / Sliding Averaging : look_ahead = 3 (mín=3, sempre ímpar) : slide = 0.7 (suaviza) : ~5min
Resultado: B11slide.shp
B) Diminuir quantidade de nós :
v.generalize.simplify : Douglas Reduction : Threshold=1 Red=50% : ~5min
Cuidar se bordas longas mão ficam muito retas.
Resultado: B11simplify.shp
12)ELIMINAR ERROS GEOMÉTRICOS DE SHP VETORIZADO A PARTIR DE RASTER:
Erros tipo auto-intersecção de polígnos em quina (self-intersect; ver nota acima na etapa filtragem de raster)
Adicionar o SHP vetorizado a partir do raster em "VIRTUAL LAYER" com seleção "GEOMETRY":
Polígonos com geometrias válidas, sem overlapping / self-intersecting rings
(https://gis.stackexchange.com/questions/234757/how-to-fix-polygon-with-ring-self-intersection-in-qgis)
Operações:
A) QGIS: Layer/add layer/add virtual layer:
SYNTAXE:
SELECT landuse,natural, ST_Buffer(geometry,0.0) as geometry FROM B11simplify ;
NOTA: Nome layer sem hífens, sem extensão .shp (B11simplify).
Salvar em novo SHP, com CRS= WGS84-EPSG:4326 : B11virtgeom4326.shp
B) QGIS: Plugin: Topology-checker:
Rule: must not have invalid geometries
(Rule: must not have multi-part geometries)
Resultado: 100% válido
13)REMOVER POLÍGONOS NAS BORDAS DO LAYER:
Eliminar todos polígonos de borda, para evitar truncados.
Salvar em novo SHP
Resultado: B11-TOTAL-4326.shp : 45823 features
Opcionalmente restringir polígonos a município, com máscara:
-QGIS:Spatial query : within
Salvar cada município em novo SHP em WGS84-EPSG:4326 :
B11-BJ-4326.shp : 8483 features
B11-MA-4326.shp : 2690 features
No JOSM
1)ABRIR SHP:
-Sem baixar dados do OSM
-salvar como .osm
Já entra como polígono ou multipol(outer/inner).
2)SIMPLIFICAR:
Select: type:way / Simplify way
3) VALIDAÇÃO:
Somente dos novos dados.
Passar validador padrão completo.
4) TESTES DE CONFLITO COM O EXISTENTE NO OSM:
A executar, se aprovado o método.
-Baixar dados do OSM;
-verificar existência de objetos de landcover previamente mapeados no OSM (como landuse=* ou natural=*, etc): se houver conflito com existentes, não substituir ou alterar objetos existentes; ao contrário, remover no novo conjunto os polígonos que os intersectarem.
VALIDAR: até validação completa dos novos objetos
Resultado Final
Caso: Bioma Mata Atlântica / Município: Bom Jesus, RS
347.763 objects : Rel:671/Ways:11.153/Nodes:335.939
Área do município: 2.626 km2
Densidade de nós gerada: 128 nós/km2
Arquivo .osm: 33,3 MB
Resultado: 100% VALIDADO
Link para verificação do arquivo validado .osm: https://www.dropbox.com/s/ukrhx67x8jx1so2/B11-Bioma-MAtlantica-BJ-valid.osm.zip?dl=0
(B11-Bioma-MAtlantica-BJ-valid.osm.zip)
Para os demais biomas, seus arquivos validados e descrição completa dos passos, acessar página de registros:
https://wiki.openstreetmap.org/wiki/User:SergioAJV/Sentinel-2_vectorizing_tests