osmc:Metodologia/Algoritmo SQL/Issues
Issues e discussões conceituais, antes de se partir para a implementação.
Issue 01 - Uso ótimo do hInt64
Nesta issue buscamos garantir melhor ocupação nos 57 disponíveis de geocódigo hInt64. Os blocos hierárquicos de bits seriam os seguintes, concatenados:
$country$L0$grid
- O bloco
$gridnão pode ser otimizado. Possui como principal requsito a garantia de um número par de bits, já que inicia em 1 m2 e não pode (?) finalizar em retângulo, sua L0 precisa ser também quadrada.
- o bloco
$L0contém a cobertura de nível L0 do país. Como o máximo são 16 células L0, podemos manter com reserva seus 4 bits. Podem surgir casos onde zero bits são necessários, ou casos onde mais de 16 células são necessárias, para usar algumas células de tamanho variável.
- o bloco
$contryjá foi otimizado. Os códigos de país já foram reduzidos de 10 para 8 bits. Por exemplo Brasil era id_country=76, agora será id_country=1.
A decisão de descartar ou não 1 bit recai sobre o bloco L0. Temos 8 bits de country e o máximo de Lmax níveis hierárquicos depois de L0 para chegar na célula de 1 m2 do país, portanto, Lmax*2 bits.
Pensando em dígitos hexadecimais, são 2 dígitos para $country, Lmax/2 dígitos para $grid e o restante para totalizar 14 no L0 (4 bits do dígito hexa vezes 14 resulta em 56 bits). No caso do Brasil são mais 20 níveis depois do L0, ou seja, vai até o L20, portanto Lmax*2 = 20*2 = 40 bits. 56-40-8 = 8 bits para L0. No caso da Colômbia são 19 níveis.
Bloco L0
Se a cobertura do país for exatamente uma célula quadrada, pode ser ausente, mas em geral sua geometria precisa ser otimizada, para mais de 45% da área de L0 ser ocupada pelo polígono do território. Sem essa exigência de ocupação as aplicações de abrangência global perdem performance. As interseções entre países precisam ser minimizadas.
WITH m as (
select pais, case when count(*)>16 then 16 else count(*) end*max(area_km2) as mx
from report.v001_osmc_coverage_l0_list group by 1
)
SELECT v.pais, count(*) n, sum(v.area_km2) as km2, round(100.0*sum(v.area_km2/m.mx))::text||'%' perc_coberto
FROM report.v001_osmc_coverage_l0_list v INNER JOIN m ON m.pais=v.pais
GROUP BY 1;
pais| n | km2 | perc_coberto BR | 18 | 8729499 | 50% CM | 14 | 466041 | 48% CO | 16 | 1886770 | 44%
O máximo de células é 16, o mínimo 2 (em países com fronteiras instáveis o mais seguro é forçar reserva). Em países com ilhas e células de cobertura com baixo percentual de uso (e possibilidade de representação com uma só célula-filha) pode-se adotar mais que 16 células de cobertura, emulando um nível L0 com mosaico de pedaços, como no caso do Brail. Vamos discutir os três casos a seguir.
A questão principal é se vai permitir ou não quantidade ímpar de células.
No caso do Brasil a codificação L0 requer 2 conjuntos de ilhas, depois das 16 coberturas territoriais. Basta 1 bit para codificar se é uma "falsa célula" ou não. O extremo sul do Brasil, que recebeu o rótulo f, teve sua área subdivida em mosaico para comportar as "falsas células" 10=g e 11=h. Na prática são todas f, e o geocódigo só "assume sua representação concreta" a partir do segundo dígito: fT (e seu litoral fY) para o território extremo sul, fN para as ilhas oceânicas sob jurisdição do estado de Espirito Santo, e fP para o arquipelogo de São Pedro e São Paulo.
pais | b16h | area_km2 ------+------+---------- BR | 00 | 275241 BR | 01 | 757637 BR | 02 | 645014 BR | 03 | 132780 BR | 04 | 619969 BR | 05 | 1075008 BR | 06 | 1099512 BR | 07 | 1054743 BR | 08 | 103225 |
pais | b16h | area_km2 ------+------+---------- BR | 09 | 333725 BR | 0a | 1099383 BR | 0b | 723005 BR | 0c | 19403 BR | 0d | 714206 BR | 0e | 54071 BR | 0f | 17440 BR | 10 | 1519 BR | 11 | 3618 |
Issue 02 - Função estimadora de cobertura interior e quadrado envolvente
Para fazer estimativas sobre células de km2. A função f(x) que retorna potências side e next:
- side: estimativa para conferir se a área de uma célula recortada (ou seja célula de cobertura) é ou não uma célula interior (sem recortes), apenas pelo valor de sua área.
- next: para chutar o quadrado envolvente inicial de um país (pode dispensar cobertura), ou de um município.
Abaixo a função com exemplo da área do Brasil.
-- A = 8510417.771 km2:
with t as (select sqrt(8510417.771 *1000^2) as x)
select p as side_pow_exp, 2^p as side_pow_val,
round( (100.0*(2^p-x)/2^p)::numeric , 2)::real side_diff_perc,
p+1 as next_pow_exp, 2^(p+1) as next_pow_val,
round( (100.0*(2^(p+1)-x)/2^(p+1))::numeric , 2)::real next_diff_perc
from (select round(log(x)/log(2)) p,x from t) t2;
-- side_pow_exp | side_pow_val | side_diff_perc | next_pow_exp | next_pow_val | next_diff_perc
-- -------------+--------------+----------------+--------------+--------------+----------------
-- 21 | 2097152 | -39.11 | 22 | 4194304 | 30.45
-- A = 757637 km2:
side_pow_exp | side_pow_val | side_diff_perc |
--------------+--------------+----------------|
20 | 1048576 | 16.99 |
-- A = 1099512 km2:
20 | 1048576 | 0 |
No caso da área do Brasil (de 8510417.771 km2) temos o next como estimativa de quadrado envolvente, com lados de 2^22 km:
- quadrado com 2^22 metros de lado contém o Brasil.
- Primeiro dígito base4 é a cobertura de 4 quadrados com 2^21 metros cada
- Primeiro dígito base16 é a cobertura de 16 quadrados com 2^20 metros de lado cada.
Nas áreas de células conferimos qual tem side_diff_perc=0, que foi apenas o caso de A = 1099512 km2.
Exemplos, onde seriam esperados expoentes pares, da ordem de 20 para cobertura do Brasil e 18 para Camarões.
| pais | is_contained | cbits | b16 | area_km2 | exp e diff% |
|---|---|---|---|---|---|
| BR = 76 = 0001001100 | f | 0001001100.00000110 | 06 | 1099512 | 20 com 0% |
| BR = 76 = 0001001100 | t | 0001001100.00000000 | 00 | 275241 | 19 com -0.07% |
| BR = 76 = 0001001100 | t | 0001001100.00000001 | 01 | 757637 | 20 com 16.99% |
| ... | |||||
| CM = 120 = 0001111000 | f | 0001111000.1001 | 9 | 68719 | 18 com 0% |
| CM = 120 = 0001111000 | t | 0001111000.0001 | 1 | 2469 | 16 com 24.18% |
| ... |
A interseção da geometria da célula com o polígono do país tem a sua área indicada pela coluna area_km2. A coluna is_contained é um flag, verdadeiro quando a célula tem interseção com as bordas do país. Analisando as áreas do Brasil, o flag foi falso apenas para os de área, cuja raíz quadada é o tamanho do lado do quadrado, e uma das potências de 2 utilizadas, . Analogamente em Camarões (CM), .
Issue 03 - Redução do excesso de pontos por efeito-varal
Conforme "How to reduce the clothesline-effect in ST_Transform?", pode-se minimizar o efeito-varal. A ideia aqui nesta issue é usasar simplificação de curva, para controle e eventual complementação do algoritmo.
select hlevel, count(*) n_cells, round(sqrt(max(st_area(geom)))/1000.0) side_km,
round(avg(ST_NPoints(geom))) geom_npts,
round(avg(ST_CharactDiam(geom))) as geom_chdiam, -- relative ref.
round(avg(ST_NPoints(geom4326))) g4326_npts,
round(avg(ST_NPoints(st_simplify(geom4326,0.00000005)))) g4326simpl_npts
-- 0.00000005 is SRID-metric absolute ref
from grid_cm.all_levels
group by 1 order by 1;
hlevel | n_cells | side_km | npts_geom | npts_geom | npts_simpl
--------+---------+---------+-----------+-----------+------------
0 | 14 | 262 | 5 | 749 | 293
0.5 | 28 | 185 | 5 | 769 | 280
1 | 56 | 131 | 5 | 749 | 152
1.5 | 112 | 93 | 5 | 769 | 141
2 | 224 | 66 | 5 | 749 | 77
2.5 | 448 | 46 | 5 | 769 | 71
3 | 896 | 33 | 5 | 749 | 39
3.5 | 1792 | 23 | 5 | 769 | 36
4 | 3584 | 16 | 5 | 749 | 20
4.5 | 7168 | 12 | 5 | 769 | 19
5 | 14336 | 8 | 5 | 749 | 11
5.5 | 28672 | 6 | 5 | 769 | 11
Conclusão: refazer a nossa função de ST_Transform_resilient() para os casos de cache, talvez não de alta performance
- no final incluir ST_SimplifyPreserveTopology, portanto mais esse parametro simplify_tolerance precisa ser incluso. Se o uso se restringe a células, pode ser ST_Simplify, mais rápida. Decidir se existe como criar uma tolerância relativa ao perímetro, abarcando a intenção da densidade. Provavelmente um fator da densidade*perimetro resolva. O foco são as células da grade.
- com base no parâmetro métrico simplify_tolerance decidir se não faz nada, visto que depois do simplify voltaria ao número de pontos original. É de fato o esperado para a grade de hlevel de 6 em diante no exemplo. Foco nas células.
Na resposta da questão quebrar em duas, "general" depois "grids" (exemplificar com níveis inteiros pulando os meios). Nas grades supor que só um eixo precisa de pontos ... Ver eles no QGIS ou pelo menos nos dados de 11 pontos. O segmentize aplicado a uma só direção é complexo, não vale a pena.
DROP FUNCTION IF EXISTS st_transform_resilient
;
CREATE FUNCTION public.st_transform_resilient(
g geometry,
srid integer,
size_fraction float DEFAULT 0.05,
tolerance float DEFAULT 0 -- 0.00000005
) RETURNS public.geometry
LANGUAGE sql IMMUTABLE
AS $f$
SELECT CASE
WHEN COALESCE(tolerance,0)>0 THEN -- ST_Simplify(celula) ou
ST_SimplifyPreserveTopology(geom,tolerance) -- generico
ELSE geom
END
FROM (
-- discuss ideal at https://gis.stackexchange.com/q/444441/7505
SELECT CASE
WHEN size>0.0 THEN ST_Transform( ST_Segmentize(g,size) , srid )
ELSE ST_Transform(g,srid)
END geom
FROM (
SELECT CASE
WHEN size_fraction IS NULL THEN 0.0
WHEN size_fraction<0 THEN -size_fraction
ELSE ST_CharactDiam(g) * size_fraction
END
) t1(size)
) t2
$f$;