osmc:Metodologia/Algoritmo SQL

Detalhamento metodológico, descrevendo os algoritmos implementados em linguagem SQL. Para obtenção do código-fonte:

1. Versão desejada se homologada: git.
2. 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;
• 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 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

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 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:

1. Com seletores e algoritmos de encode/decode.
2. 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);