osmc:Metodologia/Algoritmo SQL/Lib: mudanças entre as edições
(→Células-folha virtuais: correcao) |
|||
(32 revisões intermediárias pelo mesmo usuário não estão sendo mostradas) | |||
Linha 6: | Linha 6: | ||
[[Arquivo:DNGS-coords1.png|miniaturadaimagem|440x440px|Sistemas de coordenadas previstos pelo padrão DNGS: LatLong, XY e IJ. A superfície ''S''<sup>2</sup> é o globo terrestre, ''R''<sup>2</sup> 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''<sup>2</sup>, mas mesma área em XY.]]<!-- mais detalhes em https://math.stackexchange.com/a/849521/70274 --> | [[Arquivo:DNGS-coords1.png|miniaturadaimagem|440x440px|Sistemas de coordenadas previstos pelo padrão DNGS: LatLong, XY e IJ. A superfície ''S''<sup>2</sup> é o globo terrestre, ''R''<sup>2</sup> 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''<sup>2</sup>, mas mesma área em XY.]]<!-- mais detalhes em https://math.stackexchange.com/a/849521/70274 --> | ||
O projeto como um todo faz uso de | 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]]. | *'''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: | * '''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).<br/>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''. | ||
** '''face_id''': chave primária de 12 bits | |||
** '''cell_id''': identificador local de célula, dado por ''bit string'' com tamanho variável. | * '''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. <br/>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 [[Código_Natural/Representação_interna#hInt|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. | 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. | ||
Linha 22: | Linha 23: | ||
<pre> | <pre> | ||
COORDINATES AND STANDARD PARAMETERS: | COORDINATES AND STANDARD PARAMETERS: | ||
canonical jiL: [cover_id, j,i, L] | canonical jiL: [cover_id, j,i, L] | ||
complete jiL: [cover_id, j,i, L, s, X0, Y0, is_half] | complete jiL: [cover_id, j,i, L, s, X0, Y0, is_half] | ||
Linha 29: | Linha 31: | ||
canonical xyS: [x,y,s,is_half] | canonical xyS: [x,y,s,is_half] | ||
complete xyS: [x,y,s,is_half, cover_id, j,i, L, X0, Y0] | 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) | 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. | |||
</pre> | </pre> | ||
* '''JI0''': cada ''face_id'' pode ser representado por índices ''j0'' e ''i0'' da grade ''L0'', que por sua vez podem levar à origem <math>(x0_{j0},y0_{i0})</math> de cada célula da cobertura, nas coordenadas da projeção SRID. Na prática o que temos também é <math>(x0_{face\_id},y0_{face\_id})</math> a origem em termos de ''face_id''. | * '''JI0''': cada ''face_id'' pode ser representado por índices ''j0'' e ''i0'' da grade ''L0'', que por sua vez podem levar à origem <math>(x0_{j0},y0_{i0})</math> de cada célula da cobertura, nas coordenadas da projeção SRID. Na prática o que temos também é <math>(x0_{face\_id},y0_{face\_id})</math> 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''.<br/> 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''. | * '''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''. | ||
Linha 43: | Linha 48: | ||
<!-- [[Arquivo:Matrix-JI-coventions.png|220px|center]] --> | <!-- [[Arquivo:Matrix-JI-coventions.png|220px|center]] --> | ||
Por fim, na convenção cartesiana ''XY'' a origem é (0,0) e a equação <math>X_i=(i-1)\cdot s</math> pode ser simplificada para <math>X_i=i\cdot s</math> 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'''''. | Por fim, na convenção cartesiana ''XY'' a origem é (0,0) e a equação <math>X_i=(i-1)\cdot s</math> pode ser simplificada para <math>X_i=i\cdot s</math> 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 [[wikipedia:Z-order curve|Curva-Z]] ilustrada abaixo, nascem na matriz com índices invertidos, ''JI'', e só [[wikipedia:Lazy evaluation|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". | |||
[[Arquivo:MortonCurve IJ-JI.png|600px|center]] | [[Arquivo:MortonCurve IJ-JI.png|600px|center]] | ||
Linha 53: | Linha 60: | ||
* A projeção é de igual-área, ou bastante próxima disso. | * 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'': <math>S_{Lmax}=1</math>. | * O sistema de grades de nível ''L'' é um mosaico quadrilátero, que tem o metro como maior resolução, no nível ''Lmax'': <math>S_{Lmax}=1</math>. | ||
* A grade de nível zero, ''L''=0, é um conjunto de até 16 quadrados de área <math>A_0=2^{2\cdot Lmax}</math>, que cobre todo o domínio territorial nacional. Portanto quadrados de lado <math>S_0=2^{Lmax}</math>, ou seja, obriga-se o uso de uma potência de 2 na cobertura nacional. | * A grade de nível zero, ''L''=0, é um conjunto de até 16 quadrados de área <math>A_0=2^{2\cdot Lmax}</math>, que cobre todo o domínio territorial nacional. Portanto quadrados de lado <math>S_0=2^{Lmax}</math>, ou seja, obriga-se o uso de uma [[Discrete_National_Grid_Systems/pt#Convenções_fixadas_pelos_requisitos|potência de 2 na cobertura nacional]]. | ||
* A resolução da grade de nível ''L'' é <math>S_{L}=2^{Lmax-L}</math>, interpretada como tamanho do lado da célula quadrada, de área <math>A_{L}={S_{L}}^2 =2^{2\cdot(Lmax-L)}</math>. | * A resolução da grade de nível ''L'' é <math>S_{L}=2^{Lmax-L}</math>, interpretada como tamanho do lado da célula quadrada, de área <math>A_{L}={S_{L}}^2 =2^{2\cdot(Lmax-L)}</math>. | ||
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, <math>16 \cdot 2^{2\cdot Lmax} \gtrapprox AN</math>, portanto <math>\lfloor 0,5\cdot log_2(AN) -2 \rfloor \gtrapprox Lmax \gtrapprox \lceil 0,5\cdot log_2(AN) -2 \rceil</math>. | 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, <math>16 \cdot 2^{2\cdot Lmax} \gtrapprox AN</math>, portanto <math>\lfloor 0,5\cdot log_2(AN) -2 \rfloor \gtrapprox Lmax \gtrapprox \lceil 0,5\cdot log_2(AN) -2 \rceil</math>. 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 <math>S_{L}</math> funciona também para níveis-meio, <math>Lhalf \in \{ \forall x | x = \lceil x \rceil - 0,5 \}</math>, 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 fórmula de <math>S_{L}</math> funciona também para níveis-meio, <math>Lhalf \in \{ \forall x | x = \lceil x \rceil - 0,5 \}</math>, 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, <math>S_{Lhalf}=\sqrt{A_{Lhalf}}</math>; e a interpretação da [[Generalized_Geohash/pt#Representação_geométrica|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'', <math>A_{Lhalf}=2\cdot{A_{Lprox}}</math>. Portanto, fazendo <math>Lprox={\lceil Lhalf\rceil}</math> temos: <br /> <math>S_{Lhalf}=\sqrt{2\cdot{A_{Lprox}}} = \sqrt{2} \cdot \sqrt{{{S_{\lceil Lhalf \rceil}}^2}} = 2^{Lmax-\lceil Lhalf\rceil + 0,5}</math>. | |||
A cada <code>face_id</code>, 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'': <br /> <math>nCells_L = (S_0/S_{L})^2 = (2^{Lmax-Lmax+L})^2=2^{2L}</math> <br/>É invariante por face e por país, e válido para "níveis meio". | A cada <code>face_id</code>, 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'': <br /> <math>nCells_L = (S_0/S_{L})^2 = (2^{Lmax-Lmax+L})^2=2^{2L}</math> <br/>É invariante por face e por país, e válido para "níveis meio". | ||
Linha 69: | Linha 77: | ||
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>{jiL}=({\text{id}},j,i)_L</math>. 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 <math>xy0_{\text{id}}=(x0,y0)_{\text{id}}</math> da origem XY0 da célula. 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}=({id},j,i)_L=({id},\lfloor\frac{x-x0}{S_{L}}\rfloor ~,~ \lfloor\frac{y-y0}{S_{L}}\rfloor)_L</math> | |||
** No nível ''L'' meio (''Lhalf'') temos retângulos de lados <math>S_{Lprox}</math> no eixo Y e lados <math>2\cdot S_{Lprox}</math> no eixo X, onde <math>{Lprox}=\lceil Lhalf \rceil</math>. A quebra em partes inteiras resultará em <math>{jiL}=({\text{id}},j,i)_{Lhalf}=({id}, \lfloor \frac{x-x0}{2\cdot S_{Lprox}} \rfloor ~,~ \lfloor \frac{y-y0}{S_{Lprox}} \rfloor)_{Lhalf}</math> | |||
* "jiL to xyLRef": a reconstrução é direta, mas vai apresentar fórmulas diferentes para o quadrado e o retângulo. | |||
** ''L'' inteiro: <math>{xyLRef}=(j\cdot S_{L}+x0 ~,~ i\cdot S_{L}+y0)</math> | |||
** ''L'' meio: <math>{xyLRef}=(j\cdot 2\cdot S_{Lprox}+x0 ~,~ i\cdot S_{Lprox}+y0)</math> | |||
Devido à diferença de cálculo entre ''L'' inteiro e meio, e ao fato de <math>2^x</math> poder ser otimizado por rotação de bits, <code><nowiki>2<<(x-1)</nowiki></code>, as funções otimizadas adotam ''flag'' <code>is_half</code> 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== | ==Core Lib== | ||
Linha 120: | Linha 139: | ||
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. | 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. | ||
<syntaxhighlight lang="sql" style="font-size: 80%;"> | <syntaxhighlight lang="sql" style="font-size: 80%;"> | ||
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 | ||
$f$ LANGUAGE SQL IMMUTABLE; | $f$ LANGUAGE SQL IMMUTABLE; | ||
Linha 143: | Linha 147: | ||
SELECT array_agg( vbit_to_int(x_i) ORDER BY i) FROM unnest(x) WITH ORDINALITY t(x_i,i) | SELECT array_agg( vbit_to_int(x_i) ORDER BY i) FROM unnest(x) WITH ORDINALITY t(x_i,i) | ||
$f$ LANGUAGE SQL IMMUTABLE; | $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$ | CREATE FUNCTION vbit_interleave(x varbit, y varbit) returns varbit as $f$ | ||
SELECT bitcat_agg( substring(x,i,1) || substring(y,i,1) ) | 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; | $f$ LANGUAGE SQL IMMUTABLE; | ||
CREATE FUNCTION vbit_DeInterleave(x varbit) returns varbit[] as $f$ | CREATE FUNCTION vbit_DeInterleave(x varbit) returns varbit[] as $f$ | ||
SELECT array[ | 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) | FROM generate_series(1,bit_length(x),2) t(i) | ||
$f$ LANGUAGE SQL IMMUTABLE; | $f$ LANGUAGE SQL IMMUTABLE; | ||
CREATE FUNCTION | CREATE FUNCTION vbit_DeInterleave_to_ints(x varbit) returns int[] as $wrap$ | ||
SELECT vbit_to_int( vbit_DeInterleave(x) ) | SELECT vbit_to_int( vbit_DeInterleave(x) ) | ||
$wrap$ LANGUAGE SQL IMMUTABLE; | $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 | |||
</syntaxhighlight> | </syntaxhighlight> | ||
Linha 171: | Linha 184: | ||
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 (XY). | 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 (XY). | ||
A inversa, <code>vbit_DeInterleave</code> 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"[https://web.archive.org/web/20240811233204/https://mfem.org/howto/ncmesh/ ref]. | |||
=== Otimização do cálculo de cobertura=== | === Otimização do cálculo de cobertura=== | ||
Linha 348: | Linha 363: | ||
------ | ------ | ||
== 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: | |||
<syntaxhighlight lang="sql" style="font-size: 80%;"> | |||
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 | |||
; | |||
</syntaxhighlight> | |||
<pre> | |||
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% | |||
</pre> | |||
* 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 [https://www.poder360.com.br/governo/plano-de-guedes-tem-extincao-de-769-municipios-e-elimina-quase-20-mil-cargos/ 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 <code>isolabel_ext like 'CM-%-%'</code> e seu Lmax é 18, portanto <code>ceil( (2*18.0-log(area)/log(2))/2.0 )</code>. 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''. <!-- | |||
L5.5 = cobre Soa | |||
L8 = Side: 1.02 km | |||
L10.5 = Side: 181.02 m | |||
L13 = Side: 32 m | |||
L15.5 = Side: 5.66 m | |||
L18 = 1m² | |||
WITH s AS ( | |||
SELECT *, ceil( (36.0-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 'CM-%-%' | |||
) t0 | |||
), t AS ( | |||
SELECT COUNT(*) ntot, sum(area) areatot, sum(perim) prmtot FROM s | |||
) | |||
SELECT CASE WHEN level >=5 THEN '5.5' when level<3 THEN '~'||level::text ELSE '3' 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 | |||
; | |||
--><pre> | |||
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% | |||
</pre> | |||
A Colômbia (<code>CO</code>), 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''). | |||
<syntaxhighlight lang="sql" style="font-size: 80%;"> | |||
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] | |||
</syntaxhighlight> | |||
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. | |||
<syntaxhighlight lang="sql" style="font-size: 80%;"> | |||
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 | |||
</syntaxhighlight> | |||
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)". | |||
<syntaxhighlight lang="sql" style="font-size: 80%;"> | |||
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 | |||
; | |||
</syntaxhighlight> | |||
==Testes da lib== | ==Testes da lib== | ||
Primeiro demonstrações (de que as funções funcionam), testes para validar hipóteses e requisitos básicos; depois [https://pt.stackoverflow.com/a/13530/4186 testes de regressão]. | Primeiro demonstrações (de que as funções funcionam), testes para validar hipóteses e requisitos básicos; depois [https://pt.stackoverflow.com/a/13530/4186 testes de regressão]. | ||
===Dumps === | |||
Exemplos de dumps binários: | |||
<pre> | |||
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 | |||
</pre> | |||
Questão do dump/restore SQL: um simples NULL não carrega no psql. <code>> pg_dump -a --inserts databasename > exportfilename.sql</code>. | |||
Precisa testar se funciona com pg_restore. Comparar com tamanho gz e ver se restore permite gz. Exemplo: | |||
<pre> | |||
pg_dump -Z6 -h localhost test_db -f test_db.tar | |||
</pre> | |||
===Validação primária=== | ===Validação primária=== | ||
Testes para validação, usando funções ou tabelas já homologadas: | Testes para validação, usando funções ou tabelas já homologadas: |
Edição atual tal como às 21h25min de 15 de agosto 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 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 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.
- 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 , í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.
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: .
- 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 . 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 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", com valor obtido pela raiz quadrada da área do retângulo, ; 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, . 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: 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 da origem XY0 da célula. 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 (Lhalf) temos retângulos de lados no eixo Y e lados no eixo X, onde . A quebra em partes inteiras resultará em
- "jiL to xyLRef": a reconstrução é direta, mas vai apresentar fórmulas diferentes para o quadrado e o retângulo.
- L inteiro:
- L meio:
Devido à diferença de cálculo entre L inteiro e meio, e ao fato de 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
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á 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).
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
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;
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.