osmc:Metodologia/Algoritmo SQL/Issues: mudanças entre as edições
mSem resumo de edição |
m (→Issue 05 - Cálculo de dígitos base32: typo) |
||
(3 revisões intermediárias pelo mesmo usuário não estão sendo mostradas) | |||
Linha 201: | Linha 201: | ||
Migrar função para o git oficial! | Migrar função para o git oficial! | ||
== Issue 04 - Lib de grades == | |||
Desenvolvendo em [[osmc:Metodologia/Algoritmo SQL/Issue04]]. | |||
== Issue 05 - Cálculo de dígitos base32 == | |||
Por imposição do [[Discrete National Grid Systems/pt|padrão DNGS]] temos <math>S_{L}=2^{Lmax-L}</math>. No caso do Brasil ''Lmax''=20, no caso de Camarões ''Lmax=18''. Já havíamos notado isso ao apresentar a [[osmc:Metodologia/Algoritmo_SQL/Lib#Core|Lib Core]]. | |||
No caso de geocódigo logístico, que requer cobertura municipal, supondo que o nível da grade adotada para cobertura seja ''Lcover'', podemos obter a aproximação ao metro por um múltiplo de 2.5, que é a quantidade de níveis que se sobe a cada 5 bits do dígito da base32. Partindo de <math>1=2^{Lmax - Lcover - N \cdot 2.5}</math> obtemos por <math>log_2(1)=0</math> a fórmula <math>N= (Lmax-Lcover)/2.5</math>, onde ''N'' é o número de dígitos depois do primeiro para se chegar ao metro. Num caso típico a cobertura municipal no Brasil (''Lmax''=20) é de nível ''L7.5'', resultando em ''N''=5, ou seja, com 6 dígitos teremos o metro. Para uma cobertura de resolução menor, ''L7'', o valor já não é exato, ''N''=5.2. O número de dígitos para chegar no metro portanto será <math>Ndig1m = 1+\lceil (Lmax-Lcover)/2.5 \rceil</math>. | |||
<syntaxhighlight lang="sql" style="font-size: 80%;"> | |||
CREATE FUNCTION grid_br.logistic_Ndig1m(intlevel_cover int) RETURNS real[] AS $f$ | |||
SELECT array[ Lcover, ceil(Ndig1m), ceil(Ndig1m)-Ndig1m, | |||
2^(Lmax-Lcover-ceil(Ndig1m)*2.5), | |||
CASE WHEN ceil(Ndig1m)-Ndig1m=0 THEN NULL ELSE 2^(Lmax-Lcover-floor(Ndig1m)*2.5) END ] | |||
FROM ( | |||
SELECT *, (Lmax-Lcover)/2.5 as Ndig1m | |||
FROM (SELECT 20 as Lmax, intlevel_cover/10.0 as Lcover) t1 | |||
) t2 | |||
$f$ language SQL IMMUTABLE; | |||
-- gerando a tabela do Brasil: | |||
select x[1] as level_cover, x[2] as Ndigits_for_1m, round(x[4],2) as side_size, | |||
round(x[5],2) as size_before | |||
from ( | |||
select grid_br.logistic_Ndig1m(intLevel) as x | |||
from generate_series(40,90,5) t1(intLevel) | |||
) t2; | |||
</syntaxhighlight> | |||
A escolha de nivel ''Lcover'' de cada município no caso do Brasil vai ter o seguinte perfil dado pela função <code>logistic_Ndig1m()</code>. Como o metro é a maior resolução significativa, é interesante não permitir resoluções finais inferiores a meio metro: por exemplo quando a cobertura ótima ''Lcover'' for 4 ou 4.5, usar no máximo 6 dígitos ao invés de 7, já que 2 metros e 1.4 metros são resoluções razoáveis para endereços. Idem ''Lcover''6.5 e 7 etc. | |||
{| class="wikitable" | |||
|+ Número de dígitos base32 para se chegar no metro, conforme ''Lcover'' do município | |||
|- | |||
|''Lcover'' || Ndigits for 1m || ''side_size'' || ''size_before'' | |||
|- | |||
| 4 || 7 || 0.35 || 2.00 | |||
|- | |||
| 4.5 || 7 || 0.25 || 1.41 | |||
|- | |||
| 5 || 6 || 1.00 || | |||
|- | |||
| 5.5 || 6 || 0.71 || 4.00 | |||
|- | |||
| 6 || 6 || 0.50 || 2.83 | |||
|- | |||
| 6.5 || 6 || 0.35 || 2.00 | |||
|- | |||
| 7 || 6 || 0.25 || 1.41 | |||
|- | |||
| 7.5 || 5 || 1.00 || | |||
|- | |||
| 8 || 5 || 0.71 || 4.00 | |||
|- | |||
| 8.5 || 5 || 0.50 || 2.83 | |||
|- | |||
| 9 || 5 || 0.35 || 2.00 | |||
|} |
Edição atual tal como às 23h47min de 25 de junho de 2024
Issues e discussões conceituais, antes de se partir para a implementação. Ver também Discussão:DNGS/Decisões_soberanas#Issues, relativas a decisões humanas, sem algoritmo previsto.
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
$grid
nã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
$L0
conté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
$contry
já 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 ST_Transform_Resilient(
g geometry, -- the input geometry
srid integer, -- target SRID, to transform g
size_fraction float DEFAULT 0.05, -- 1/density. Density of points per charactDiam (or negative for absolute fraction).
tolerance float DEFAULT 0 -- E.g. on srid=4326 use 0.00000005.
-- ZERO=0.000000001. Smallest=0.000000005. Start cell with 0.00000002.
) RETURNS geometry
language SQL IMMUTABLE
AS $f$
SELECT CASE
WHEN COALESCE(size_fraction,0.0)>0.0 AND COALESCE(tolerance,0)>0 THEN
ST_SimplifyPreserveTopology(geom,tolerance) -- ST_Simplify enough for grid cells
ELSE geom
END
FROM (
SELECT CASE
WHEN size>0.0 THEN ST_Transform( ST_Segmentize(g,size) , srid )
ELSE ST_Transform(g,srid)
END geom, size
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$;
COMMENT ON FUNCTION ST_Transform_Resilient IS
'Avoid clothesline-effect. See problem/solution discussed at https://gis.stackexchange.com/q/444441/7505'
;
Migrar função para o git oficial!
Issue 04 - Lib de grades
Desenvolvendo em osmc:Metodologia/Algoritmo SQL/Issue04.
Issue 05 - Cálculo de dígitos base32
Por imposição do padrão DNGS temos . No caso do Brasil Lmax=20, no caso de Camarões Lmax=18. Já havíamos notado isso ao apresentar a Lib Core.
No caso de geocódigo logístico, que requer cobertura municipal, supondo que o nível da grade adotada para cobertura seja Lcover, podemos obter a aproximação ao metro por um múltiplo de 2.5, que é a quantidade de níveis que se sobe a cada 5 bits do dígito da base32. Partindo de obtemos por a fórmula , onde N é o número de dígitos depois do primeiro para se chegar ao metro. Num caso típico a cobertura municipal no Brasil (Lmax=20) é de nível L7.5, resultando em N=5, ou seja, com 6 dígitos teremos o metro. Para uma cobertura de resolução menor, L7, o valor já não é exato, N=5.2. O número de dígitos para chegar no metro portanto será .
CREATE FUNCTION grid_br.logistic_Ndig1m(intlevel_cover int) RETURNS real[] AS $f$
SELECT array[ Lcover, ceil(Ndig1m), ceil(Ndig1m)-Ndig1m,
2^(Lmax-Lcover-ceil(Ndig1m)*2.5),
CASE WHEN ceil(Ndig1m)-Ndig1m=0 THEN NULL ELSE 2^(Lmax-Lcover-floor(Ndig1m)*2.5) END ]
FROM (
SELECT *, (Lmax-Lcover)/2.5 as Ndig1m
FROM (SELECT 20 as Lmax, intlevel_cover/10.0 as Lcover) t1
) t2
$f$ language SQL IMMUTABLE;
-- gerando a tabela do Brasil:
select x[1] as level_cover, x[2] as Ndigits_for_1m, round(x[4],2) as side_size,
round(x[5],2) as size_before
from (
select grid_br.logistic_Ndig1m(intLevel) as x
from generate_series(40,90,5) t1(intLevel)
) t2;
A escolha de nivel Lcover de cada município no caso do Brasil vai ter o seguinte perfil dado pela função logistic_Ndig1m()
. Como o metro é a maior resolução significativa, é interesante não permitir resoluções finais inferiores a meio metro: por exemplo quando a cobertura ótima Lcover for 4 ou 4.5, usar no máximo 6 dígitos ao invés de 7, já que 2 metros e 1.4 metros são resoluções razoáveis para endereços. Idem Lcover6.5 e 7 etc.
Lcover | Ndigits for 1m | side_size | size_before |
4 | 7 | 0.35 | 2.00 |
4.5 | 7 | 0.25 | 1.41 |
5 | 6 | 1.00 | |
5.5 | 6 | 0.71 | 4.00 |
6 | 6 | 0.50 | 2.83 |
6.5 | 6 | 0.35 | 2.00 |
7 | 6 | 0.25 | 1.41 |
7.5 | 5 | 1.00 | |
8 | 5 | 0.71 | 4.00 |
8.5 | 5 | 0.50 | 2.83 |
9 | 5 | 0.35 | 2.00 |