osmc:Metodologia/Algoritmo SQL/Lib: mudanças entre as edições
m (→Fórmulas e definições: link) |
(→Fórmulas e definições: add discretização) |
||
Linha 73: | Linha 73: | ||
A fórmula de ''S'' também nos permite calcular, conforme [[osmc:Metodologia/Algoritmo_SQL/Issues#Issue_05_-_Cálculo_de_dígitos_base32|issue 05]], o número de dígitos base32 para se chegar no metro, a partir do nível ''Lcover'' do município: <math>Ndig1m = 1+\lceil (Lmax-Lcover)/2.5 \rceil</math>. O padrão DNGS não impõe a necessidade de interoperabilidade na bas32, mas se for uma decisão soberana, a escolha de ''Lcover'' nos municípios deve obedecer a essa restrição. | A fórmula de ''S'' também nos permite calcular, conforme [[osmc:Metodologia/Algoritmo_SQL/Issues#Issue_05_-_Cálculo_de_dígitos_base32|issue 05]], o número de dígitos base32 para se chegar no metro, a partir do nível ''Lcover'' do município: <math>Ndig1m = 1+\lceil (Lmax-Lcover)/2.5 \rceil</math>. O padrão DNGS não impõe a necessidade de interoperabilidade na bas32, mas se for uma decisão soberana, a escolha de ''Lcover'' nos municípios deve obedecer a essa restrição. | ||
Fórmulas de discretização: dado um ponto <math>{xy}=pt(x,y)</math> com valores contínuos, obter valores discretos <math>{xyLRef}</math>relativos ao canto inferior da célula de nível ''L'' que contém ''pt''. Essa mesma célula pode ser referenciada pelos índices ''j'' e ''i'' dentro da respectiva célula-mãe de cobertura, ou seja, dentro de ''cover_id'', portanto coordenadas <math>{ijL}=({cover_id},j,i)_L</math>. A discretização segue esse processo: collapso do valor contínuo "XY to jiL", depois a reconstrução como valor discreto, "jiL to xyLRef". | |||
* "XY to jiL": a função ''xy_to_cover()'' faz de forma otimizada a identificação da célula de cobertura, devolvendo as coordenadas <math>xy0_{cover_id}=(x0,y0)</math> da origem XY0 da cobertura. O próximo passo é encontrar a posição dentro da ''cover_id'', ou seja, <math>(x-x0,y-y0)</math>. Por fim, teremos duas situações, a de grade quadrada, com ''L'' inteiro e a de grade retangular, com "''L'' meio": | |||
** No nível ''L'' inteiro temos células de lado <math>S_{L}=2^{Lmax-L}</math>, e podemos quebrar em partes inteiras <math>{jiL}=(j,i)_L=(\frac{x-x0}{S_{L}},\frac{y-y0}{S_{L}})</math> | |||
** No nível ''L'' meio temos retângulos de lados <math>S_{L}</math> no eixo Y e <math>2\cdot S_{L}</math> no eixo X. A quebrar em partes inteiras resultará em <math>{jiL}=(j,i)_L=(\frac{x-x0}{2\cdot S_{L}},\frac{y-y0}{S_{L}})</math> | |||
* "jiL to xyRef": a reconstrução é direta, mas vai apresentar fórmulas diferentes para o quadrado e o retângulo. | |||
** ''L'' inteiro: .. | |||
** ''L'' meio: .. | |||
==Core Lib== | ==Core Lib== |
Edição das 05h32min de 16 de julho de 2024
Biblioteca para simplificar a geração do sistema de grades de um país. Exemplificando implementação dentro do schema grid_br
das grades do Brasil.
A seguir avaliando o uso default de XYref no lugar de XYcenter, para o desenho da célula. Algumas funções de flexibilização precisam ser agrupadas como "helper functions", enquanto as demais como "core functions". Por exemplo as convenções de array de inteiros baseadas em nível (xyL e ijL) são core, enquanto as baseadas em size-side (xyS e ijS) são helper.
Sistemas de coordenadas
![](/w/images/thumb/b/b0/DNGS-coords1.png/420px-DNGS-coords1.png)
O projeto como um todo faz uso de 2 sistemas de coordenadas globais:
- LatLong: coordenadas de latitude e longitude. A rigor coordenadas WGS84 no formato default Geo URI.
- global_cell_id: global mosaic of DNGS cells. O mosaico global é composto de faces quadrilateras, identificadas pelo prefixo de célula:
- face_id: chave primária de 12 bits (SRID, country_id, country_cover_id), com 8 bits para SRID×country_id e 4 bits para country_cover_id.
- cell_id: identificador local de célula, dado por bit string com tamanho variável. Começa pelo vazio, que resulta na própria country_cover_id, como grade nível L0. Depois a cada 2 bits temos um nível hierárquico (L1, L2, ...) até o nível Lmax ou completar 56 bits (limite para hInt). Em caso de fronteira apenas porções do território nacional são válidas, mas essa validade é presumida.
- XY: a cada conjunto de países com mesmo SRID pode-se adotar as mesmas coordenadas planas, de modo que XY é opcional a LatLong como entrada em tal contexto.
Não confundir com o plano canônico XY limitado a cada face.
Na maior parte das aplicações o global_cell_id pode ser simplificado para um identificador local da nação, com face_id=country_cover_id, ou seja, omitindo SRID e country_id quando são constantes. Todavia o ganho com isso, por não alterar o total máximo de 64 bits, não justifica a perda de interoperabilidade, sugere-se manter pelo menos internamente no banco de dados.
Para as aplicações globais foram definidos os métodos seletor de projeção e seletor de jurisdição, que devolvem a face_id e seus parâmetros.
No padrão DNGS polígonos, linhas e multipontos só podem ser definidos no interior de uma nação.
O algoritmo de GGeohash encode da face_id fornece o cell_id dentro da respectiva face. Esse tipo de algoritmo pode fazer uso outros dois sistemas de coordenadas:
COORDINATES AND STANDARD PARAMETERS: canonical jiL: [cover_id, j,i, L] complete jiL: [cover_id, j,i, L, s, X0, Y0, is_half] canonical xyL: [x,y,L] canonical xyLRef: [x,y,L] complete xyLRef: [x,y,L, cover_id, j,i, s, X0, Y0, is_half] canonical xyS: [x,y,s,is_half] complete xyS: [x,y,s,is_half, cover_id, j,i, L, X0, Y0] PS: xyL is request-only, countinuous values. xyLRef is discrete. xyS is request or generic (high performance) is_half is true for degenerated grid by 0.5*s, NULL for stdandard grid , false for degenerated by 2*s. Use false for standard degenerated.
- JI0: cada face_id pode ser representado por índices j0 e i0 da grade L0, que por sua vez podem levar à origem de cada célula da cobertura, nas coordenadas da projeção SRID. Na prática o que temos também é a origem em termos de face_id.
- XY da face_id no nível L: o sistema de coordenadas planas SRID é reescrito para cada face_id, para que inicie em (0,0). Segue a "convenção cartesiana escolar", de X horizontal e Y vertical.
Quando discretizado é denominado XYLref, complementar de JIL. - JIL, índices J e I da célula-filha da face_id no nível L: no interior da face, cada célula de nível L pode ser localizada por seus índices j e i. O encode/decode de cell_id é realizado primariamente com essas coordenadas. Geometricamente cada jiL pode ser associado univocamente ao retângulo da sua célula, ou a um único XYref (canto inferior esquerdo) ou ao seu XYcenter.
O sistema JI difere da tradição matemática de indexação de matrizes, com elementos , índice i variando nas linhas e j variando nas colunas. A origem tradicional é no canto superior esquerdo, enquanto a origem XY do eixo cartesiano é no canto inferior.
Para equiparar a orientação de I×J a X×Y é necessário aplicar a rotação de 90 graus no sentido anti-horário sobre I×J, e trocar a ordem das coordenadas, ou seja, convencionar J×I.
Por fim, na convenção cartesiana XY a origem é (0,0) e a equação pode ser simplificada para quando mudamos o domínio de I, do intervalo inteiro [1, m] para [0, m-1]; e o domínio de J para [0, n-1]. Com esse ajuste final temos a convenção DNGS para índices JI. O resultado desta convenção JI pode ser apreciado pela representação da Curva-Z:
Fórmulas e definições
Por imposição do padrão DNGS temos:
- As células são quadrilateras, idealmente quadradas.
- A projeção é de igual-área, ou bastante próxima disso.
- O sistema de grades de nível L é um mosaico quadrilátero, que tem o metro como maior resolução, no nível Lmax: .
- A grade de nível zero, L=0, é um conjunto de até 16 quadrados de área , que cobre todo o domínio territorial nacional. Portanto quadrados de lado , ou seja, obriga-se o uso de uma potência de 2 na cobertura nacional.
- A resolução da grade de nível L é , interpretada como tamanho do lado da célula quadrada, de área .
No caso do Brasil Lmax=20, no caso de Camarões Lmax=18. Essa escolha em geral pode ser aproximada pela área nacional, AN, ou seja, , portanto .
A fórmula de funciona também para níveis-meio, , por exemplo L=1,5. São níveis com grade de geometria degenerada, retangular. Duas importantes suposições levam à demonstração de valide:
- A interpretção de S como "tamanho do lado genérico S de um retângulo" e seu valor como raiz quadrada da área do retângulo, . O entendimento da construção geométrica da célula de nível-meio como união de 2 células do próximo nível inteiro, Lprox, . Portanto, fazendo temos:
.
A cada face_id
, ou seja, a cada célula da cobertura L0 do país, podemos calcular o número de células-filhas no nível L:
É invariante por face e por país, e válido para "níveis meio".
O padrão DNGS também prevê a cobertura municipal e regras para a construção da cobertura e do código logístico:
O número de células por município de área M, no nível L, pode ser obtido por . A cobertura uniforme municipal de c celulas do nivel K tem área .
Portanto número de células na área C, no nível L>K, será
O valor de C sempre será maior do que M, de modo que o fator C/M será sempre maior que 1, ou . Dessa desigualdade de áreas, substituindo pelos valores correspondentes de número de células, podemos provar que .
A fórmula de S também nos permite calcular, conforme issue 05, o número de dígitos base32 para se chegar no metro, a partir do nível Lcover do município: . O padrão DNGS não impõe a necessidade de interoperabilidade na bas32, mas se for uma decisão soberana, a escolha de Lcover nos municípios deve obedecer a essa restrição.
Fórmulas de discretização: dado um ponto com valores contínuos, obter valores discretos relativos ao canto inferior da célula de nível L que contém pt. Essa mesma célula pode ser referenciada pelos índices j e i dentro da respectiva célula-mãe de cobertura, ou seja, dentro de cover_id, portanto coordenadas . A discretização segue esse processo: collapso do valor contínuo "XY to jiL", depois a reconstrução como valor discreto, "jiL to xyLRef".
- "XY to jiL": a função xy_to_cover() faz de forma otimizada a identificação da célula de cobertura, devolvendo as coordenadas da origem XY0 da cobertura. O próximo passo é encontrar a posição dentro da cover_id, ou seja, . Por fim, teremos duas situações, a de grade quadrada, com L inteiro e a de grade retangular, com "L meio":
- No nível L inteiro temos células de lado , e podemos quebrar em partes inteiras
- No nível L meio temos retângulos de lados no eixo Y e no eixo X. A quebrar em partes inteiras resultará em
- "jiL to xyRef": a reconstrução é direta, mas vai apresentar fórmulas diferentes para o quadrado e o retângulo.
- L inteiro: ..
- L meio: ..
Core Lib
Tanto o padrão DGGS como o DNGS são relativos a "Discrete Grid Systems", ou seja, um conjunto hierárquico de grades, cada qual no seu nível L, e assumindo apenas valores discretizados de posição. As "funções core" do padrão DNGS garantem a ida e volta entre os diversos sistemas de coordenadas necessários à manipulação dos dados, antes e depois da discretização.
Discretizadores: transformam as coordenadas contínuas XY, de posição no plano projetado, em coordenadas discretas IJ, de localização na grade hierárquica. Quem amarra a posição hierárquica — grid hierarchical level do sistema de grades — com localização XY é o tamanho de lado S (side size) da célula de nível L.
-- Core functions:
drop FUNCTION if exists grid_br.xyS_collapseTo_ijS(int,int,int,boolean);
drop FUNCTION if exists grid_br.xyL_collapseTo_ijL(int[]);
CREATE FUNCTION grid_br.xyS_collapseTo_ijS(y int, x int, s int default 1048576, is_half boolean default false) RETURNS int[] AS $f$
SELECT array[ (y-6727000)/CASE WHEN is_half AND s>=2 THEN s/2 ELSE s END, (x-2715000)/s, s ]
$f$ language SQL IMMUTABLE;
CREATE FUNCTION grid_br.xyL_collapseTo_ijL(xyL int[]) RETURNS int[] AS $f$
SELECT array[ $ijS[1], $ijS[2], intlevel ]
FROM (
SELECT grid_br.xyS_collapseTo_ijS( $1[1], $1[2], (2^(20-intlevel/10.0))::int, (intlevel%2)=1 ) AS ijS
FROM ( SELECT COALESCE($1[3],0) ) t1(intlevel)
) t2
$f$ language SQL IMMUTABLE;
Construção da 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)
;
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.'
;
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 (XY).
Otimização do cálculo de cobertura
Ver Tratamento das configurações onde já foi discutido e solucionado o tema. Aqui retomando para questões de otimização.
A definição nacional, no caso do Brasil é
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 fT fP fN
A entrada é o ponto pt=(x,y). O primeiro passo é obter as coordenadas ij0 de pt na cobertura. Temos dois possíveis algoritmos:
- A função discretizadora
ijL0=grid_br.xyS_collapseTo_ijL(x,y)
, onde o tamanho de célula s0 é o default.
- Uma árvore de decisão com 4 nívieis de IFs (análogo a índice BBOX) pode nos fornecer rapidamente o valor ijL0.
O que é mais rápido, depende de benchmark na linguagem adotada (C ou SQL). Por hora ficamos com o primeiro algoritmo.
Como ijL0 é um caso muito especial de ij podemos obter o seu valor num formato mais conveniente. Exemplos:
ij0 := i::text||j::text
poderia usar um objeto JSONb como array associativa, mas não é tão rápido quanto array SQL. Seria algo hardcoded numa função,('{"40":[123,456,"0"],"41":[0,1,"1"],"13":[22,33,"e"]}'::jsonb)->ij0
.
i*10+j
resulta em array de 44 posições. Seriam arrays leves; e delas, a partir de ij0, teremos:- coordenadas iniciais x0 e y0, por arrays
x0_from_ij0
ex0_from_ij0
; - índice cbits0, pela array
cbits0_from_ij0
.
- coordenadas iniciais x0 e y0, por arrays
Por fim as arrays podem ser "hardcoded" nas funções usuárias, otimizando a obtenção dos valores desejado.
Lembrar que para obter a array the uma multiarray não é simples em SQL, exemplo de obtenção do indice i: SELECT array(select unnest( ('{{1,2,3},{4,5,6},{7,8,9}}'::int[][])[i:i] ));
Algoritmo e funções finais de resolução
Algoritmo principal, tendo como entradas: pt e nível L. Ponto pt em coordenadas planas, portanto x e y.
- Célula da cobertura nacional:
ij0=grid_br.xyS_collapseTo_ijS(x,y);
com s0 é default, ou função otimizada.x0=x0_from_ij0[ij0]; y0=y0_from_ij0[ij0];
ou função de ij0 retorando xy0.cbits0 = grid_br.IJ0_to_vbitL0( ij0, false )
- Código cbits e geometria da célula do nível L:
ijL=grid_br.xyL_collapseTo_ijL(x-xy0[1], y-xy0[2], L);
cbits = cbits0 || ints_to_interleavedbits(ijL)
grid_br.xyS_draw_anycell( grid_br.ijL_to_xyL(ijL) )
Com um passo a mais 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 especial.
drop FUNCTION if exists grid_br.xyS_to_cbits(int,int,int,boolean)
;
-- revisar se é XY!
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( grid_br.xyS_collapseTo_ijL(x-xy0[1], y-xy0[2], 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
- Pendente confirmação e inclusão de ASSERT nas founções. (ou correção da regra para que a função estabeleça a convenção)
As bordas são conjuntos de pontos, formados pela linha divisória entre 2 células e pelo ponto de fronteira entre 4 células. Apesar de serem infinitesimais, as bordas não podem ser "pontos sem dono". Cada ponto de borda precisa pertencer a uma das células, de um dos lados da borda. São convenções técnicas que definem esse pertencimento:
- Por convenção DNGS as linhas afetadas pelo ponto de referência pertencem à célula referida. Todos os pontos de fronteira entre 4 células, portanto, possuem um dono.
Exemplo: o ponto XY=(0.0,0.0) pertence à célula IJ=(0,0). A função ijS_to_xySref() fui concebida de modo a garantir essa convenção (e vice-versa). - Todos os demais pontos, sobre as linhas de borda entre 2 células, pertencem à célula de cantos inferior e esquerdo sobre as respectivas linhas. São as linhas dos eixos locais XY, com origem no canto inferior esquerdo, dada pelo ponto de referência.
Na ilustração abaixo 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" do Brasil, onde não existem "células contíguas", o dono é o país estrangeiro, portanto, a rigor, a linha inteira (incluindo o canto do quadrado) não existe na grade.
Geração do sistema de grades
Algoritmo genérico, para gerar grades nacionais ou municipais. Existem dois algoritmos mais gerais de geração: "scan IJ" e scan NatCodes".
A função grid_br.parents_to_children_cuting
implementa de maneira simples o algoritmo "scan NatCodes", que pode ser expresso como recorrência:
- Gera células-filhas da cobertura.
- Remove filhas que ficaram fora da interseção da geometria do país (ou município).
- Assume filhas como cobetura e volta ao 1.
-- primeiro esboço da função
CREATE or replace FUNCTION grid_br.parents_to_children_cuting(
p_intlevel int, -- Target level, where level=intlevel/10.0. Use parent as levelZero? no: grid_br level.
p_parent_list varbit[], -- Natcod parents. A list of distinct grid-cell IDs.
p_cutBy_geom geometry -- when noT null
-- decidir se faz cut de fato ou se apenas confere interseção
) RETURNS setof table_name
language SQL VOLATILE
AS $f$
with maioria da cobertura me dá o nivel corrente
WITH generated AS (
SELECT *,
CASE WHEN p_cutBy_geom IS NULL THEN true ELSE ST_Contains(p_cutBy_geom,geom) END
AS is_contained
FROM (
SELECT cbits, null as geom -- se foi cortado e st_area() é a area do nível, entao está contido.
FROM (select cbits FROM unnest(p_parent_list)) t1
WHERE p_parent_list IS NOT NULL -- guard, for stop when null
) t2
),
cut AS (
SELECT * FROM (
SELECT cbits, is_contained,
CASE
WHEN is_contained OR p_cutBy_geom IS NULL THEN geom
ELSE ST_Intersection(p_cutBy_geom,geom)
END AS geom
FROM generated
) t3
WHERE geom IS NOT NULL
)
SELECT * FROM cut
UNION ALL
SELECT * FROM grid_br.parents_to_children_cuting(
$1,
CASE WHEN bit_length(cbits) p_intlevel+1>20 THEN NULL ELSE (SELECT array_agg(cbits) FROM cut) END,
$3
) t4
$f$;
Casos de uso da grade completa e parcial
A persistência do sistema de grades tem um custo muito alto, mas em geral, mesmo sem geometria, a lista de identificadores de célula terá um papel importante como chave-primária de conjuntos de atributos, no que consiste o sistema de informação. Existem duas estratégias para se persistir informações:
- geo-field: atributos armazenados em todas as células, na grade inteira, eventualmente limitada pela resolução. Geocampos de alta resolução podem ser armazenados de forma separada, como raster ou geomvals.
- geo-objects: pontos, linhas ou areas, cada tipo tem seu tratamento dimensional distinto. Uma célula de área pode ser usada como cobertura, a mesma célula se considerada ponto terá a sua área considerada apenas como incerteza em torno do centro.
O sistema completo de grades é mais importante no caso de geo-fields. Os primeiros níveis podem ser geo-fields orientadores, para descobrir onde se encontram geo-objects de interesse. Municípios, devido ao interesse por gestão polcal podem formarmar seus geofields locais, um só banco de dados para a getão municipal (não o país inteiro), viabilizando o armazenamento, que requer capacidade exponencial de disco em função do nível.
Helper lib
Funções complementares às "core functions", a maioria implementações "wrap" para simplesmente compatibilizar tipos e estruturas. Ainda assim a ortogonalidade é requerida apenas nos tipos principais (ex. baseados no intLevel), não precisam garantir a conversão direta de tipos secundários (ex. cell side size).
drop FUNCTION if exists grid_br.xyS_collapseTo_ijS(int[],boolean);
drop FUNCTION if exists grid_br.xylevel_collapseTo_ijS(int,int,int);
CREATE FUNCTION grid_br.xyL_collapseTo_ijS(y int, x int, intlevel int default 0) RETURNS int[] AS $wrap$
SELECT grid_br.xyS_collapseTo_ijS($1, $2, (2^(20-intlevel/10.0))::int, (intlevel%2)=1)
$wrap$ LANGUAGE SQL IMMUTABLE
;
CREATE FUNCTION grid_br.xyL_collapseTo_ijS(xyl int[]) RETURNS int[] AS $wrap$
SELECT grid_br.xyS_collapseTo_ijS(xyl[1], xyl[2], xyl[3])
$wrap$ LANGUAGE SQL IMMUTABLE
;
-----
drop FUNCTION if exists grid_br.ijS_to_xySref(int,int,int,boolean);
drop FUNCTION if exists grid_br.ijS_to_xySref(int[],boolean);
CREATE or replace FUNCTION grid_br.ijS_to_xySref(i int, j int, s int, is_half boolean default false) RETURNS int[] AS $f$
SELECT array[ 6727000 + i*CASE WHEN is_half AND s>=2 THEN s/2 ELSE s END, 2715000 + j*s, s ]
$f$ LANGUAGE SQL IMMUTABLE;
CREATE FUNCTION grid_br.ijS_to_xySref(xyl int[], is_half boolean default false) RETURNS int[] AS $wrap$
SELECT grid_br.ijS_to_xySref(xyl[1], xyl[2], xyl[3],is_half)
$wrap$ LANGUAGE SQL IMMUTABLE;
Testes da lib
Primeiro demonstrações (de que as funções funcionam), testes para validar hipóteses e requisitos básicos; depois testes de regressão.
Validação primária
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;
Grade da Colômbia, 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(xys[2], xys[1]), 952019) as geom
from (
select code_b16h, ijs[1] i, ijs[2] j, grid_br.ijS_to_xySref(ijs) as xys, 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.