osmc:Metodologia/Algoritmo SQL/Lib

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

Sistemas de coordenadas previstos pelo padrão DNGS: LatLong, XY e IJ. A superfície S² é o globo terrestre, R² o "plano LatLong WGS84" e XY o plano projetado Albers. A passagem para IJ é a discretização. Um mesmo angulo sólido α em locais diferentes do globo terá áreas diferentes em R², mas mesma área em XY.

O projeto como um todo faz uso de 3 sistemas de coordenadas globais:

• LatLong: coordenadas de latitude e longitude. A rigor coordenadas WGS84 no formato default Geo URI.

• XY contínuo e discreto: todo país adota uma projeção SRID, de modo que LatLong será transformado em XY, ainda sem discretização. Depois de discretizado é denominado XYLref, pois a discretização depende do nível L e do ponto "XYL de referência" (canto inferior esquerdo da célula).
  PS: 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 a zona plana XYL limitada a cada face_id.

• 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, com 8 bits para SRID×country_id distintos e 4 bits para country_cover_id (identificador de "face plana" na forma de quadrante da cobertura nacional).
  • cell_id: identificador local de célula, dado por bit string com tamanho variável. Dentro da face cada coordenada discreta jiL é também única e definida pelo algoritmo de encode/decode.
    Detalhes. A representação bit string 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.

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 $(x0_{j0},y0_{i0})$ de cada célula da cobertura, nas coordenadas da projeção SRID. Na prática o que temos também é $(x0_{face\_id},y0_{face\_id})$ a origem em termos de face_id.
• XYLref 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 XY, com X horizontal e Y vertical. Quando discretizado é denominado XYLref, complementar de JIL.
  O ponto central da célula, XYLcnt, é usado apenas em aplicações externas.
• 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 $a_{ij}$, índice i variando nas linhas e j variando nas colunas. A origem tradicional $a_{11}$ é 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 $X_{i}=(i-1)\cdot s$ pode ser simplificada para $X_{i}=i\cdot s$ 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.

Alguns algoritmos, como a Curva-Z ilustrada abaixo, nascem na matriz com índices invertidos, JI, e só posteriormente, na construção geométrica, adquirem o aspecto rotacionado. Para se manter o "algoritmo da curva-Z" na grade retangular deveria haver alternância entre JI e IJ. Na prática preferimos manter JI na grade retangular, adotando um algoritmo de curva modificado, denominado "curva Z degenerada".

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: $S_{Lmax}=1$.
• A grade de nível zero, L=0, é um conjunto de até 16 quadrados de área $A_{0}=2^{2\cdot Lmax}$, que cobre todo o domínio territorial nacional. Portanto quadrados de lado $S_{0}=2^{Lmax}$, ou seja, obriga-se o uso de uma potência de 2 na cobertura nacional.
• A resolução da grade de nível L é $S_{L}=2^{Lmax-L}$, interpretada como tamanho do lado da célula quadrada, de área $A_{L}={S_{L}}^{2}=2^{2\cdot (Lmax-L)}$.

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, $16\cdot 2^{2\cdot Lmax}\gtrapprox AN$, portanto $\lfloor 0,5\cdot log_{2}(AN)-2\rfloor \gtrapprox Lmax\gtrapprox \lceil 0,5\cdot log_{2}(AN)-2\rceil$. Uma estimativa mais efetiva pode ser conseguida substituindo-se AN pela área da BBOX quadrada que contém AN. Um território alongado, como o do Chile, vai ocupar só uma fração da área da sua BBOX.

A fórmula de $S_{L}$ funciona também para níveis-meio, $Lhalf\in \{\forall x|x=\lceil x\rceil -0,5\}$, 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", com valor obtido pela raiz quadrada da área do retângulo, $S_{Lhalf}={\sqrt {A_{Lhalf}}}$; e a interpretação 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, $A_{Lhalf}=2\cdot {A_{Lprox}}$. Portanto, fazendo $Lprox={\lceil Lhalf\rceil }$ temos:
  $S_{Lhalf}={\sqrt {2\cdot {A_{Lprox}}}}={\sqrt {2}}\cdot {\sqrt {{S_{\lceil Lhalf\rceil }}^{2}}}=2^{Lmax-\lceil Lhalf\rceil +0,5}$.

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:
  $nCells_{L}=(S_{0}/S_{L})^{2}=(2^{Lmax-Lmax+L})^{2}=2^{2L}$
