osmc:Metodologia/Algoritmo 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
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:
osmCode_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; e grid_l0_cell_sci_base 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 são copiadas/coladas na unção de inclusão desses dados na sua tabela. Exemplos:
- Inclusão da configuração do Brasil https://git.osm.codes/GGeohash/blob/main/src/step04def-ini.sql#L95
- Inclusão da configuração de Camarões https://git.osm.codes/GGeohash/blob/main/src/step04def-ini.sql#L96
- Inclusão da configuração da Colômbia https://git.osm.codes/GGeohash/blob/main/src/step04def-ini.sql#L98
A configuração é efetivada pela chamada com o país específico, por exemplo Brasil SELECT osmc.L0cover_upsert('BR');
A função 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 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 = 76 = 0001001100 | f | 0001001100.00000110 | 06 | 1099512 |
BR = 76 = 0001001100 | t | 0001001100.00000000 | 00 | 275241 |
BR = 76 = 0001001100 | t | 0001001100.00000001 | 01 | 757637 |
BR = 76 = 0001001100 | t | 0001001100.00000010 | 02 | 645014 |
... | ||||
BR = 76 = 0001001100 | t | 0001001100.00001110 | 0e | 54071 |
BR = 76 = 0001001100 | t | 0001001100.00001111 | 0f | 17440 |
BR = 76 = 0001001100 | t | 0001001100.00010000 | 10 | 1519 |
BR = 76 = 0001001100 | t | 0001001100.00010001 | 11 | 3618 |
CM = 120 = 0001111000 | f | 0001111000.1001 | 9 | 68719 |
CM = 120 = 0001111000 | t | 0001111000.0001 | 1 | 2469 |
CM = 120 = 0001111000 | t | 0001111000.0010 | 2 | 24181 |
... | ||||
CM = 120 = 0001111000 | t | 0001111000.1101 | d | 47183 |
CM = 120 = 0001111000 | t | 0001111000.1110 | e | 5164 |
CO = 170 = 0010101010 | t | 0010101010.0000 | 0 | 8083 |
CO = 170 = 0010101010 | t | 0010101010.0001 | 1 | 141953 |
CO = 170 = 0010101010 | t | 0010101010.0010 | 2 | 57590 |
... | ||||
CO = 170 = 0010101010 | t | 0010101010.1110 | e | 84471 |
CO = 170 = 0010101010 | t | 0010101010.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 de país nessa versão ainda eram códigos ISO ocupando 10 bits. Por exemplo Brasil 0001001100
é 76. Na versão mais atual são códigos internos (sem compromisso com ISO) de 8 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 tem interseção com as bordas 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 . . . grid_l0_ghost: {02:fT fY, 44:fP, 24: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
-- ... (natcod.baseh_to_vbit('{0,1,2,3,4,5,6,7,8,9,a,b,c,d,e,.,.,.}'::text[],16))
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
- DNGS - Discrete National Grid Systems/pt
- osmc:Metodologia/Algoritmo SQL/Issues
- osmc:Api
- Países e respectivos SQLs: osmc:BR/SQL, osmc:CO/SQL, osmc:CM/SQL.
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);