osmc:Metodologia/Algoritmo SQL/Lib: mudanças entre as edições
m (rev agg) |
|||
| Linha 123: | Linha 123: | ||
drop function if exists vbit_DeInterleave; | drop function if exists vbit_DeInterleave; | ||
; | DROP AGGREGATE if exists bitcat_agg(varbit); | ||
CREATE AGGREGATE bitcat_agg(varbit) ( | |||
STYPE = varbit | |||
,SFUNC = bitcat -- \do+ "||" | |||
); | |||
CREATE FUNCTION vbit_to_int(x varbit) returns int as $f$ | CREATE FUNCTION vbit_to_int(x varbit) returns int as $f$ | ||
SELECT ( substring(0::bit(32),bit_length(x)+1) || x )::bit(32)::int | SELECT ( substring(0::bit(32),bit_length(x)+1) || x )::bit(32)::int | ||
| Linha 137: | Linha 142: | ||
CREATE FUNCTION vbit_interleave(x varbit, y varbit) returns varbit as $f$ | CREATE FUNCTION vbit_interleave(x varbit, y varbit) returns varbit as $f$ | ||
SELECT | SELECT bitcat_agg( substring(x,i,1) || substring(y,i,1) ) | ||
FROM generate_series(1,bit_length(x)) t(i) | FROM generate_series(1,bit_length(x)) t(i) | ||
$f$ LANGUAGE SQL IMMUTABLE; | $f$ LANGUAGE SQL IMMUTABLE; | ||
| Linha 162: | Linha 167: | ||
Cada célula da cobertura ''L0'' já tem seu código de 4 bits (dígito hexadecimal), em seguida qualquer outra célula de qualquer nível hierárquico terá um código ''cbits'' do nível ''L'' com ''L*2'' bits obtidos das coordenadas IJ através da função <code>ints_to_interleavedbits(i,j)</code>. O tamanho de lado da célula do nível ''L'', no caso do Brasil, será <math>2^{20-L}</math> metros. A célula de 1 m é justamente a célula do nível ''L20'', portanto terá <math>20*2=40</math> bits de comprimento, concatenados aos 4 bits de L0 e ao prefixo do país. | Cada célula da cobertura ''L0'' já tem seu código de 4 bits (dígito hexadecimal), em seguida qualquer outra célula de qualquer nível hierárquico terá um código ''cbits'' do nível ''L'' com ''L*2'' bits obtidos das coordenadas IJ através da função <code>ints_to_interleavedbits(i,j)</code>. O tamanho de lado da célula do nível ''L'', no caso do Brasil, será <math>2^{20-L}</math> metros. A célula de 1 m é justamente a célula do nível ''L20'', portanto terá <math>20*2=40</math> bits de comprimento, concatenados aos 4 bits de L0 e ao prefixo do país. | ||
A degeneração geométrica, de quadrado para retângulo, é relativa ao segundo argumento de <code>vbit_interleave(x,y)</code>. Como a função é sempre chamada com a mesma ordem dos argumentos, sempre teremos ou só retangulos orizontais (XY) ou só verticais (YX). | |||
=== Algoritmo e funções finais de resolução === | === Algoritmo e funções finais de resolução === | ||
Edição das 06h20min de 20 de junho de 2024
Lib de Grades
Biblioteca para simplificar a geração de grades. Exemplificando implementação dentro do schema grid_br das grades do Brasil. Convensionou-se em GGeohash usar YX no lugar de XY.
A seguir avaliando o uso default de XYref no lugar de XYcenter, para o desenho da célula.
drop FUNCTION if exists grid_br.xyS_collapseTo_ijS(int,int,int);
drop FUNCTION if exists grid_br.xyS_collapseTo_ijS(int[]);
drop FUNCTION if exists grid_br.xylevel_collapseTo_ijS(int,int,int);
drop FUNCTION if exists grid_br.xylevel_collapseTo_ijS(int[]);
CREATE FUNCTION grid_br.xyS_collapseTo_ijS(y int, x int, s int default 1048576) RETURNS int[] AS $f$
SELECT array[ (y-6727000)/s, (x-2715000)/s, s ]
$f$ LANGUAGE SQL IMMUTABLE;
CREATE FUNCTION grid_br.xyS_collapseTo_ijS(yxs int[]) RETURNS int[] AS $wrap$
SELECT grid_br.xyS_collapseTo_ijS(yxs[1], yxs[2], yxs[3])
$wrap$ LANGUAGE SQL IMMUTABLE
;
CREATE FUNCTION grid_br.xylevel_collapseTo_ijS(y int, x int, level int default 0) RETURNS int[] AS $wrap$
SELECT grid_br.xyS_collapseTo_ijS($1, $2, (2^(20-level))::int)
$wrap$ LANGUAGE SQL IMMUTABLE
;
CREATE FUNCTION grid_br.xylevel_collapseTo_ijS(yxl int[]) RETURNS int[] AS $wrap$
SELECT grid_br.xyS_collapseTo_ijS(yxl[1], yxl[2], yxl[3])
$wrap$ LANGUAGE SQL IMMUTABLE
;
-----
drop FUNCTION if exists grid_br.ijS_to_xySref(int,int,int);
drop FUNCTION if exists grid_br.ijS_to_xySref(int[]);
CREATE or replace FUNCTION grid_br.ijS_to_xySref(i int, j int, s int) RETURNS int[] AS $f$
SELECT array[ 6727000 + i*s, 2715000 + j*s, s ]
$f$ LANGUAGE SQL IMMUTABLE;
CREATE FUNCTION grid_br.ijS_to_xySref(yxl int[]) RETURNS int[] AS $wrap$
SELECT grid_br.ijS_to_xySref(yxl[1], yxl[2], yxl[3])
$wrap$ LANGUAGE SQL IMMUTABLE;
Testes para validação, usando funções ou tabelas já homologadas:
select code_b16h, ijs[1] i, ijs[2] j, count(*) n from (
select a.code_b16h, grid_br.xylevel_collapseTo_ijS(st_y(pt.geom)::int,st_x(pt.geom)::int,hlevel::int) as ijs
from grid_br.all_levels a, lateral ST_DumpPoints( ST_GeneratePoints(a.geom, 200) ) as pt
where hlevel=0
) t group by 1,2,3 order by 1;
Comparando com a definição de grade do Brasil, de zero a "e":
grid_l0_cell_ij: 40 41 42 43 30 31 32 33 34 21 22 23 11 12 13 grid_l0_cell_b16: 0 1 2 3 4 5 6 7 8 9 a b c d e code_b16h | i | j | n -----------+---+---+----- 0 | 4 | 0 | 200 1 | 4 | 1 | 200 2 | 4 | 2 | 200 3 | 4 | 3 | 200 4 | 3 | 0 | 200 5 | 3 | 1 | 200 6 | 3 | 2 | 200 7 | 3 | 3 | 200 8 | 3 | 4 | 200 9 | 2 | 1 | 200 a | 2 | 2 | 200 b | 2 | 3 | 200 c | 1 | 1 | 200 d | 1 | 2 | 200 e | 1 | 3 | 200
Conferindo a localização do ponto xyRef:
drop view test_l0_refpoint;
create view test_l0_refpoint as
select *, ST_SetSRID(ST_POINT(yxs[2], yxs[1]), 952019) as geom
from (
select code_b16h, ijs[1] i, ijs[2] j, grid_br.ijS_to_xySref(ijs) as yxs, count(*) n
from (
select a.code_b16h, grid_br.xylevel_collapseTo_ijS(st_y(pt.geom)::int,st_x(pt.geom)::int,hlevel::int) as ijs
from grid_br.all_levels a, lateral ST_DumpPoints( ST_GeneratePoints(a.geom, 20) ) as pt
where hlevel=0 and code_b16h !~ '^f'
) t1 group by 1,2,3,4 order by 1
) t2;
Portanto grid_br.ijS_to_xySref(grid_br.xylevel_collapseTo_ijS(ponto,nivel)) fornece automaticamente a célula L0 de referência onde será aplicado algoritmo GGeohash.
Construtor de geometria
Como vimos a referência é o ponto inferior esquerdo do contorno da célula, basta criar os demais 3 pontos de um quadrado. Na função abaixo a menor célula tem 1 metro pois necessitaria fração para menos que 1 m.
DROP FUNCTION if exists grid_br.xyS_draw_anycell(int,int,int)
;
DROP FUNCTION if exists grid_br.xyS_draw_anycell(int[])
;
CREATE FUNCTION grid_br.xyS_draw_anycell(
y int, -- Yref
x int, -- Xref (XY=canto inferior esquerdo)
s int -- side size
) RETURNS geometry AS $f$
SELECT CASE WHEN s=0 THEN NULL ELSE ST_GeomFromText( format(
'POLYGON((%s %s,%s %s,%s %s,%s %s,%s %s))',
x,y, x+s,y, x+s,y+s, x,y+s, x,y
), 952019) END
$f$ LANGUAGE SQL IMMUTABLE;
COMMENT ON FUNCTION grid_br.xyS_draw_anycell(int,int,int)
IS 'Draws a square-cell using the requested ref-point, with the requested side size.'
;
CREATE FUNCTION grid_br.xyS_draw_anycell(int[]) RETURNS geometry AS $wrap$
SELECT grid_br.xyS_draw_anycell($1[1],$1[2],$1[3])
$wrap$ LANGUAGE SQL IMMUTABLE;
Construtor do identificador cbits
Existem várias formas de obter o identificador baseado em curva de morton. O identificador IJ da célula, como vimos, é fácil de ser obtido no GGeohash.
drop function if exists vbit_to_int(varbit);
drop function if exists vbit_to_int(varbit[]);
drop function if exists vbit_interleave;
drop function if exists ints_to_interleavedbits;
drop function if exists vbit_DeInterleave;
DROP AGGREGATE if exists bitcat_agg(varbit);
CREATE AGGREGATE bitcat_agg(varbit) (
STYPE = varbit
,SFUNC = bitcat -- \do+ "||"
);
CREATE FUNCTION vbit_to_int(x varbit) returns int as $f$
SELECT ( substring(0::bit(32),bit_length(x)+1) || x )::bit(32)::int
$f$ LANGUAGE SQL IMMUTABLE;
CREATE FUNCTION vbit_to_bigint(x varbit) returns int as $f$
SELECT ( substring(0::bit(64),bit_length(x)+1) || x )::bit(64)::bigint
$f$ LANGUAGE SQL IMMUTABLE;
CREATE FUNCTION vbit_to_int(x varbit[]) returns int[] as $f$
SELECT array_agg( vbit_to_int(x_i) ORDER BY i) FROM unnest(x) WITH ORDINALITY t(x_i,i)
$f$ LANGUAGE SQL IMMUTABLE;
CREATE FUNCTION vbit_interleave(x varbit, y varbit) returns varbit as $f$
SELECT bitcat_agg( substring(x,i,1) || substring(y,i,1) )
FROM generate_series(1,bit_length(x)) t(i)
$f$ LANGUAGE SQL IMMUTABLE;
CREATE FUNCTION ints_to_interleavedbits(x int, y int, len int default 32) returns varbit as $f$
SELECT vbit_interleave( substring(x::bit(32),33-len), substring(y::bit(32),33-len) )
$f$ LANGUAGE SQL IMMUTABLE;
-- SELECT ints_to_interleavedbits(47,19,6); -- 100110101111
CREATE FUNCTION vbit_DeInterleave(x varbit) returns varbit[] as $f$
SELECT array[ string_agg(substring(x,i,1)::text,'')::varbit,
string_agg(substring(x,i+1,1)::text,'')::varbit
]
FROM generate_series(1,bit_length(x),2) t(i)
$f$ LANGUAGE SQL IMMUTABLE;
-- SELECT vbit_to_int(vbit_DeInterleave(b'100110101111')); -- 47, 19
CREATE FUNCTION vbit_DeInterleave_to_int(x varbit) returns int[] as $wrap$
SELECT vbit_to_int( vbit_DeInterleave(x) )
$wrap$ LANGUAGE SQL IMMUTABLE;
Os indexadores IJ são a forma mais compacta (normalizada) de representação da grade, e independente da projeção ou métrica usadas. A função ints_to_interleavedbits(i,j) nos fornece exatamente o indexador desejado.
Cada célula da cobertura L0 já tem seu código de 4 bits (dígito hexadecimal), em seguida qualquer outra célula de qualquer nível hierárquico terá um código cbits do nível L com L*2 bits obtidos das coordenadas IJ através da função ints_to_interleavedbits(i,j). O tamanho de lado da célula do nível L, no caso do Brasil, será metros. A célula de 1 m é justamente a célula do nível L20, portanto terá bits de comprimento, concatenados aos 4 bits de L0 e ao prefixo do país.
A degeneração geométrica, de quadrado para retângulo, é relativa ao segundo argumento de vbit_interleave(x,y). Como a função é sempre chamada com a mesma ordem dos argumentos, sempre teremos ou só retangulos orizontais (XY) ou só verticais (YX).
Algoritmo e funções finais de resolução
Algoritmo principal:
- pt na projeção do Brasil, ou seja, com coordenadas planas x e y.
ij0=grid_br.xyS_collapseTo_ijS(x,y); xy0=grid_br.ijS_to_xySref(ij0);cbits0 = grid_br.IJ0_to_L0( ij0, false )cbits = cbits0 || ints_to_interleavedbits(x-x0, y-y0, 20-level_desejado)
É um pouco mais complexo, para contemplar os casos de pontos sobre cobertura fantasma. Existe uma condição de validade e um ajuste do ponto ao nível:
SE lenght(cbits0)>4 e level_desejado<1.5 THEN NULL; ELSE recalcula xy0 dentro da célula.
drop FUNCTION if exists grid_br.xyS_to_cbits(int,int,int,boolean)
;
-- cuidado é YX! ver subtração
CREATE FUNCTION grid_br.xyS_to_cbits(
y int,
x int,
s int,
use_country_prefix boolean false
) RETURNS varbit AS $f$
-- falta conferir se lenght(cbits0)>4 AND s<3
SELECT CASE WHEN use_country_prefix THEN b'0001001100' ELSE b'' END
|| cbits0 || ints_to_interleavedbits(x-xy0[1], y-xy0[2], 20-s)
FROM (
SELECT grid_br.ijS_to_xySref(ij0) as xy0, grid_br.IJ0_to_L0(ij0,false) as cbits0
FROM (SELECT grid_br.xyS_collapseTo_ijS(x,y)) t0(ij0)
) t1
$f$ LANGUAGE SQL IMMUTABLE;
Demais convenções
...
Problemas de exatidão e problemas topológicos, relativos a descontinuidades ou fronteiras, apesar de não parecerem relevantes do ponto de vista prático, precisam ser formalmente tratados para evitar bugs inesperados quando "der azar" de surgir na prática um desses pontos especiais.
Célula dona da borda
As bordas, ou seja, a linha divisória entre duas células, não pode ser uma "linha sem dono", precisa pertencer a uma das células de um dos lados da borda. Por convenção DNGS as linhas afetadas pelo ponto de referência pertencem à célula referida.
Na ilsutração mostramos que dois lados da célula permanecem com o rótulo da referência, e outros dois lados pertencem às respectivas células contíguas. Por exemplo a célula "6" da ilustração possui dois lados próprios, o inferior e o esquerdo, e dois lados impróprios, das células "2" e "7".
No caso especial de fronteira nacional, como as células "0" ou "3", onde não existem "células contíguas", o dono é o país estrangeiro, portanto, a rigor, a linha não existe na grade.
teste
Fonte 80%
select lkjsdklsj+2;
Fonte normal:
select lkjsdklsj+2;