É 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 $nCellsM_{L}=M/{S_{L}}^{2}=M\cdot 2^{2\cdot (L-Lmax)}$. A cobertura uniforme municipal de c celulas do nivel K tem área $C=c\cdot {S_{K}}^{2}=c\cdot 2^{2\cdot (Lmax-K)}$.
Portanto número de células na área C, no nível L>K, será
  $nCellsC_{L}=C/{S_{L}}^{2}=c\cdot 2^{2\cdot (Lmax-K)}\cdot 2^{2\cdot (L-Lmax)}=c\cdot 2^{2\cdot (Lmax-K)+2\cdot (L-Lmax)}=c\cdot 2^{2\cdot (L-K)}$
O valor de C sempre será maior do que M, de modo que o fator C/M será sempre maior que 1, ou $C/M\geq 1$. Dessa desigualdade de áreas, substituindo pelos valores correspondentes de número de células, podemos provar que $nCellsC_{L}\geq nCellsM_{L}$.

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: $Ndig1m=1+\lceil (Lmax-Lcover)/2.5\rceil$. 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 ${xy}=pt(x,y)$ com valores contínuos, obter valores discretos ${xyLRef}$ 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 ${jiL}=({\text{id}},j,i)_{L}$. A discretização segue esse processo: colapso 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 cover_id, devolvendo as coordenadas $xy0_{\text{id}}=(x0,y0)_{\text{id}}$ da origem XY0 da célula. O próximo passo é encontrar a posição dentro da cover_id, ou seja, $(x-x0,y-y0)$. 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 $S_{L}=2^{Lmax-L}$, e podemos quebrar em partes inteiras ${jiL}=({id},j,i)_{L}=({id},\lfloor {\frac {x-x0}{S_{L}}}\rfloor ~,~\lfloor {\frac {y-y0}{S_{L}}}\rfloor )_{L}$
  • No nível L meio (Lhalf) temos retângulos de lados $S_{Lprox}$ no eixo Y e lados $2\cdot S_{Lprox}$ no eixo X, onde ${Lprox}=\lceil Lhalf\rceil$. A quebra em partes inteiras resultará em ${jiL}=({\text{id}},j,i)_{Lhalf}=({id},\lfloor {\frac {x-x0}{2\cdot S_{Lprox}}}\rfloor ~,~\lfloor {\frac {y-y0}{S_{Lprox}}}\rfloor )_{Lhalf}$
• "jiL to xyLRef": a reconstrução é direta, mas vai apresentar fórmulas diferentes para o quadrado e o retângulo.
  • L inteiro: ${xyLRef}=(j\cdot S_{L}+x0~,~i\cdot S_{L}+y0)$
  • L meio: ${xyLRef}=(j\cdot 2\cdot S_{Lprox}+x0~,~i\cdot S_{Lprox}+y0)$

Devido à diferença de cálculo entre L inteiro e meio, e ao fato de $2^{x}$ poder ser otimizado por rotação de bits, 2<<(x-1), as funções otimizadas adotam flag is_half conjugado a Lprox inteiro no lugar de L flutuante. Com isso garante-se também o uso de aritmética inteira em todas as etapas do processamento de grade, resultando em maior confiabilidade e reversibilidade, além de performance.

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

A curva de Morton emerge do entrelaçamento dos bits. Exemplo:
ints_to_interleavedbits(47,19,6).

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.

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_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 AGGREGATE bitcat_agg(varbit) ( STYPE = varbit, SFUNC = bitcat ); -- op. \do+ "||"
----

CREATE FUNCTION vbit_interleave(x varbit, y varbit) returns varbit as $f$
  SELECT CASE WHEN x=b'' OR y=b'' THEN x||y ELSE (
    SELECT  bitcat_agg( substring(x,i,1) || substring(y,i,1) )
    FROM generate_series(1,bit_length(x)) t(i)
  ) END
$f$ LANGUAGE SQL IMMUTABLE;

CREATE FUNCTION vbit_DeInterleave(x varbit) returns varbit[] as $f$
  SELECT  array[  bitcat_agg(substring(x,i,1)),
                 bitcat_agg(substring(x,i+1,1))
          ]
  FROM generate_series(1,bit_length(x),2) t(i)
$f$ LANGUAGE SQL IMMUTABLE;

CREATE FUNCTION vbit_DeInterleave_to_ints(x varbit) returns int[] as $wrap$
  SELECT vbit_to_int( vbit_DeInterleave(x) )
$wrap$ LANGUAGE SQL IMMUTABLE;
-- SELECT vbit_DeInterleave_to_ints(b'100110101111'); -- 47, 19

