SQL, extract municipio

Todos os municípios podem ser exportados a partir da carga do "planet", usando o aplicativo de terminal osm2pgsql, conforme descrito em Do zero. A seguir exemplos de localização e verificação, garantindo a qualidade da extração.

Exemplo de RR-BoaVista no stable 2018 JSON

Pelo datasets.ok.org.br/city-codes temos confirmação de que o código IBGE é 1400100 e a entrada Wikidata é Q181056, onde fica destacado que o OSM-relation-ID é 326286. Ainda assim, convém no SQL conferir se não existem outros polígonos de boundary:administrative com mesmo código IBGE. Se fossemos recuperar direto pelo OSM-relation-ID da Wikidata, convém lembrar que é negativo na representação osm_id. Outro problema é que queremos o "polígono puro", mas a relation inclui o ponto do admin_centre (falta conferir se o osm2pgsql separou as duas coisas (geometria pura de ST_Polygon).

 SELECT osm_id, ST_GeometryType(way) as geom_type, tags 
 FROM planet_osm_polygon 
 WHERE tags->>'boundary'='administrative' and tags->>'IBGE:GEOCODIGO'='1400100';

A listagem confirma apenas um item e com o perfil nas tags e na identificação (osm_id = -326286) conforme esperado.

Bom lembrar também que no Brasil todos os seus ~5,5 mil municípios são admin_level:8. Assim um jeito simples de conferir se estão todos lá e sem duplicação é:

CREATE VIEW vw_planet_osm_polygon_count_city  AS
 SELECT tags->>'IBGE:GEOCODIGO' as ibge_code, ST_GeometryType(way) as geom_type, count(*) n
 FROM planet_osm_polygon 
 WHERE tags->>'boundary'='administrative' AND tags->>'admin_level'='8' AND tags?'IBGE:GEOCODIGO'
 GROUP BY 1,2
 ORDER BY 1,2
;
SELECT count(*) FROM vw_planet_osm_polygon_count_city; --  5570
SELECT *  FROM vw_planet_osm_polygon_count_city WHERE n>1; --  0 rows

Ok está consistemte.

Exemplo de RR-BoaVista no stable 2020 hstore

Para o caso de não se ter adotado a conversão (recomendada) de HStore para JSONb, e aproveitando dados de 2020 para comparar com 2018.

 SELECT osm_id, ST_GeometryType(way) as geom_type, tags 
 FROM planet_osm_polygon 
 WHERE tags->'boundary'='administrative' and tags->'IBGE:GEOCODIGO'='1400100';
-- ok confirmado osm_id -326286

CREATE VIEW vw_planet_osm_polygon_count_city  AS
 SELECT tags->'IBGE:GEOCODIGO' as ibge_code, ST_GeometryType(way) as geom_type, count(*) n
 FROM planet_osm_polygon 
 WHERE tags->'boundary'='administrative' AND tags->'admin_level'='8' AND tags?'IBGE:GEOCODIGO'
 GROUP BY 1,2
 ORDER BY 1,2
;
SELECT count(*) FROM vw_planet_osm_polygon_count_city; --  mudou de 5570 para 5568
SELECT *  FROM vw_planet_osm_polygon_count_city WHERE n>1; --  0 rows

A mudança se refere ao balanço final de exclusões/inclusões oficiais de municípios no Brasil (? conferir!).

Polígono as-is e simplificado GeoJSON do município

O mais importante, por fim, é publicar o polígono do município no git conforme formato mais recomendado. Já temos a tabela e o osm_id, portanto basta exportar o resultado de ST_AsGeoJSON() do polígono.

 SELECT ST_AsGeoJSON(way,7) as municipio_326286_geojson
 FROM planet_osm_polygon 
 WHERE osm_id=-326286;

O limite no número de casas decimais se justifica por seu significado, ou seja, não precisamos mais que um metro de precisam nos limites de município, já que nem os limites dados por hidrografia ou por definição oficial (IBGE) chega a ter essa precisão. Pode-se testar a sensibilidade da geometria por:

 SELECT ST_AsGeoJSON(way,7)=ST_AsGeoJSON(way,9) AS compare FROM planet_osm_polygon WHERE osm_id=-326286;

Por fim, a simplificação dos pontos preservando a geometria original confere também estabilidade adicional. O melhor algoritmo para isso é ST_SimplifyPreserveTopology(way,tol). A verificação de sensibilidade ao parâmetro tol, entretanto, deve ser realizada com ST_Equals(g_original,g_simplificada), não com =, conforme esta aqui... Portanto devemos exportar como geometria stable:

 SELECT ST_AsGeoJSON( ST_SimplifyPreserveTopology(way,0), 7 ) geom 
 FROM planet_osm_polygon WHERE osm_id=-326286;

Verificando a sensibilidade em todas as geometrias de municípios do Brasil, nos parâmetros $digs e $tol:

 SELECT  ST_AsGeoJSON(way,$digs)=ST_AsGeoJSON(way,9) as geojs_cmp, 
         ST_Equals( way, ST_SimplifyPreserveTopology(way,$tol) ) as simp_cmp, 
         count(*) n
 FROM planet_osm_polygon 
 WHERE tags->'boundary'='administrative' AND tags->'admin_level'='8' AND tags?'IBGE:GEOCODIGO'
 GROUP BY 1,2
 ORDER BY 1,2

• Zero ocorrências para ambos: $digs=7 e $tol=0
• Limite, onde começam a haver ocorrências: $digs=6 (mais de 90%) e $tol=0.00000000001 (menos de 50%).

Portanto os parâmetros escolhidos (7 e 0) são adequados para o presente e o futuro. Outra opção segura seria adotar $tol=0.000000000005 para reduzir sensibilidade a "falsas modificações" na geometria.