osmc:Metodologia/Algoritmo SQL

De Documentação

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:

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:

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

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:

Report v002.png

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).

BRgrid-L0-QGISv1.png

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.

GridAFA-sci-BR-L1.png
-- 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).

GridAFA-sci-BR-L2.png
-----
-- 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

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);
Disponível em “https://wiki.addressforall.org/w/index.php?title=osmc:Metodologia/Algoritmo_SQL&oldid=4069”