CREATE or replace FUNCTION ints_to_interleavedbits(
  x int, y int, len int default 32, is_half boolean default null
) returns varbit
as $f$
  SELECT vbit_interleave( substring(x::bit(32),33-len), substring( y::bit(32), 33-len+CASE WHEN is_half THEN 1 ELSE 0 END ) )
$f$ LANGUAGE SQL IMMUTABLE;
-- SELECT ints_to_interleavedbits(47,19,6);      -- 100110101111
-- SELECT ints_to_interleavedbits(47,19,6,true); --  11001011111

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á $2^{20-L}$ metros. A célula de 1 m é justamente a célula do nível L20, portanto terá $20*2=40$ 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).

A inversa, vbit_DeInterleave pode requerer a concatenação do bit zero na frente para manter os pares. A inclusão de zero no inicio ou no fim determina qual a dimensão (eixo X ou eixo Y) será distorcida. Essa operação é chamada de "refinamento anisotrópico"ref.

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 e x0_from_ij0;
  • índice cbits0, pela array cbits0_from_ij0.

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.

1. 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 )
2. 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

Grade e seus rótulos de célula, em púrpura claro. Os pontos de referência XYref foram destacados com estrelas vermelhas, e o rótulo da "célula dona do ponto" em preto. Linhas de borda, rotuladas também em preto pelas células com cantos inferior e esquerdo sobre a 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:

1. 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).
2. 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:

1. Gera células-filhas da cobertura.
2. Remove filhas que ficaram fora da interseção da geometria do país (ou município).
3. 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;

---

Grades municipais logísticas

Convencionamos que a grade logística padrão DNGS, com dígitos base 32, precisa apenas 6 dígitos para "chegar na porta de casa", o que em municípios pequenos é a localização com incerteza de 1m², e nos maiores de 32 m². Cada dígito da base 32 equivale a um acréscimo de 2.5 níveis. Como o metro é Lmax, o próximo seria Lmax-2.5 depois Lmax-5.0, Lmax-7.5, ... e assim por diante.

No Brasil, partindo da grade de 1m² (L20) temos mais 5 níveis: L20, L17.5, L15, L12.5, L10, L7.5. Quando a área é maior podemos adotar L5 na cobertura, mas para manter o total de 6 níveis é preciso interromper antes do 1m²: L17.5, L15, L12.5, L10, L7.5, L5. Ou seja, as coberturas interoperáveis são L5 e L7.5.

A estimativa do nível da célula de cobertura do município pode ser obtida da área do município, Lmax-log2(area)/2. Disso podemos traçar um perfil do país para entender quantos seriam interoperáveis. Aplicando ao caso do Brasil:

WITH s AS (
  SELECT *, ceil( (2*20-log(area)/log(2))/2.0 ) as level
  FROM (
     SELECT *, st_area(geom,true) as area, ST_Perimeter(geom,true) as perim
     FROM optim.jurisdiction_geom where isolabel_ext like 'BR-%-%'
  ) t0
), t AS (
  SELECT COUNT(*) ntot, sum(area) areatot, sum(perim) prmtot FROM s 
)
 SELECT CASE WHEN level >=7 THEN '7.5' when level <=4 THEN '~'||level::text ELSE '5' END as Level,
       count(*) n_municipios,
       round(100.0*count(*)::float/t.ntot)::text||'%' as n_perc,
       round(100.0*sum(area)/areatot)::text||'%' as area_perc,
       round(100.0*sum(perim)/prmtot)::text||'%' as perim_perc
 FROM s, t
 GROUP BY 1, areatot, ntot, prmtot
 ORDER BY 1
;

 level | n_municipios | n_perc | area_perc | perim_perc 
-------+--------------+--------+-----------+------------
 ~2    |            9 | 0%     | 10%       | 2%
 ~3    |           57 | 1%     | 21%       | 7%
 ~4    |          283 | 5%     | 26%       | 17%
 5     |         3300 | 59%    | 39%       | 60%
 7.5   |         1921 | 34%    | 4%        | 13%

• Portanto 93% dentro do padrão interoperável (forçados aos níveis L5 ou L7.5).
• Apenas ~350 municípios gigantes demandando melhor acomodação da sua grade de cobertura, sem "forçar o nível".
  • Por serem gigantes, existe maior probabilidade deles sofrem desmembramento (representam mais de 50% da área nacional), do que incorporação.
