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