osmc:Metodologia/Algoritmo SQL: mudanças entre as edições
(Configuração do país) |
m (→Tratamento das configurações: fT fP fN) |
||
(29 revisões intermediárias por 2 usuários não estão sendo mostradas) | |||
Linha 1: | Linha 1: | ||
Detalhamento metodológico, descrevendo os algoritmos implementados em linguagem SQL. | Detalhamento metodológico, descrevendo os algoritmos implementados em linguagem SQL. | ||
Para obtenção do código-fonte: | |||
# Versão desejada se homologada: git. | |||
# Versão corrente do servidor: se em dúvida, usar <code>pg_dump --schema-only -f dl05_dump_test.sql dl05</code>. | |||
== Configuração do país == | == Configuração do país == | ||
Linha 7: | Linha 10: | ||
* Configurações da Colômbia https://git.osm.codes/CO_new/blob/main/conf.yaml#L34 | * Configurações da Colômbia https://git.osm.codes/CO_new/blob/main/conf.yaml#L34 | ||
Os parâmetros principais do YAML, | Os parâmetros principais do YAML da Colômbia, por exemplo, são: | ||
* | <pre> | ||
* | dngsCode_definition: | ||
srid: 102022 | |||
grid_xy_origin: 3678500 970000 | |||
grid_cell_side: 524288 | |||
grid_l0_cell_ij: 02 03 10 11 12 13 20 21 22 23 30 31 32 40 41 42 | |||
grid_l0_cell_idxb16: 8 a 1 3 9 b 4 6 c e 5 7 d 0 2 f | |||
</pre> | |||
onde: | |||
* [https://postgis.net/docs/using_postgis_dbmanagement.html#spatial_ref_sys ''srid''] é a projeção igual-area do país; | |||
* ''grid_xy_origin'': é a origem no plano projetado; | |||
* ''grid_cell_side'': é o tamanho (potência de 2) da célula ''L0'' adotada; | |||
* ''grid_l0_cell'': a lista de coordenadas ''ij'' da matriz de células, indicando aquelas selecionados para a cobertura; | |||
* ''grid_l0_cell_idxb16'': a sua correspondente indexação, já expressa em base16. | |||
A parametrização YAML poderia ser lida diretamente do arquivo (pendente implementar essa estratégia); todavia, atualmente, seus parâmetros estão embutidos na função wrapper <code>[https://github.com/osm-codes/GGeohash/blob/main/src/step04def-ini.sql#L67-L82 osmc.L0cover_upsert]</code>. | |||
A configuração é efetivada pela chamada com o país específico, por exemplo Brasil <code>SELECT osmc.L0cover_upsert('BR');</code> | A configuração é efetivada pela chamada com o país específico, por exemplo Brasil <code>SELECT osmc.L0cover_upsert('BR');</code> | ||
A função de que | A função <code>osmc.L0_upsert()</code> atualiza a tabela ''osmc.coverage'' com linhas ''is_country'' e sem popular as colunas ''cindex'', ''is_overlay'' ou ''kx_prefix'' (utilizadas apenas por coberturas municipais). A tabela abaixo, gerada pelo relatório ''osmc_report.v001_osmc_coverage_l0_list'' apresentado em seguida, ilustra os casos típicos do Brasil, Camarões e Colômbia: | ||
{| class="wikitable" | |||
! pais !! is_contained !! cbits !! b16 !! area_km2 | |||
|- | |||
| BR = 1 = 00000001 || f || 00000001.00000000 || 00 || 275241 | |||
|- | |||
| BR = 1 = 00000001 || f || 00000001.00000001 || 01 || 757637 | |||
|- | |||
| BR = 1 = 00000001 || f || 00000001.00000010 || 02 || 645014 | |||
|- | |||
| ... | |||
|- | |||
| BR = 1 = 00000001 || f || 00000001.00010000 || 10 || 1519 | |||
|- | |||
| BR = 1 = 00000001 || f || 00000001.00010001 || 11 || 3618 | |||
|- | |||
| BR = 1 = 00000001 || t || 00000001.00000110 || 06 || 1099512 | |||
|- | |||
| CM = 3 = 00000011 || f || 00000011.0001 || 1 || 16259 | |||
|- | |||
| CM = 3 = 00000011 || f || 00000011.0010 || 2 || 1261 | |||
|- | |||
| CM = 3 = 00000011 || f || 00000011.0011 || 3 || 43764 | |||
|- | |||
| ... | |||
|- | |||
| CM = 3 = 00000011 || f || 00000011.1110 || e || 49039 | |||
|- | |||
| CM = 3 = 00000011 || f || 00000011.1111 || f || 14919 | |||
|- | |||
| CM = 3 = 00000011 || t || 00000011.1010 || a || 68719 | |||
|- | |||
| CO = 2 = 00000010 || f || 00000010.0000 || 0 || 8083 | |||
|- | |||
| CO = 2 = 00000010 || f || 00000010.0001 || 1 || 230987 | |||
|- | |||
| CO = 2 = 00000010 || f || 00000010.0010 || 2 || 57590 | |||
|- | |||
| ... | |||
|- | |||
| CO = 2 = 00000010 || f || 00000010.1101 || d || 196017 | |||
|- | |||
| CO = 2 = 00000010 || f || 00000010.1110 || e || 84471 | |||
|- | |||
| CO = 2 = 00000010 || f || 00000010.1111 || f || 61264 | |||
|} | |||
A coluna ''cbits'' é composta de duas partes, o código do país e o código da sua célula ''L0''. O código é traduzido em hexadecimal (base16) na coluna ''b16'', e suas células formam a cobertura do país, ou seja, a grade de nível zero. Nos [[osmc:Convenções/Grade científica multifinalitária|geocódigos de notação científica]] apare cerão sempre como primeiro dígito. | |||
Os códigos dos países são números inteiros de 8 bits (atualmente, sem compromisso com códigos ISO de 10 bits), garantindo melhor ocupação nos 57 bits (56 bits para acomodar 14 dígitos hexadecimais) disponíveis de geocódigo [[Código_Natural/Representação_interna#hInt|hInt64]] (essa representação alternativa é mais eficiente e versátil do que o tipo de dado ''varbit'' (bit string) do PostgreSQL. | |||
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 está contida no polígono do país. Ver [[osmc:Metodologia/Algoritmo_SQL/Issues#Issue_02_-_Função_estimadora_de_cobertura_interior_e_quadrado_envolvente|Issue 02 para detalhes sobre a dedução do caso cheio pela área]]. | |||
<syntaxhighlight lang="sql"> | |||
CREATE VIEW report.v001_osmc_coverage_l0_list AS | |||
SELECT isolabel_ext ||' = '||(b'0000000000000000000000'||cb10)::bit(32)::int ||' = '||cb10 as pais, | |||
is_contained, | |||
cb10::text||'.'||cb_suffix as cbits, natcod.vbit_to_baseh(cb_suffix,16) as b16, | |||
round(st_area(geom)/1000.0^2) as area_km2 | |||
FROM ( | |||
select *, substring(cbits,1,10) as cb10, substring(cbits,11) as cb_suffix | |||
from osmc.coverage | |||
) t | |||
WHERE is_country ORDER BY 1,2,cbits | |||
; | |||
</syntaxhighlight> | |||
=== Seletor de projeção === | |||
[[Arquivo:DNGS-BBOXes-exemplo2.png|thumb|380px|BBOX WGS84 sobre o globo, quando apenas Brasil e Colômbia adotavam o padrão DNGS. A BBOX original de cada país é decomposta em BBOXes "puras" e "de fronteira". As puras permitem rápida decisão (alta performance), enquanto as de fronteira, ''border boxes'' em rosa, requerem avaliação de pertinência ao polígono de fronteira.]] | |||
O padrão DNGS, conforme seu [[Discrete_National_Grid_Systems/pt#Seletor_de_jurisdição|algoritmo seletor de jurisdição]], requer que qualquer ponto sobre o globo, expresso em WGS84 sem projeção, tenha seu geocódigo resolvido conforme a jurisdição da sua posição. O algoritmo decide qual projeção usar. A verificação em grandes BBOXes, sem bordas, tem alta performance, enquanto que a decisão sobre polígonos de borda tem menor performance. | |||
A grade ''L0'' serve de referência para a criação de recortes BBOX WGS84 sobre o globo, onde as escolhas de projeção são reproduzidas de forma mais grosseira que o mosaico ''L0'', e, nas bordas, recortadas. Supondo ''pt'' a geometria de um ponto: | |||
: <code>SELECT srid FROM osmc.coverage_srid <br/>WHERE pt && geom<br/> AND (is_not_border OR ST_Intersects(pt,geom))</code> | |||
: indexação prévia com <br/><code>CREATE INDEX osmc_coverage_srid_idx1<br/>ON osmc.coverage_srid(is_not_border,geom) <br/>USING GIST (geom)</code> | |||
: a performance pode ser avaliada por <br/><code>VACUUM ANALYZE osmc_coverage_srid_idx1</code> | |||
O conhecimento prévio do SRID (pode ser por contexto) permite por fim fazer a busca na tabela <code>osmc.coverage</code>. | |||
=== Seletor de jurisdição e L0 === | |||
Uma vez selecionada a projeção, o ponto é transformado conforme SRID e confrontado apenas contra as BBOXes de células L0 daquela projeção. No caso especial de ser uma projeção adotada por mais de um país (por exemplo Itália e Vaticano), e o ponto cair em uma interseção de BBOXes, a seleção do recorte de ''L0'' correto dependerá do ponto estar contido no polígono do recorte. | |||
Justamente para contemplar os casos especiais a seleção de jurisdição acaba usando uma condição e indexação mais complexa, similar à seletora de projeção: | |||
: <code>SELECT cbits as cbits_l0 FROM osmc.coverage <br/>WHERE srid=$srid AND pt && geom<br/> AND (is_not_border OR ST_Intersects(pt,geom))</code> | |||
: indexação prévia com <br/><code>CREATE INDEX osmc_coverage_idx1<br/>ON osmc.coverage(srid,is_not_border,geom) <br/>USING GIST (geom)</code> | |||
O valor ''cbits_l0'' é o identificador de célula ''L0'' do país correto, e o identificador de país é seu prefixo. A partir do ponto ''pt'' (com geometria já no SRID do país) e ''cbits_l0'' obtemos por algoritmo GGeohash o valor ''cbits_pt'' da célula pontual, <code>cbits_pt=encode(pt,cbits_l0)</code>. | |||
'''NOTA'''. Em aplicações que não partem de um cenário global, ou de pontos que não precisam ser selecionados como interiores ou não ao território nacional, pode-se adotar um seletor mais simples: <code>grid_br.xyS_collapseTo_ijS(pt)</code> por default retorna as coordenadas IJ de L0. | |||
: <code>SELECT cbits as cbits_l0 FROM osmc.coverage <br/>WHERE srid=$srid AND pt && geom</code> | |||
: indexação prévia com <br/><code>CREATE INDEX osmc_coverage_idx1<br/>ON osmc.coverage(srid,geom) <br/>USING GIST (geom)</code> | |||
Por fim, como o país não pode ter mais do que 16 células L0, uma cadeia de condições IF ou CASE é suficiente e provavelmente mais rápido do que consultar o banco de dados. Usando <code>ST_X(pt)</code> e <code>ST_Y(pt)</code> e os valores de BBOX, podemos otimizar o processo de escolha de ''L0'' e ''encode'' do ponto. | |||
=== Tratamento das configurações === | |||
Usando o exemplo do Brasil, "grid_lc" indica que a cobertura é uma matriz de 5x5 (linhas por colunas), e "grid_l0_cell" quais são as células IJ selecionadas pelo índice "grid_l0_cell_idx": | |||
<pre> | |||
grid_lc: 5 5 | |||
grid_l0_cell: 40 41 42 43 30 31 32 33 34 21 22 23 11 12 13 02 44 24 | |||
grid_l0_cell_idx: 0 1 2 3 4 5 6 7 8 9 a b c d e fT fP fN | |||
</pre> | |||
As configurações podem ser reescritas com auxilio de uma query: | |||
<syntaxhighlight lang="sql" style="font-size: 80%;"> | |||
WITH m AS ( -- background matrix | |||
select i,j from generate_series(0,4) t1(i), generate_series(0,4) t2(j) | |||
), c AS ( -- cover cells | |||
select i,j, (natcod.baseh_to_vbit('{0,1,2,3,4,5,6,7,8,9,a,b,c,d,e,fT,fP,fN}'::text[],16))[(p+2)/3] as cbits | |||
from ( | |||
select i,j, strpos('40 41 42 43 30 31 32 33 34 21 22 23 11 12 13 02 44 24', i::text||j::text)::float as p | |||
from m | |||
) tt where p>0 | |||
) SELECT array_agg(a order by i) -- IJ_to_cover matrix | |||
FROM ( | |||
SELECT m.i, array_agg( c.cbits order by m.j) a | |||
FROM c RIGHT JOIN m ON c.i=m.i AND c.j=m.j | |||
GROUP BY 1 ORDER BY 1 | |||
) t4 | |||
; | |||
</syntaxhighlight> | |||
A matriz resultante <code>IJ_to_cover={{NULL,NULL,1111101,NULL,NULL}, ..., {0000,0001,0010,0011,1111011}}</code> é a indexadora de ''L0'', a ser utilizada nas funções de encode. Conforme o tipo de resolução é possível ignorar as fantasmas (tamanho maior que 4 bits) ou utilizá-las como as demais, com o cuidado de não subdividir antes do nível ''L1.5''. | |||
:<small>'''NOTA TÉCNICA'''. A rigor, por terem duas células de cobertura associadas ao quadrante ''IJ=02'', a resolução seria em duas etapas: usando "{0,1,2,3,4,5,6,7,8,9,a,b,c,d,e,.,.,.}" na cobertura, o resultado seria <code>IJ_to_cover={{NULL,NULL,"",NULL,NULL}, ...,{0000,0001,0010,0011,""}}</code>, com NULLs e bit strings vazias, "". Estabelecendo com "" a condição para resolução secundária, mas isso traria maior complexidade, mais fácil por hora ignorar a célula "fY" do mar territorial. </small> <!-- Os NULLs indicam que a célula da matriz não é usada. Os valores "" indicam que a resolução será por ''grid_l0_ghost'': é a condição para a verificação secundária, esses quadrantes no nível ''L0'' não possuem uma resolução direta, apenas no nível ''L1.5''. --> | |||
Temos por fim a seguinte '''função de resolução ''L0''''': | |||
<syntaxhighlight lang="sql" style="font-size: 80%;"> | |||
DROP FUNCTION grid_br.IJ0_to_L0(int,int,boolean) | |||
; | |||
DROP FUNCTION grid_br.IJ0_to_L0(int[],boolean) | |||
; | |||
CREATE FUNCTION grid_br.IJ0_to_L0(i int, j int, ignore_ghost boolean default true) RETURNS varbit AS $f$ | |||
SELECT CASE WHEN ignore_ghost AND bit_length(cbits)>4 THEN NULL ELSE cbits END | |||
FROM ( | |||
SELECT ('{{NULL,NULL,1111101,NULL,NULL},{NULL,1100,1101,1110,NULL},{NULL,1001,1010,1011,1111010},{0100,0101,0110,0111,1000},{0000,0001,0010,0011,1111011}}'::varbit[])[i+1][j+1] | |||
) t(cbits) | |||
$f$ LANGUAGE SQL IMMUTABLE; | |||
CREATE FUNCTION grid_br.IJ0_to_L0(ij int[], ignore_ghost boolean default true) RETURNS varbit AS $wrap$ | |||
SELECT grid_br.IJ0_to_L0($1[1], $1[2], $2) | |||
$wrap$ LANGUAGE SQL IMMUTABLE; | |||
</syntaxhighlight> | |||
Portanto <code>grid_br.IJ0_to_L0( grid_br.xyS_collapseTo_ijS(x,y) )</code> determina a cobertura L0. Com ''cbits'' determinamos quais valores Xmax e Ymax subtrair [[osmc:Metodologia/Algoritmo SQL/Issue04|cálculo GGeohash de uma célula de nível superior a zero]]. | |||
<!-- | |||
Situações: | |||
* mais geral, sem nação definida: ponto sem projeção, dado em latitude-longitude. O [[osmc:Metodologia#Seletor_de_jurisdição|seletor de jurisdição baseado em BBOX WGS84]] que determina a projeção a usar. | |||
* menos exigente, com ponto já na projeção nacional porém sem exigência de saber se fora ou dentro do território: aplicativos dedicados não esperam receber pontos exteriores, portanto não carece conferir se, mesmo em células ''L0'' de borda, é ou não interior. O seletor de célula se reduz ao BBOX da projeção, dispensa uso do polígono. | |||
* mais exigente: em aplicações onde certeza de ser ponto interior é necessária. Aí entra o uso do flag, pois se for selecionada uma célula ''L0'' com flag de borda, o seu contorno poligonal também precisará ser verificado. Tipicamente 3 passos na comparação do ponto ''pt'' com a célula ''cell'': <code>pt && cell AND (not_border OR ST_Contains(pt,cell))</code>. A otimização desse tipo de query requer de indexador específico. | |||
São esperadas [[osmc:Metodologia#Projeção_e_cobertura_nacionais|potências de 2]] nas células cheias de ''is_contained''. Para o Brasil é esperada a potência 20 no tamanho de lado do quadrado: <math>\sqrt{1099512\ km^2} = 1048576 \ m = 2^{20}</math>. Para Camarões (CM) potência 18: <math>\sqrt{68719\ km^2} = 262143.1 \ m \approx 262144 \ m = 2^{18}</math>. | |||
--> | |||
== Gerando a grade L0 de cobertura do país == | |||
Existem duas possibilidades de trabalho com a grade: | |||
# Com seletores e algoritmos de ''encode''/''decode''. | |||
# Com geometria das grades de células maiores, partindo de ''L0'', e depois, acima de certo nível, usando funções de geração de grade virtual. | |||
Como vimos, as configurações são transferidas para as mesmas tabelas onde foram geradas e arquivadas as geometrias L0. Podemos simplesmente recuperá-las. Tabela atual ''osmc.coverage'': | |||
<pre> | |||
Column | Type | Collation | Nullable | Default | |||
---------------+-------------+-----------+----------+--------- | |||
cbits | bit varying | | | | |||
isolabel_ext | text | | | | |||
cindex | text | | | | |||
bbox | integer[] | | | | |||
status | smallint | | | 0 | |||
is_country | boolean | | | false | |||
is_contained | boolean | | | false | |||
is_overlay | boolean | | | false | |||
kx_prefix | text | | | | |||
geom | geometry | | | | |||
geom_srid4326 | geometry | | | | |||
</pre> | |||
<syntaxhighlight lang="sql"> | |||
DROP view if exists report.v002_osmc_coverage_l0_geoms | |||
; | |||
CREATE VIEW report.v002_osmc_coverage_l0_geoms AS | |||
select row_number() OVER (ORDER BY cbits) AS gid, *, | |||
isolabel_ext ||'+'||cbits_b16 as afacode | |||
from ( | |||
select isolabel_ext,cbits, | |||
natcod.vbit_to_baseh(substring(cbits,11)::bit varying, 16) AS cbits_b16, | |||
geom_srid4326 | |||
from osmc.coverage where is_country and isolabel_ext IN ('BR','CO','CM') | |||
) t | |||
; | |||
</syntaxhighlight>Resultado da view no QGIS: | |||
[[Arquivo:Report v002.png|centro|semmoldura|620px]] | |||
===Sem recorte=== | |||
Pelo AFAcode ou sua representação binária (coluna ''cbits'' de ''osmc.coverage''), a célula JSON é obtida através de uma função pública, <code>api.osmcode_decode_scientific_absolute(codigo,'BR',18)</code>. Para fornecer a célula simples ela foi reescrita como <code>osmc.decode_scientific_absolute_geoms(codigo,'BR',18)</code> (retorna tabela PostGIS e portanto não vale como API). | |||
[[Arquivo:BRgrid-L0-QGISv1.png|centro|semmoldura|620px]] | |||
Funções ''osmc.decode_scientific_absolute_geoms'' e ''osmc.L0cover_br_geoms'' em https://github.com/osm-codes/GGeohash/blob/main/src/step03def-lib.sql#L1520 | |||
===Recorte L1 scientifica === | |||
Com os geradores abaixo obtemos a cobetrura de células com ~524 km de lado, que permite a comparação com as células originais do IBGE, conforme [[osmc:BR#Grade_de_cobertura]]. | |||
[[Arquivo:GridAFA-sci-BR-L1.png|centro|semmoldura|620x620px]] | |||
<syntaxhighlight lang="sql"> | |||
-- permanece "f*" por ser menor que hmrv: | |||
-- osmc.decode_scientific_absolute_geoms('fT,fY,fP,fN', 'BR', 18) | |||
-- geração da sequencia completa das demais coberturas candidatas: | |||
select array_agg(i||j) from unnest('{0,1,2,3,4,5,6,7,8,9,a,b,c,d,e}'::text[]) t1(i), unnest('{h,m,r,v}'::text[]) t2(j); | |||
-- Visualização no QGIS: | |||
drop table if exists tmp_br_cover_L1_scientific | |||
; | |||
create table tmp_br_cover_L1_scientific AS | |||
select row_number() over() as gid, * | |||
from ( | |||
SELECT 76::bit(10) || (natcod.baseh_to_vbit('0'||lower(code),16)) as cbits, | |||
code as b16_label, | |||
'0'||lower(code) as prefix, geom | |||
FROM osmc.decode_scientific_absolute_geoms( | |||
'0h,1h,2h,3h,4h,5h,6h,7h,8h,9h,ah,bh,ch,dh,eh,0m,1m,2m,3m,4m,5m,6m,7m,8m,9m,am,bm,cm,dm,em,0r,1r,2r,3r,4r,5r,6r,7r,8r,9r,ar,br,cr,dr,er,0v,1v,2v,3v,4v,5v,6v,7v,8v,9v,av,bv,cv,dv,ev' | |||
,'BR',18) | |||
where st_intersects(geom, st_transform((select geom from optim.jurisdiction_geom where isolabel_ext='BR'),952019) ) | |||
UNION ALL | |||
SELECT 76::bit(10) || (natcod.baseh_to_vbit('0'||lower(code),16)) as cbits, | |||
code as b16_label, | |||
'0'||lower(code) as prefix, geom | |||
FROM osmc.decode_scientific_absolute_geoms('fT,fY,fP,fN','BR',18) | |||
order by 1 | |||
) t3 | |||
; | |||
</syntaxhighlight> | |||
===Recorte L2 scientifica === | |||
Com os geradores abaixo obtemos a cobetrura de hexadecimais de 2 dígitos (células de ~262 km de lado). | |||
[[Arquivo:GridAFA-sci-BR-L2.png|centro|semmoldura|620x620px]] | |||
<syntaxhighlight lang="sql"> | |||
----- | |||
-- Geração da cobertura de dois dígitos hexa: | |||
-- fe,fa,f6,f4,ff,fb,f7,f5 obtidos de: | |||
SELECT array_agg(natcod.vbit_to_str(substring(cbits||x,15),'16h')) as b16_labels | |||
FROM ( | |||
SELECT (18+(row_number() over()))*16 as gid, x, | |||
76::bit(10) || (natcod.baseh_to_vbit('0'||lower(code),16)) as cbits | |||
FROM osmc.decode_scientific_absolute_geoms('fT,fY,fP,fN', 'BR', 18), | |||
LATERAL unnest(array[b'0',b'1']) t1(x) | |||
) t2 | |||
; | |||
-- geração da sequencia completa das demais coberturas candidatas: | |||
select array_agg(i||j) from unnest('{0,1,2,3,4,5,6,7,8,9,a,b,c,d,e}'::text[]) t1(i), unnest('{0,1,2,3,4,5,6,7,8,9,a,b,c,d,e,f}'::text[]) t2(j); | |||
-- Visualização no QGIS: | |||
drop table if exists tmp_br_cover_L2_scientific | |||
; | |||
create table tmp_br_cover_L2_scientific AS | |||
select row_number() over() as gid, * | |||
from ( | |||
SELECT 76::bit(10) || (natcod.baseh_to_vbit('0'||lower(code),16)) as cbits, | |||
code as b16_label, | |||
'0'||lower(code) as prefix, geom | |||
FROM osmc.decode_scientific_absolute_geoms( | |||
'00,01,02,03,04,05,06,07,08,09,0a,0b,0c,0d,0e,0f,10,11,12,13,14,15,16,17,18,19,1a,1b,1c,1d,1e,1f,20,21,22,23,24,25,26,27,28,29,2a,2b,2c,2d,2e,2f,30,31,32,33,34,35,36,37,38,39,3a,3b,3c,3d,3e,3f,40,41,42,43,44,45,46,47,48,49,4a,4b,4c,4d,4e,4f,50,51,52,53,54,55,56,57,58,59,5a,5b,5c,5d,5e,5f,60,61,62,63,64,65,66,67,68,69,6a,6b,6c,6d,6e,6f,70,71,72,73,74,75,76,77,78,79,7a,7b,7c,7d,7e,7f,80,81,82,83,84,85,86,87,88,89,8a,8b,8c,8d,8e,8f,90,91,92,93,94,95,96,97,98,99,9a,9b,9c,9d,9e,9f,a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,aa,ab,ac,ad,ae,af,b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,ba,bb,bc,bd,be,bf,c0,c1,c2,c3,c4,c5,c6,c7,c8,c9,ca,cb,cc,cd,ce,cf,d0,d1,d2,d3,d4,d5,d6,d7,d8,d9,da,db,dc,dd,de,df,e0,e1,e2,e3,e4,e5,e6,e7,e8,e9,ea,eb,ec,ed,ee,ef' | |||
,'BR',18) | |||
where st_intersects(geom, st_transform((select geom from optim.jurisdiction_geom where isolabel_ext='BR'),952019) ) | |||
UNION ALL | |||
SELECT 76::bit(10) || (natcod.baseh_to_vbit('0'||lower(code),16)) as cbits, | |||
code as b16_label, | |||
'0'||lower(code) as prefix, geom | |||
FROM osmc.decode_scientific_absolute_geoms('fe,fa,f6,f4,ff,fb,f7,f5','BR',18) | |||
order by 1 | |||
) t3 | |||
; | |||
</syntaxhighlight> | |||
==Ver também== | |||
* [[osmc:Metodologia/Algoritmo SQL/Lib]] | |||
* Países e respectivos SQLs: [[osmc:BR]]/[[osmc:BR/SQL|SQL]], [[osmc:CO]]/[[osmc:CO/SQL|SQL]], [[osmc:CM]]/[[osmc:CM/SQL|SQL]]. | |||
* [[osmc:Metodologia/Algoritmo SQL/Issues]] | |||
* DNGS - [[Discrete National Grid Systems/pt]] | |||
* [[osmc:Api]] | |||
------ | |||
Outras dicas: | |||
<syntaxhighlight lang="sql"> | <syntaxhighlight lang="sql"> | ||
-- Simulando a chamada online da cobertura do município de Boa Vista: | |||
select i, natcod.vbit_to_str(i::bit(5),'32nvu') as x_b32nvu, | |||
substr( api.osmcode_decode_postal( natcod.vbit_to_str(i::bit(5),'32nvu') ,'BR-RR-BoaVista')::text, 1,120) as geojson_cit120 | |||
from generate_series(0,31) t(i); | |||
</syntaxhighlight> | </syntaxhighlight> |
Edição atual tal como às 05h04min de 2 de julho de 2024
Detalhamento metodológico, descrevendo os algoritmos implementados em linguagem SQL. Para obtenção do código-fonte:
- Versão desejada se homologada: git.
- Versão corrente do servidor: se em dúvida, usar
pg_dump --schema-only -f dl05_dump_test.sql dl05
.
Configuração do país
Depois de percorrido todo o roteiro de decisões soberanas, o país está apto a implementar seu geocódigo oficial. As decisões são resumidas num arquivo YAML. Exemplos:
- Configurações do Brasil https://git.osm.codes/BR_new/blob/main/conf.yaml#L34
- Configurações de Camarões https://git.osm.codes/CM/blob/main/conf.yaml#L37
- Configurações da Colômbia https://git.osm.codes/CO_new/blob/main/conf.yaml#L34
Os parâmetros principais do YAML da Colômbia, por exemplo, são:
dngsCode_definition: srid: 102022 grid_xy_origin: 3678500 970000 grid_cell_side: 524288 grid_l0_cell_ij: 02 03 10 11 12 13 20 21 22 23 30 31 32 40 41 42 grid_l0_cell_idxb16: 8 a 1 3 9 b 4 6 c e 5 7 d 0 2 f
onde:
- srid é a projeção igual-area do país;
- grid_xy_origin: é a origem no plano projetado;
- grid_cell_side: é o tamanho (potência de 2) da célula L0 adotada;
- grid_l0_cell: a lista de coordenadas ij da matriz de células, indicando aquelas selecionados para a cobertura;
- grid_l0_cell_idxb16: a sua correspondente indexação, já expressa em base16.
A parametrização YAML poderia ser lida diretamente do arquivo (pendente implementar essa estratégia); todavia, atualmente, seus parâmetros estão embutidos na função wrapper osmc.L0cover_upsert
.
A configuração é efetivada pela chamada com o país específico, por exemplo Brasil SELECT osmc.L0cover_upsert('BR');
A função osmc.L0_upsert()
atualiza a tabela osmc.coverage com linhas is_country e sem popular as colunas cindex, is_overlay ou kx_prefix (utilizadas apenas por coberturas municipais). A tabela abaixo, gerada pelo relatório osmc_report.v001_osmc_coverage_l0_list apresentado em seguida, ilustra os casos típicos do Brasil, Camarões e Colômbia:
pais | is_contained | cbits | b16 | area_km2 |
---|---|---|---|---|
BR = 1 = 00000001 | f | 00000001.00000000 | 00 | 275241 |
BR = 1 = 00000001 | f | 00000001.00000001 | 01 | 757637 |
BR = 1 = 00000001 | f | 00000001.00000010 | 02 | 645014 |
... | ||||
BR = 1 = 00000001 | f | 00000001.00010000 | 10 | 1519 |
BR = 1 = 00000001 | f | 00000001.00010001 | 11 | 3618 |
BR = 1 = 00000001 | t | 00000001.00000110 | 06 | 1099512 |
CM = 3 = 00000011 | f | 00000011.0001 | 1 | 16259 |
CM = 3 = 00000011 | f | 00000011.0010 | 2 | 1261 |
CM = 3 = 00000011 | f | 00000011.0011 | 3 | 43764 |
... | ||||
CM = 3 = 00000011 | f | 00000011.1110 | e | 49039 |
CM = 3 = 00000011 | f | 00000011.1111 | f | 14919 |
CM = 3 = 00000011 | t | 00000011.1010 | a | 68719 |
CO = 2 = 00000010 | f | 00000010.0000 | 0 | 8083 |
CO = 2 = 00000010 | f | 00000010.0001 | 1 | 230987 |
CO = 2 = 00000010 | f | 00000010.0010 | 2 | 57590 |
... | ||||
CO = 2 = 00000010 | f | 00000010.1101 | d | 196017 |
CO = 2 = 00000010 | f | 00000010.1110 | e | 84471 |
CO = 2 = 00000010 | f | 00000010.1111 | f | 61264 |
A coluna cbits é composta de duas partes, o código do país e o código da sua célula L0. O código é traduzido em hexadecimal (base16) na coluna b16, e suas células formam a cobertura do país, ou seja, a grade de nível zero. Nos geocódigos de notação científica apare cerão sempre como primeiro dígito.
Os códigos dos países são números inteiros de 8 bits (atualmente, sem compromisso com códigos ISO de 10 bits), garantindo melhor ocupação nos 57 bits (56 bits para acomodar 14 dígitos hexadecimais) disponíveis de geocódigo hInt64 (essa representação alternativa é mais eficiente e versátil do que o tipo de dado varbit (bit string) do PostgreSQL.
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 está contida no polígono do país. Ver Issue 02 para detalhes sobre a dedução do caso cheio pela área.
CREATE VIEW report.v001_osmc_coverage_l0_list AS
SELECT isolabel_ext ||' = '||(b'0000000000000000000000'||cb10)::bit(32)::int ||' = '||cb10 as pais,
is_contained,
cb10::text||'.'||cb_suffix as cbits, natcod.vbit_to_baseh(cb_suffix,16) as b16,
round(st_area(geom)/1000.0^2) as area_km2
FROM (
select *, substring(cbits,1,10) as cb10, substring(cbits,11) as cb_suffix
from osmc.coverage
) t
WHERE is_country ORDER BY 1,2,cbits
;
Seletor de projeção
O padrão DNGS, conforme seu algoritmo seletor de jurisdição, requer que qualquer ponto sobre o globo, expresso em WGS84 sem projeção, tenha seu geocódigo resolvido conforme a jurisdição da sua posição. O algoritmo decide qual projeção usar. A verificação em grandes BBOXes, sem bordas, tem alta performance, enquanto que a decisão sobre polígonos de borda tem menor performance.
A grade L0 serve de referência para a criação de recortes BBOX WGS84 sobre o globo, onde as escolhas de projeção são reproduzidas de forma mais grosseira que o mosaico L0, e, nas bordas, recortadas. Supondo pt a geometria de um ponto:
SELECT srid FROM osmc.coverage_srid
WHERE pt && geom
AND (is_not_border OR ST_Intersects(pt,geom))- indexação prévia com
CREATE INDEX osmc_coverage_srid_idx1
ON osmc.coverage_srid(is_not_border,geom)
USING GIST (geom) - a performance pode ser avaliada por
VACUUM ANALYZE osmc_coverage_srid_idx1
O conhecimento prévio do SRID (pode ser por contexto) permite por fim fazer a busca na tabela osmc.coverage
.
Seletor de jurisdição e L0
Uma vez selecionada a projeção, o ponto é transformado conforme SRID e confrontado apenas contra as BBOXes de células L0 daquela projeção. No caso especial de ser uma projeção adotada por mais de um país (por exemplo Itália e Vaticano), e o ponto cair em uma interseção de BBOXes, a seleção do recorte de L0 correto dependerá do ponto estar contido no polígono do recorte.
Justamente para contemplar os casos especiais a seleção de jurisdição acaba usando uma condição e indexação mais complexa, similar à seletora de projeção:
SELECT cbits as cbits_l0 FROM osmc.coverage
WHERE srid=$srid AND pt && geom
AND (is_not_border OR ST_Intersects(pt,geom))- indexação prévia com
CREATE INDEX osmc_coverage_idx1
ON osmc.coverage(srid,is_not_border,geom)
USING GIST (geom)
O valor cbits_l0 é o identificador de célula L0 do país correto, e o identificador de país é seu prefixo. A partir do ponto pt (com geometria já no SRID do país) e cbits_l0 obtemos por algoritmo GGeohash o valor cbits_pt da célula pontual, cbits_pt=encode(pt,cbits_l0)
.
NOTA. Em aplicações que não partem de um cenário global, ou de pontos que não precisam ser selecionados como interiores ou não ao território nacional, pode-se adotar um seletor mais simples: grid_br.xyS_collapseTo_ijS(pt)
por default retorna as coordenadas IJ de L0.
SELECT cbits as cbits_l0 FROM osmc.coverage
WHERE srid=$srid AND pt && geom- indexação prévia com
CREATE INDEX osmc_coverage_idx1
ON osmc.coverage(srid,geom)
USING GIST (geom)
Por fim, como o país não pode ter mais do que 16 células L0, uma cadeia de condições IF ou CASE é suficiente e provavelmente mais rápido do que consultar o banco de dados. Usando ST_X(pt)
e ST_Y(pt)
e os valores de BBOX, podemos otimizar o processo de escolha de L0 e encode do ponto.
Tratamento das configurações
Usando o exemplo do Brasil, "grid_lc" indica que a cobertura é uma matriz de 5x5 (linhas por colunas), e "grid_l0_cell" quais são as células IJ selecionadas pelo índice "grid_l0_cell_idx":
grid_lc: 5 5 grid_l0_cell: 40 41 42 43 30 31 32 33 34 21 22 23 11 12 13 02 44 24 grid_l0_cell_idx: 0 1 2 3 4 5 6 7 8 9 a b c d e fT fP fN
As configurações podem ser reescritas com auxilio de uma query:
WITH m AS ( -- background matrix
select i,j from generate_series(0,4) t1(i), generate_series(0,4) t2(j)
), c AS ( -- cover cells
select i,j, (natcod.baseh_to_vbit('{0,1,2,3,4,5,6,7,8,9,a,b,c,d,e,fT,fP,fN}'::text[],16))[(p+2)/3] as cbits
from (
select i,j, strpos('40 41 42 43 30 31 32 33 34 21 22 23 11 12 13 02 44 24', i::text||j::text)::float as p
from m
) tt where p>0
) SELECT array_agg(a order by i) -- IJ_to_cover matrix
FROM (
SELECT m.i, array_agg( c.cbits order by m.j) a
FROM c RIGHT JOIN m ON c.i=m.i AND c.j=m.j
GROUP BY 1 ORDER BY 1
) t4
;
A matriz resultante IJ_to_cover={{NULL,NULL,1111101,NULL,NULL}, ..., {0000,0001,0010,0011,1111011}}
é a indexadora de L0, a ser utilizada nas funções de encode. Conforme o tipo de resolução é possível ignorar as fantasmas (tamanho maior que 4 bits) ou utilizá-las como as demais, com o cuidado de não subdividir antes do nível L1.5.
- NOTA TÉCNICA. A rigor, por terem duas células de cobertura associadas ao quadrante IJ=02, a resolução seria em duas etapas: usando "{0,1,2,3,4,5,6,7,8,9,a,b,c,d,e,.,.,.}" na cobertura, o resultado seria
IJ_to_cover={{NULL,NULL,"",NULL,NULL}, ...,{0000,0001,0010,0011,""}}
, com NULLs e bit strings vazias, "". Estabelecendo com "" a condição para resolução secundária, mas isso traria maior complexidade, mais fácil por hora ignorar a célula "fY" do mar territorial.
Temos por fim a seguinte função de resolução L0:
DROP FUNCTION grid_br.IJ0_to_L0(int,int,boolean)
;
DROP FUNCTION grid_br.IJ0_to_L0(int[],boolean)
;
CREATE FUNCTION grid_br.IJ0_to_L0(i int, j int, ignore_ghost boolean default true) RETURNS varbit AS $f$
SELECT CASE WHEN ignore_ghost AND bit_length(cbits)>4 THEN NULL ELSE cbits END
FROM (
SELECT ('{{NULL,NULL,1111101,NULL,NULL},{NULL,1100,1101,1110,NULL},{NULL,1001,1010,1011,1111010},{0100,0101,0110,0111,1000},{0000,0001,0010,0011,1111011}}'::varbit[])[i+1][j+1]
) t(cbits)
$f$ LANGUAGE SQL IMMUTABLE;
CREATE FUNCTION grid_br.IJ0_to_L0(ij int[], ignore_ghost boolean default true) RETURNS varbit AS $wrap$
SELECT grid_br.IJ0_to_L0($1[1], $1[2], $2)
$wrap$ LANGUAGE SQL IMMUTABLE;
Portanto grid_br.IJ0_to_L0( grid_br.xyS_collapseTo_ijS(x,y) )
determina a cobertura L0. Com cbits determinamos quais valores Xmax e Ymax subtrair cálculo GGeohash de uma célula de nível superior a zero.
Gerando a grade L0 de cobertura do país
Existem duas possibilidades de trabalho com a grade:
- Com seletores e algoritmos de encode/decode.
- Com geometria das grades de células maiores, partindo de L0, e depois, acima de certo nível, usando funções de geração de grade virtual.
Como vimos, as configurações são transferidas para as mesmas tabelas onde foram geradas e arquivadas as geometrias L0. Podemos simplesmente recuperá-las. Tabela atual osmc.coverage:
Column | Type | Collation | Nullable | Default ---------------+-------------+-----------+----------+--------- cbits | bit varying | | | isolabel_ext | text | | | cindex | text | | | bbox | integer[] | | | status | smallint | | | 0 is_country | boolean | | | false is_contained | boolean | | | false is_overlay | boolean | | | false kx_prefix | text | | | geom | geometry | | | geom_srid4326 | geometry | | |
DROP view if exists report.v002_osmc_coverage_l0_geoms
;
CREATE VIEW report.v002_osmc_coverage_l0_geoms AS
select row_number() OVER (ORDER BY cbits) AS gid, *,
isolabel_ext ||'+'||cbits_b16 as afacode
from (
select isolabel_ext,cbits,
natcod.vbit_to_baseh(substring(cbits,11)::bit varying, 16) AS cbits_b16,
geom_srid4326
from osmc.coverage where is_country and isolabel_ext IN ('BR','CO','CM')
) t
;
Resultado da view no QGIS:
Sem recorte
Pelo AFAcode ou sua representação binária (coluna cbits de osmc.coverage), a célula JSON é obtida através de uma função pública, api.osmcode_decode_scientific_absolute(codigo,'BR',18)
. Para fornecer a célula simples ela foi reescrita como osmc.decode_scientific_absolute_geoms(codigo,'BR',18)
(retorna tabela PostGIS e portanto não vale como API).
Funções osmc.decode_scientific_absolute_geoms e osmc.L0cover_br_geoms em https://github.com/osm-codes/GGeohash/blob/main/src/step03def-lib.sql#L1520
Recorte L1 scientifica
Com os geradores abaixo obtemos a cobetrura de células com ~524 km de lado, que permite a comparação com as células originais do IBGE, conforme osmc:BR#Grade_de_cobertura.
-- permanece "f*" por ser menor que hmrv:
-- osmc.decode_scientific_absolute_geoms('fT,fY,fP,fN', 'BR', 18)
-- geração da sequencia completa das demais coberturas candidatas:
select array_agg(i||j) from unnest('{0,1,2,3,4,5,6,7,8,9,a,b,c,d,e}'::text[]) t1(i), unnest('{h,m,r,v}'::text[]) t2(j);
-- Visualização no QGIS:
drop table if exists tmp_br_cover_L1_scientific
;
create table tmp_br_cover_L1_scientific AS
select row_number() over() as gid, *
from (
SELECT 76::bit(10) || (natcod.baseh_to_vbit('0'||lower(code),16)) as cbits,
code as b16_label,
'0'||lower(code) as prefix, geom
FROM osmc.decode_scientific_absolute_geoms(
'0h,1h,2h,3h,4h,5h,6h,7h,8h,9h,ah,bh,ch,dh,eh,0m,1m,2m,3m,4m,5m,6m,7m,8m,9m,am,bm,cm,dm,em,0r,1r,2r,3r,4r,5r,6r,7r,8r,9r,ar,br,cr,dr,er,0v,1v,2v,3v,4v,5v,6v,7v,8v,9v,av,bv,cv,dv,ev'
,'BR',18)
where st_intersects(geom, st_transform((select geom from optim.jurisdiction_geom where isolabel_ext='BR'),952019) )
UNION ALL
SELECT 76::bit(10) || (natcod.baseh_to_vbit('0'||lower(code),16)) as cbits,
code as b16_label,
'0'||lower(code) as prefix, geom
FROM osmc.decode_scientific_absolute_geoms('fT,fY,fP,fN','BR',18)
order by 1
) t3
;
Recorte L2 scientifica
Com os geradores abaixo obtemos a cobetrura de hexadecimais de 2 dígitos (células de ~262 km de lado).
-----
-- Geração da cobertura de dois dígitos hexa:
-- fe,fa,f6,f4,ff,fb,f7,f5 obtidos de:
SELECT array_agg(natcod.vbit_to_str(substring(cbits||x,15),'16h')) as b16_labels
FROM (
SELECT (18+(row_number() over()))*16 as gid, x,
76::bit(10) || (natcod.baseh_to_vbit('0'||lower(code),16)) as cbits
FROM osmc.decode_scientific_absolute_geoms('fT,fY,fP,fN', 'BR', 18),
LATERAL unnest(array[b'0',b'1']) t1(x)
) t2
;
-- geração da sequencia completa das demais coberturas candidatas:
select array_agg(i||j) from unnest('{0,1,2,3,4,5,6,7,8,9,a,b,c,d,e}'::text[]) t1(i), unnest('{0,1,2,3,4,5,6,7,8,9,a,b,c,d,e,f}'::text[]) t2(j);
-- Visualização no QGIS:
drop table if exists tmp_br_cover_L2_scientific
;
create table tmp_br_cover_L2_scientific AS
select row_number() over() as gid, *
from (
SELECT 76::bit(10) || (natcod.baseh_to_vbit('0'||lower(code),16)) as cbits,
code as b16_label,
'0'||lower(code) as prefix, geom
FROM osmc.decode_scientific_absolute_geoms(
'00,01,02,03,04,05,06,07,08,09,0a,0b,0c,0d,0e,0f,10,11,12,13,14,15,16,17,18,19,1a,1b,1c,1d,1e,1f,20,21,22,23,24,25,26,27,28,29,2a,2b,2c,2d,2e,2f,30,31,32,33,34,35,36,37,38,39,3a,3b,3c,3d,3e,3f,40,41,42,43,44,45,46,47,48,49,4a,4b,4c,4d,4e,4f,50,51,52,53,54,55,56,57,58,59,5a,5b,5c,5d,5e,5f,60,61,62,63,64,65,66,67,68,69,6a,6b,6c,6d,6e,6f,70,71,72,73,74,75,76,77,78,79,7a,7b,7c,7d,7e,7f,80,81,82,83,84,85,86,87,88,89,8a,8b,8c,8d,8e,8f,90,91,92,93,94,95,96,97,98,99,9a,9b,9c,9d,9e,9f,a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,aa,ab,ac,ad,ae,af,b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,ba,bb,bc,bd,be,bf,c0,c1,c2,c3,c4,c5,c6,c7,c8,c9,ca,cb,cc,cd,ce,cf,d0,d1,d2,d3,d4,d5,d6,d7,d8,d9,da,db,dc,dd,de,df,e0,e1,e2,e3,e4,e5,e6,e7,e8,e9,ea,eb,ec,ed,ee,ef'
,'BR',18)
where st_intersects(geom, st_transform((select geom from optim.jurisdiction_geom where isolabel_ext='BR'),952019) )
UNION ALL
SELECT 76::bit(10) || (natcod.baseh_to_vbit('0'||lower(code),16)) as cbits,
code as b16_label,
'0'||lower(code) as prefix, geom
FROM osmc.decode_scientific_absolute_geoms('fe,fa,f6,f4,ff,fb,f7,f5','BR',18)
order by 1
) t3
;
Ver também
- osmc:Metodologia/Algoritmo SQL/Lib
- Países e respectivos SQLs: osmc:BR/SQL, osmc:CO/SQL, osmc:CM/SQL.
- osmc:Metodologia/Algoritmo SQL/Issues
- DNGS - Discrete National Grid Systems/pt
- osmc:Api
Outras dicas:
-- Simulando a chamada online da cobertura do município de Boa Vista:
select i, natcod.vbit_to_str(i::bit(5),'32nvu') as x_b32nvu,
substr( api.osmcode_decode_postal( natcod.vbit_to_str(i::bit(5),'32nvu') ,'BR-RR-BoaVista')::text, 1,120) as geojson_cit120
from generate_series(0,31) t(i);