• A chance de um Estado (ou o Senado) determinar uma grande alteração de limites geográficos é pequena, e ainda assim, a chance de ser entre municípios interoperáveis é de n_perc ~90%.
  • A proposta mais promissora não progrediu, e não existe previsão de outra. Consistia em agregar municípios pequenos (boa parte dos n_perc 34% com 4% da área) ao vizinho mais rico.
• As alterações menores (ex. por melhora de precisão no IBGE ou OSM), medidas por taxa perimetral alterada anualmente, são naturalmente absorvidas pela cobertura vigente, num pior caso por inclusão de cobertura de reserva... Só com muito azar a mudança calha de impactar dois municípios vizinhos: novamente haverá 90% de chances de serem interoperáveis.

Para outros países basta trocar a condição WHERE sobre isolabel_ext e o valor Lmax de 20 para o do país. Camarões é selecionado pela condição isolabel_ext like 'CM-%-%' e seu Lmax é 18, portanto ceil( (2*18.0-log(area)/log(2))/2.0 ). Ajustar também o CASE dos níveis forçados. Partindo da grade de 1m² (L18) temos mais 5 níveis: L18, L15.5, L13, L10.5, L8, L5.5.

Quando a área é maior podemos adotar L3 na cobertura, mas para manter o total de 6 níveis é preciso interromper antes do 1m²: L15.5, L13, L10.5, L8, L5.5, L3. Ou seja, as coberturas interoperáveis são L3 e L5.5.

 level | n_municipios | n_perc | area_perc | perim_perc 
-------+--------------+--------+-----------+------------
 ~2    |           24 | 7%     | 41%       | 21%
 3     |          236 | 66%    | 56%       | 69%
 5.5   |          100 | 28%    | 3%        | 11%

A Colômbia (CO), com Lmax=19 resulta em: ...

Sensibilidade à indexação

Municípios podem ser indexados de diferentes formas, igualmente eficientes, e/ou com critérios de eficiência diferentes. Dois algoritmos podem diferenciar por exemplo com relação à decisão de "poeira", quanto ao tamanho da área que será tomada como insignificantemente para merecer uma cobertura sobre ela. Também com relação à noção de "resiliência da cobertura adotada", forçando um reserva maior ou menor (e portanto afetando o número ou liberdade para uso de overlays urbanos).

Quanto ao uso de overlays e reservas, esse é o critério de ordem: na primeira versão vem primeiro cobertura-base depois cobertua overlay, ordenadas por varbit nos respectivos blocos. Depois de oficializado, vale a ordem estrita da array. As eventuais substituições quebrando a rodem varbit, as inclusões de overlay e reserva sempre no final, ordenadas dentro do possível. O controle de versões e transparência pública em git são também requisitos.

Se mudamos por exemplo a ordem das coberturas, os geocódigos mudam, tornando as duas versões de geocódigo totalmente incompatíveis e com isso gerando grande impacto. Se mudamos o critério de "poeira", eliminando uma célula, a ordenação também será afetada.

Devido a esse risco de variabilidade é necessário impor um padrão nacional. Bons critérios para a escolha desse padrão são a simplicidade e eficiência algorítmicas. A garantia de resiliência também: deixar que a "quase-poeira" seja coberta por overlays ao invés de cobertura-base.

Ver também:

• osmc:Convenções/Coberturas municipais
• osmc:Convenções/Coberturas municipais/Algoritmos#Cobertura-base

Células-folha virtuais

Devido ao crescimento exponencial e à opção de representação das células menores como pontos, podemos dispensar o "cache" dessas células. Cada célula-mãe geraria seus 32 pontos internos e com identificadores consistentes. Uma array de "shift from XYref" da célula-mãe permite calcular rapidamente a posição de cada dígito suplementar da base32 (child_suffix_int).

CREATE VIEW prod.cm_logistic_vw_sheets0test AS
select *, gid::bit(64) as gid_vbit, case when l=180 then 1 else 32 end as s
from (
  select (l.gid<<5)+child_suffix_int as gid,
    l.jurisd_local_id, l.gid as parent_gid, l.l +25 as l,
    code_b32nvu||('{0,1,2,3,4,5,6,7,8,9,B,C,D,F,G,H,J,K,L,M,N,P,Q,R,S,T,U,V,W,X,Y,Z}'::text[])[child_suffix_int+1] as code_b32nvu, geom as orig_geom
  from (select * from prod.cm_logistic where l=130 and code_b32nvu like '9%' limit 2) l, generate_series(0,31) t0(child_suffix_int)
) t1
limit 100;

