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: 02 03 10 11 12 13 20 21 22 23 30 31 32 40 41 42
  grid_l0_cell_sci_base: 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 aparecerão sempre como primeiro dígito.

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 jurisdição e L0

Em aplicações que partem de um cenário global, ou de pontos que precisam ser selecionados como interiores ou não ao território nacional, é importante detectar não só a célula L0 mas também o polígono territorial. Situações:

  • mais geral, sem nação definida: ponto sem projeção, dado em latitude-longitude. 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: pt && cell AND (not_border OR ST_Contains(pt,cell)). A otimização desse tipo de query requer de indexador específico.


São esperadas 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: . Para Camarões (CM) potência 18: .

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

Gerando a grade L0 de cobertura do país

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 
---------------+-------------+-----------+----------+---------
          | 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
-- FUNCAO MIGROU PARA O GIT, APAGAR!
CREATE FUNCTION osmc.decode_scientific_absolute_geoms(
   p_code text, -- um ou mais (separados por virgula) afaCodes científicos separados por virgula
   p_iso text,  -- pais de contextualização do afaCode.
   p_base integer DEFAULT 18,  -- detecta antes se usa gambiarra se falsa célula ... não devia precisar.
   p_use_resilient boolean DEFAULT true
) RETURNS TABLE (
	code text,
	area real,
	side real,
	truncated_code text,
	base text,
	geom geometry,
	geom4326 geometry
) language SQL IMMUTABLE
AS $f$
    SELECT
        TRANSLATE(code_tru,'gqhmrvjknpstzy','GQHMRVJKNPSTZY') as code,
        ST_Area(v.geom) as area,
        SQRT(ST_Area(v.geom)) as side,
        truncated_code,
        osmc.string_base(p_base) as base,
        v.geom,
        CASE WHEN p_use_resilient THEN ST_Transform_resilient(v.geom,4326,0.005) ELSE ST_Transform(v.geom,4326) END as geom4326
    FROM (
      SELECT DISTINCT code16h,

      -- trunca
      CASE
        WHEN p_base <> 18 AND length(code16h) > 12 AND up_iso IN ('BR')           THEN substring(code16h,1,12)
        WHEN p_base <> 18 AND length(code16h) > 11 AND up_iso IN ('EC','CO','UY') THEN substring(code16h,1,11)
        WHEN p_base <> 18 AND length(code16h) > 10 AND up_iso IN ('CM')           THEN substring(code16h,1,10)
        WHEN p_base =  18 AND length(code)    > 11 AND up_iso IN ('BR')           THEN substring(code,1,11)
        WHEN p_base =  18 AND length(code)    > 10 AND up_iso IN ('UY')           THEN substring(code,1,10)
        ELSE (CASE WHEN p_base=18 THEN code ELSE code16h END)
      END AS code_tru,

      -- flag
      CASE
        WHEN p_base <> 18 AND length(code16h) > 12 AND up_iso IN ('BR')           THEN TRUE
        WHEN p_base <> 18 AND length(code16h) > 11 AND up_iso IN ('EC','CO','UY') THEN TRUE
        WHEN p_base <> 18 AND length(code16h) > 10 AND up_iso IN ('CM')           THEN TRUE
        WHEN p_base =  18 AND length(code)    > 11 AND up_iso IN ('BR')           THEN TRUE
        WHEN p_base =  18 AND length(code)    > 10 AND up_iso IN ('UY')           THEN TRUE
        ELSE NULL
      END AS truncated_code,

      -- vbit code16h
      CASE
        WHEN length(code16h) > 12 AND up_iso IN ('BR')           THEN natcod.baseh_to_vbit(substring(code16h,1,12),16)
        WHEN length(code16h) > 11 AND up_iso IN ('EC','CO','UY') THEN natcod.baseh_to_vbit(substring(code16h,1,11),16)
        WHEN length(code16h) > 10 AND up_iso IN ('CM')           THEN natcod.baseh_to_vbit(substring(code16h,1,10),16)
        ELSE natcod.baseh_to_vbit(code16h,16)
      END AS codebits,

      code, up_iso

      FROM
      (
        SELECT code, upper(p_iso) AS up_iso,
                CASE
                  WHEN p_base = 18 THEN osmc.decode_16h1c(code,upper(p_iso))
                  ELSE code
                END AS code16h
        FROM regexp_split_to_table(lower(p_code),',') code
      ) u
    ) c,
    LATERAL
    (
      SELECT ggeohash.draw_cell_bybox(ggeohash.decode_box2(osmc.vbit_withoutL0(codebits,c.up_iso),bbox, CASE WHEN c.up_iso='EC' THEN TRUE ELSE FALSE END),false,ST_SRID(geom)) AS geom
      FROM osmc.coverage
      WHERE is_country IS TRUE AND isolabel_ext = c.up_iso -- cobertura nacional apenas
        AND
        CASE
        WHEN up_iso IN ('CO','CM') THEN ( ( osmc.extract_L0bits(cbits,'CO')   # codebits::bit(4) ) = 0::bit(4) ) -- 1 dígito base16h
        ELSE                            ( ( osmc.extract_L0bits(cbits,up_iso) # codebits::bit(8) ) = 0::bit(8) ) -- 2 dígitos base16h
        END
    ) v
    WHERE
    CASE WHEN up_iso = 'UY' THEN c.code16h NOT IN ('0eg','10g','12g','00r','12r','0eh','05q','11q') ELSE TRUE END
$f$;

DROP FUNCTION if exists osmc.L0cover_br_geoms
;
CREATE FUNCTION osmc.L0cover_br_geoms() 
RETURNS TABLE (
    gid  int,
	cbits         varbit,
	b16_label     text,
	prefix        text,
	isolabel_ext  text,
	bbox          integer[],
	is_contained  boolean,
	geom          geometry,
	geom_cell     geometry,
	geom_srid4326 geometry
) language SQL AS $f$
  SELECT row_number() over() as gid,
         jurisd_base_id::bit(10) || (natcod.baseh_to_vbit(prefix,16)) as cbits,
         CASE WHEN left(prefix,1)='0' THEN substring(prefix,2,1) ELSE prefix END as b16_label,
         prefix,-- natcod.vbit_to_baseh(substring(cbits,11),16) as b16_label,
         'BR', bbox,
         CASE WHEN ST_ContainsProperly(geom_country,geom_cell) IS FALSE THEN TRUE ELSE FALSE END,
         geom,
         geom_cell,
         ST_Transform(geom,4326)
  FROM (
      SELECT 76 AS jurisd_base_id, prefix, bbox,geom_country,
         ST_Intersection(ggeohash.draw_cell_bybox(bbox,false,952019),geom_country) AS geom,
         ggeohash.draw_cell_bybox(bbox,false,952019) AS geom_cell
      FROM unnest(
            '{00,01,02,03,04,05,06,07,08,09,0a,0b,0c,0d,0e,0f,10,11}'::text[],
            array[20,21,22,23,15,16,17,18,19,11,12,13,6,7,8,2,24,14]
      ) t(prefix,quadrant),
      LATERAL (
         SELECT osmc.ij_to_bbox(quadrant%5,quadrant/5,2715000,6727000,1048576)
      ) u(bbox),
      LATERAL (
         SELECT ST_Transform(geom,952019)
         FROM optim.vw01full_jurisdiction_geom g
         WHERE g.isolabel_ext = 'BR' AND jurisd_base_id = 76
      ) r(geom_country)
      WHERE quadrant IS NOT NULL
  ) y

   UNION ALL

   SELECT 18+(row_number() over()) as gid,
        76::bit(10) || (natcod.baseh_to_vbit('0'||lower(code),16)) as cbits,
        code as b16_label,
        '0'||lower(code)  as prefix,
        'BR', NULL as bbox,  false as is_contained,
        NULL as geom, -- st_intersect(geom_BR,geom_cell)
	geom as geom_cell, geom4326
   FROM osmc.decode_scientific_absolute_geoms('fT,fY,fP,fN','BR',18)

  ORDER BY 1
$f$;
COMMENT ON FUNCTION osmc.L0cover_br_geoms()
  IS 'L0cover BR from configs ingest (gid<=18) into osmc.coverage using osmc.l0cover_upsert().'
;

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=4009”