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, ao invés de serem lidos diretamente do arquivo (pendente implementar essa estratégia), 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');
Deveria se chamar L0_upsert_get(contry)
. A função de que fato implementa as configurações é osmc.L0_upsert(). Importante notar que, dentro dela, temos um comando unnest(a,b)
que não deve ser confundido com o produto cruzado, o PostgreSQL alinha as duas arrays mesmo que de tamanhos diferentes, fazendo uma espécie de left join. Ilustrando com duas arrays reduzidas:
SELECT count(*) n, count(distinct prefix) n_prefix -- n=5 | n_prefix=3
FROM unnest('{40,41,42}'::int[],'{00,01,02,0a,0b}'::text[]) t(prefix,quadrant);
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 consulta abaixo ilustra os casos típicos do Brasil e Colômbia.
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,
status, round(st_area(geom)/1000000.0) 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
;
pais | is_contained | cbits | b16 | area_km2 | status |
---|---|---|---|---|---|
BR = 76 = 0001001100 | f | 0001001100.00000110 | 06 | 1099512 | 1 |
BR = 76 = 0001001100 | t | 0001001100.00000000 | 00 | 275241 | 1 |
BR = 76 = 0001001100 | t | 0001001100.00000001 | 01 | 757637 | 1 |
BR = 76 = 0001001100 | t | 0001001100.00000010 | 02 | 645014 | 1 |
... | |||||
BR = 76 = 0001001100 | t | 0001001100.00001110 | 0e | 54071 | 1 |
BR = 76 = 0001001100 | t | 0001001100.00001111 | 0f | 17440 | 1 |
BR = 76 = 0001001100 | t | 0001001100.00010000 | 10 | 1519 | 1 |
BR = 76 = 0001001100 | t | 0001001100.00010001 | 11 | 3618 | 1 |
CM = 120 = 0001111000 | f | 0001111000.1001 | 9 | 68719 | 1 |
CM = 120 = 0001111000 | t | 0001111000.0001 | 1 | 2469 | 1 |
CM = 120 = 0001111000 | t | 0001111000.0010 | 2 | 24181 | 1 |
... | |||||
CM = 120 = 0001111000 | t | 0001111000.1101 | d | 47183 | 1 |
CM = 120 = 0001111000 | t | 0001111000.1110 | e | 5164 | 1 |
CO = 170 = 0010101010 | t | 0010101010.0000 | 0 | 8083 | 1 |
CO = 170 = 0010101010 | t | 0010101010.0001 | 1 | 141953 | 1 |
CO = 170 = 0010101010 | t | 0010101010.0010 | 2 | 57590 | 1 |
... | |||||
CO = 170 = 0010101010 | t | 0010101010.1110 | e | 84471 | 1 |
CO = 170 = 0010101010 | t | 0010101010.1111 | f | 61264 | 1 |
A coluna cbits é composta de duas partes, o código do país e o código da sua célula L0 (traduzido em hexadecimal=base16). 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. 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: .
- PS: "bug de nome de coluna" a corrigir, a coluna is_contained deveria indicar que a célula está totalmente contida no polígono do país, mas está sendo usado o seu NOT.
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).
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);