-- Fornece array da conversão de L13 para XY0ref de L15.5: 
select array_agg(xdiff) xref, array_agg(ydiff) yref
from (
 select s.gid, c.code_b32nvu, child_suffix_int,
       st_xmin(c.geom) - st_xmin(s.orig_geom) as xdiff,
       st_ymin(c.geom) - st_ymin(s.orig_geom) as ydiff
 from prod.cm_logistic_vw_sheets0test s inner join prod.cm_logistic c on c.gid=s.gid
 order by c.gid limit 32
) t;

-- Fornece array da conversão de L13 para XYcenter de L15.5: 
select array_agg(xdiff) xcntr, array_agg(ydiff) ycntr
from (
 select s.gid, c.code_b32nvu, child_suffix_int,
       st_xmin(st_centroid(c.geom)) - st_xmin(s.orig_geom) as xdiff,
       st_ymin(st_centroid(c.geom)) - st_ymin(s.orig_geom) as ydiff
 from prod.cm_logistic_vw_sheets0test s inner join prod.cm_logistic c on c.gid=s.gid
 order by c.gid limit 32
) t;
-- digit_k=('{0,1,2,3,4,5,6,7,8,9,B,C,D,F,G,H,J,K,L,M,N,P,Q,R,S,T,U,V,W,X,Y,Z}'::text[])[k]
-- X_k=st_xmin(s.orig_geom)+Xcntr=('{4,4,12,12,4,4,12,12,20,20,28,28,20,20,28,28,4,4,12,12,4,4,12,12,20,20,28,28,20,20,28,28}'::int[])[k]
-- Y_k=st_ymin(s.orig_geom)+Ycntr=('{2,6,2,6,10,14,10,14,2,6,2,6,10,14,10,14,18,22,18,22,26,30,26,30,18,22,18,22,26,30,26,30}'::int[])[k]

Usando o resultado para uma view construtora de L15.5 a partir de L13, com pontos centrais de retângulos; depois repetindo todo o processo a partir de L10.5 e produzindo L13, com pontos centrais de quadrados.

CREATE VIEW cm_logistic_vw_pts_l130to155 as
 -- prod.cm_logistic_vw_pts_l130to155 AS
 select gid, jurisd_local_id, l, code_b32nvu, ST_POINT(
     x0 + ('{4,4,12,12,4,4,12,12,20,20,28,28,20,20,28,28,4,4,12,12,4,4,12,12,20,20,28,28,20,20,28,28}'::int[])[k+1], 
     y0 + ('{2,6,2,6,10,14,10,14,2,6,2,6,10,14,10,14,18,22,18,22,26,30,26,30,18,22,18,22,26,30,26,30}'::int[])[k+1], 
     32632
    ) as geom, orig_g
 from (
  select k, (l.gid<<5)+k as gid,
    l.jurisd_local_id, l.gid as parent_gid, l.l +25 as l,
    code_b32nvu||('{0,1,2,3,4,5,6,7,8,9,B,C,D,F,G,H,J,K,L,M,N,P,Q,R,S,T,U,V,W,X,Y,Z}'::text[])[k+1] as code_b32nvu, 
    st_xmin(l.geom) as x0, st_ymin(l.geom) as y0, l.geom as orig_g
  from (select * from prod.cm_logistic where l=130) l, generate_series(0,31) t0(k)
 ) t1
;

CREATE VIEW cm_logistic_vw_pts_l105to130 as
 select gid, jurisd_local_id, l, code_b32nvu, ST_POINT(
     x0 + ('{16,48,16,48,80,112,80,112,16,48,16,48,80,112,80,112,144,176,144,176,208,240,208,240,144,176,144,176,208,240,208,240}'::int[])[k+1], 
     y0 + ('{16,16,48,48,16,16,48,48,80,80,112,112,80,80,112,112,16,16,48,48,16,16,48,48,80,80,112,112,80,80,112,112}'::int[])[k+1], 
     32632
    ) as geom, orig_g
 from (
  select k, (l.gid<<5)+k as gid,
    l.jurisd_local_id, l.gid as parent_gid, l.l +25 as l,
    code_b32nvu||('{0,1,2,3,4,5,6,7,8,9,B,C,D,F,G,H,J,K,L,M,N,P,Q,R,S,T,U,V,W,X,Y,Z}'::text[])[k+1] as code_b32nvu, 
    st_xmin(l.geom) as x0, st_ymin(l.geom) as y0, l.geom as orig_g
  from (select * from prod.cm_logistic where l=105) l, generate_series(0,31) t0(k)
 ) t1
;
-- para saltar:
select array_agg(x_i/16 order by i) as x from unnest('{16,48,16,48,80,112,80,112,16,48,16,48,80,112,80,112,144,176,144,176,208,240,208,240,144,176,144,176,208,240,208,240}'::int[]) WITH ORDINALITY t(x_i,i);
--select gid, jurisd_local_id, l, code_b32nvu, ST_POINT(
--    x0 + ('{1,3,1,3,5,7,5,7,1,3,1,3,5,7,5,7,9,11,9,11,13,15,13,15,9,11,9,11,13,15,13,15}'::int[])[k+1], 
--    y0 + ('{1,1,3,3,1,1,3,3,5,5,7,7,5,5,7,7,1,1,3,3,1,1,3,3,5,5,7,7,5,5,7,7}'::int[])[k+1], 
--  where code_b32nvu like '5V98%' -- para não explodir o QGIS

Outra forma de calcular, nesse caso para chegar de L15.5 a L18. Um cuidado especial a ser dado pela função DeInterleave que, ao tratar de "refinemento anisotrópico", inverte i com j. No caso abaixo onde *s*=1 não houve necessidade de corrigir o "s*(0.5+i)".

select array_agg(xy[1]) xref, array_agg(xy[2]) yref
from (
 select *, natcod.hiddenBig_to_vBit(gid) as gid_vbit, vbit_DeInterleave_to_int(b'0'||k::bit(5)) as xy
 from (
  select (l.gid<<5)+k as gid, k,
    l.jurisd_local_id, l.gid as parent_gid, l.l +25 as l,
    code_b32nvu||('{0,1,2,3,4,5,6,7,8,9,B,C,D,F,G,H,J,K,L,M,N,P,Q,R,S,T,U,V,W,X,Y,Z}'::text[])[k+1] as code_b32nvu
    --,geom as orig_geom
  from (select * from prod.cm_logistic where l=155 and code_b32nvu like '977%' limit 1) l, generate_series(0,31) t0(k)
 ) t1
) t2;
-- xref={0,0,1,1,0,0,1,1,2,2,3,3,2,2,3,3,0,0,1,1,0,0,1,1,2,2,3,3,2,2,3,3}
-- yref={0,1,0,1,2,3,2,3,0,1,0,1,2,3,2,3,4,5,4,5,6,7,6,7,4,5,4,5,6,7,6,7}
-- swap(xref,yref) by vbit_DeInterleave_to_int(b'0'..)
-------------------------------------------
CREATE VIEW cm_logistic_vw_pts_l155to180 as
 select gid, jurisd_local_id, l, code_b32nvu, ST_POINT(
     x0 + 0.5 + ('{0,1,0,1,2,3,2,3,0,1,0,1,2,3,2,3,4,5,4,5,6,7,6,7,4,5,4,5,6,7,6,7}'::int[])[k+1], 
     y0 + 0.5 + ('{0,0,1,1,0,0,1,1,2,2,3,3,2,2,3,3,0,0,1,1,0,0,1,1,2,2,3,3,2,2,3,3}'::int[])[k+1], 
     32632
    ) as geom, orig_g
 from (
  select k, (l.gid<<5)+k as gid,
    l.jurisd_local_id, l.gid as parent_gid, l.l +25 as l,
    code_b32nvu||('{0,1,2,3,4,5,6,7,8,9,B,C,D,F,G,H,J,K,L,M,N,P,Q,R,S,T,U,V,W,X,Y,Z}'::text[])[k+1] as code_b32nvu, 
    st_xmin(l.geom) as x0, st_ymin(l.geom) as y0, l.geom as orig_g
  from (select * from prod.cm_logistic where l=155 limit 200) l, generate_series(0,31) t0(k)
 ) t1
;

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.

Dumps

Exemplos de dumps binários:

pg_dump -U postgres -f gridCaruaru_geomLatLong-2024-07-29.dump -Fc -t prod.br_logistic prod
pg_restore -U postgres -d prod /tmp/gridCaruaru_geomLatLong-2024-07-29.dump

Questão do dump/restore SQL: um simples NULL não carrega no psql. > pg_dump -a --inserts databasename > exportfilename.sql.

Precisa testar se funciona com pg_restore. Comparar com tamanho gz e ver se restore permite gz. Exemplo:

pg_dump -Z6 -h localhost test_db -f test_db.tar

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.