osmc:Metodologia/Algoritmo SQL/Lib: mudanças entre as edições
(22 revisões intermediárias pelo mesmo usuário não estão sendo mostradas) | |||
Linha 1: | Linha 1: | ||
Biblioteca para simplificar a geração do sistema de grades de um país. Exemplificando implementação dentro do ''schema'' <code>grid_br</code> das grades do Brasil. | |||
Lembretes: | |||
* Convencionou-se em GGeohash usar YX no lugar de XY. | |||
* ... | |||
A seguir avaliando o uso default de XYref no lugar de XYcenter, para o desenho da célula. | 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 == | |||
O projeto como um todo faz uso de 3 sistemas de coordenadas | |||
* LatLong: WGS84 no formato [[Geo URI]] | |||
* XY: cartesiano escolar, coordenadas planas da projeção igual-area oficial do país | |||
* IJ: XY discretizado | |||
* JI? A curva de morton rotacional IJ | |||
... | |||
== Core == | |||
'''Quantizadores''': transformam as coordenadas contínuas YX, 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 YX é o tamanho de lado ''S'' (''side size'') da célula de nível ''L''. | |||
Por imposição do [[Discrete National Grid Systems/pt|padrão DNGS]] temos <math>S_{L}=2^{Lmax-L}</math>. No caso do Brasil ''Lmax''=20, no caso de Camarões ''Lmax=18''. | |||
:Nota. A fórmula de ''S'' funciona também para níveis-meio, por exemplo ''L''=1.5. Isso se deve à suposição de que o "tamanho do lado genérico ''S'' de um retângulo" seja a raiz quadrada da área do retângulo; e pela [[Generalized_Geohash/pt#Representação_geométrica|construção geométrica dos níveis-meio]], cujas células (necessariamente de áreas iguais) são a união de 2 células do próximo nível inteiro: <br/> <math>S_{Lhalf}=\sqrt{2\cdot{Area_{\lceil Lhalf\rceil}}} = \sqrt{2} \cdot \sqrt{{{S_{\lceil Lhalf\rceil}}^2}} = 2^{Lmax-\lceil Lhalf\rceil + 0.5}</math> onde <math>Lhalf = \forall L | L = \lceil L \rceil - 0.5</math>. <br/>A mesma 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>. | |||
A decisão de usar a sequência YX e não XY precisa talvez ser revista. Na cultura escolar brasileira XY é horizontal-vertical. Na cultura das imagens de satélite e geoprocessamento XY é vertical-horizontal. Adotamos a "cultura PostGIS", das funções ''standard spatial type'' (aquelas com prefixo "ST_"). | |||
<syntaxhighlight lang="sql" style="font-size: 80%;"> | <syntaxhighlight lang="sql" style="font-size: 80%;"> | ||
-- Core functions: | |||
drop FUNCTION if exists grid_br.yxS_collapseTo_ijS(int,int,int,boolean); | |||
drop FUNCTION if exists grid_br.yxL_collapseTo_ijL(int[]); | |||
CREATE FUNCTION grid_br.yxS_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.yxL_collapseTo_ijL(yxL 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; | |||
</syntaxhighlight> | </syntaxhighlight> | ||
'''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. | 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. | ||
<syntaxhighlight lang="sql" style="font-size: 80%;"> | <syntaxhighlight lang="sql" style="font-size: 80%;"> | ||
DROP FUNCTION if exists grid_br.xyS_draw_anycell(int,int,int) | DROP FUNCTION if exists grid_br.xyS_draw_anycell(int,int,int) | ||
; | ; | ||
CREATE FUNCTION grid_br.yxS_draw_anycell( | |||
CREATE FUNCTION grid_br. | |||
y int, -- Yref | y int, -- Yref | ||
x int, -- Xref (XY=canto inferior esquerdo) | x int, -- Xref (XY=canto inferior esquerdo) | ||
Linha 108: | Linha 61: | ||
IS 'Draws a square-cell using the requested ref-point, with the requested side size.' | IS 'Draws a square-cell using the requested ref-point, with the requested side size.' | ||
; | ; | ||
</syntaxhighlight> | </syntaxhighlight> | ||
Linha 123: | Linha 73: | ||
drop function if exists vbit_DeInterleave; | drop function if exists vbit_DeInterleave; | ||
; | DROP AGGREGATE if exists bitcat_agg(varbit); | ||
CREATE AGGREGATE bitcat_agg(varbit) ( | |||
STYPE = varbit | |||
,SFUNC = bitcat -- \do+ "||" | |||
); | |||
CREATE FUNCTION vbit_to_int(x varbit) returns int as $f$ | CREATE FUNCTION vbit_to_int(x varbit) returns int as $f$ | ||
SELECT ( substring(0::bit(32),bit_length(x)+1) || x )::bit(32)::int | SELECT ( substring(0::bit(32),bit_length(x)+1) || x )::bit(32)::int | ||
Linha 137: | Linha 92: | ||
CREATE FUNCTION vbit_interleave(x varbit, y varbit) returns varbit as $f$ | CREATE FUNCTION vbit_interleave(x varbit, y varbit) returns varbit as $f$ | ||
SELECT | SELECT bitcat_agg( substring(x,i,1) || substring(y,i,1) ) | ||
FROM generate_series(1,bit_length(x)) t(i) | FROM generate_series(1,bit_length(x)) t(i) | ||
$f$ LANGUAGE SQL IMMUTABLE; | $f$ LANGUAGE SQL IMMUTABLE; | ||
Linha 163: | Linha 118: | ||
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. | ||
==Demais convenções== | A degeneração geométrica, de quadrado para retângulo, é relativa ao segundo argumento de <code>vbit_interleave(x,y)</code>. Como a função é sempre chamada com a mesma ordem dos argumentos, sempre teremos ou só retangulos orizontais (XY) ou só verticais (YX). | ||
=== Otimização do cálculo de cobertura === | |||
Ver [[osmc:Metodologia/Algoritmo_SQL#Tratamento_das_configurações|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'' na projeção do Brasil, ou seja, com coordenadas planas ''x'' e ''y''. O primeiro passo é obter as coordenadas ''ij0'' do ponto na cobertura. | |||
O tamanho de célula ''s0'' da cobertura nacional é o ''default'' em <code>ij0=grid_br.xyS_collapseTo_ijS(x,y)</code>. Alternativamente uma cadeia de ''IF''s (árvore de decisão BBOX) pode nos fornecer rapidamente o valor. O que é mais rápido, depende de ''benchmark'' em cada linguagem (C ou SQL). | |||
Como ij0 é um caso muito especial de ''ij'' podemos obter o seu valor num formato mais conveniente, por exemplo usando 3 bits (0=000, 1=001, 2=010, 3=011, 4=100, 5=101) e concatenando os valores de ''i'' e ''j'': 40=b'100000'=32 41=b'100001'=33 42=b'100010'=34 ... 13=b'001011'=11 02=b'000010'=2 44=b'100100'=36 24=b'010100'=20. <!-- | |||
select varbit_to_int(b'100000') as "40", varbit_to_int(b'100001') as "41", | |||
varbit_to_int(b'100010') as "42", varbit_to_int(b'001011') as "13", | |||
varbit_to_int(b'000010') as "02", varbit_to_int(b'100100') as "44"; | |||
40 | 41 | 42 | 13 | 02 | 44 | |||
----+----+----+----+----+---- | |||
32 | 33 | 34 | 11 | 2 | 36 | |||
--> Resulta numa array esparsa de 40 posições. Teremos: | |||
* coordenadas iniciais ''x0'' e ''y0'' dadas por ''arrays'' <code>x0_from_ij0</code> e <code>x0_from_ij0</code>. | |||
* índice cbits0 dado por ''array'' <code>cbits0_from_ij0</code>. | |||
Por fim as 'arrays'' podem ser "hardcoded" nas funções usuárias, otimizando a obtenção dos valores desejado. | |||
=== 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: | |||
#* <code>ij0=grid_br.xyS_collapseTo_ijS(x,y);</code> com ''s0'' é default, ou função otimizada. | |||
#* <code>x0=x0_from_ij0[ij0]; y0=y0_from_ij0[ij0];</code> ou função de ij0 retorando xy0. | |||
#* <code>cbits0 = grid_br.IJ0_to_vbitL0( ij0, false )</code> | |||
# Código ''cbits'' e geometria da célula do nível ''L'': | |||
#* <code>ijL=grid_br.xyL_collapseTo_ijL(x-xy0[1], y-xy0[2], L);</code> | |||
#* <code>cbits = cbits0 || ints_to_interleavedbits(ijL)</code> | |||
#* <code>grid_br.xyS_draw_anycell( grid_br.ijL_to_xyL(ijL) )</code> | |||
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: | |||
:<code>SE lenght(cbits0)>4 e level_desejado<1.5 THEN NULL; ELSE recalcula xy0 dentro da célula especial.</code> | |||
<syntaxhighlight lang="sql" style="font-size: 80%;"> | |||
drop FUNCTION if exists grid_br.xyS_to_cbits(int,int,int,boolean) | |||
; | |||
-- revisar se é YX! | |||
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; | |||
</syntaxhighlight> | |||
===Demais convenções=== | |||
... | ... | ||
Linha 176: | Linha 191: | ||
No caso especial de fronteira nacional, como as células "0" ou "3", onde não existem "células contíguas", o dono é o país estrangeiro, portanto, a rigor, a linha não existe na grade. | No caso especial de fronteira nacional, como as células "0" ou "3", onde não existem "células contíguas", o dono é o país estrangeiro, portanto, a rigor, a linha não existe na grade. | ||
== | === Geração do sistema de grades === | ||
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: | |||
* [[wikipedia:Field (geography)|''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 [https://postgis.net/docs/raster.html raster] ou [https://postgis.net/docs/geomval.html 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. | |||
Algorítmo genérico, para gerar grades nacionais ou municipais. O algoritmo de <code>grid_br.parents_to_children_cuting</code> é simples, pode ser expresso como recorê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. | |||
<syntaxhighlight lang="sql" style="font-size: 80%;"> | <syntaxhighlight lang="sql" style="font-size: 80%;"> | ||
select | -- 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$; | |||
</syntaxhighlight> | </syntaxhighlight> | ||
<syntaxhighlight lang="sql"> | == Helper lib == | ||
select | 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''). | ||
<syntaxhighlight lang="sql" style="font-size: 80%;"> | |||
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(yxl int[]) RETURNS int[] AS $wrap$ | |||
SELECT grid_br.xyS_collapseTo_ijS(yxl[1], yxl[2], yxl[3]) | |||
$wrap$ LANGUAGE SQL IMMUTABLE | |||
; | |||
----- | |||
drop FUNCTION if exists grid_br.ijS_to_xySref(int,int,int,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(yxl int[], is_half boolean default false) RETURNS int[] AS $wrap$ | |||
SELECT grid_br.ijS_to_xySref(yxl[1], yxl[2], yxl[3],is_half) | |||
$wrap$ LANGUAGE SQL IMMUTABLE; | |||
</syntaxhighlight> | |||
------ | |||
== 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]. | |||
=== Validação primária === | |||
Testes para validação, usando funções ou tabelas já homologadas: | |||
<syntaxhighlight lang="sql" style="font-size: 80%;"> | |||
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; | |||
</syntaxhighlight> | |||
Grade da Colômbia, de zero a "e": | |||
<pre> | |||
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 | |||
</pre> | |||
[[Arquivo:Grid br-L0-xyref-points.png|miniaturadaimagem|420x420px]] | |||
Conferindo a localização do ponto xyRef: | |||
<syntaxhighlight lang="sql" style="font-size: 80%;"> | |||
drop view test_l0_refpoint; | |||
create view test_l0_refpoint as | |||
select *, ST_SetSRID(ST_POINT(yxs[2], yxs[1]), 952019) as geom | |||
from ( | |||
select code_b16h, ijs[1] i, ijs[2] j, grid_br.ijS_to_xySref(ijs) as yxs, count(*) n | |||
from ( | |||
select a.code_b16h, grid_br.xylevel_collapseTo_ijS(st_y(pt.geom)::int,st_x(pt.geom)::int,hlevel::int) as ijs | |||
from grid_br.all_levels a, lateral ST_DumpPoints( ST_GeneratePoints(a.geom, 20) ) as pt | |||
where hlevel=0 and code_b16h !~ '^f' | |||
) t1 group by 1,2,3,4 order by 1 | |||
) t2; | |||
</syntaxhighlight> | </syntaxhighlight> | ||
Portanto <code>grid_br.ijS_to_xySref(grid_br.xylevel_collapseTo_ijS(ponto,nivel))</code> fornece automaticamente a célula ''L0'' de referência onde será aplicado algoritmo GGeohash. |
Edição das 08h42min de 2 de julho de 2024
Biblioteca para simplificar a geração do sistema de grades de um país. Exemplificando implementação dentro do schema grid_br
das grades do Brasil.
Lembretes:
- Convencionou-se em GGeohash usar YX no lugar de XY.
- ...
A seguir avaliando o uso default de XYref no lugar de XYcenter, para o desenho da célula. 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
O projeto como um todo faz uso de 3 sistemas de coordenadas
- LatLong: WGS84 no formato Geo URI
- XY: cartesiano escolar, coordenadas planas da projeção igual-area oficial do país
- IJ: XY discretizado
- JI? A curva de morton rotacional IJ
...
Core
Quantizadores: transformam as coordenadas contínuas YX, 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 YX é o tamanho de lado S (side size) da célula de nível L.
Por imposição do padrão DNGS temos . No caso do Brasil Lmax=20, no caso de Camarões Lmax=18.
- Nota. A fórmula de S funciona também para níveis-meio, por exemplo L=1.5. Isso se deve à suposição de que o "tamanho do lado genérico S de um retângulo" seja a raiz quadrada da área do retângulo; e pela construção geométrica dos níveis-meio, cujas células (necessariamente de áreas iguais) são a união de 2 células do próximo nível inteiro:
onde .
A mesma 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: .
A decisão de usar a sequência YX e não XY precisa talvez ser revista. Na cultura escolar brasileira XY é horizontal-vertical. Na cultura das imagens de satélite e geoprocessamento XY é vertical-horizontal. Adotamos a "cultura PostGIS", das funções standard spatial type (aquelas com prefixo "ST_").
-- Core functions:
drop FUNCTION if exists grid_br.yxS_collapseTo_ijS(int,int,int,boolean);
drop FUNCTION if exists grid_br.yxL_collapseTo_ijL(int[]);
CREATE FUNCTION grid_br.yxS_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.yxL_collapseTo_ijL(yxL 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.yxS_draw_anycell(
y int, -- Yref
x int, -- Xref (XY=canto inferior esquerdo)
s int -- side size
) RETURNS geometry AS $f$
SELECT CASE WHEN s=0 THEN NULL ELSE ST_GeomFromText( format(
'POLYGON((%s %s,%s %s,%s %s,%s %s,%s %s))',
x,y, x+s,y, x+s,y+s, x,y+s, x,y
), 952019) END
$f$ LANGUAGE SQL IMMUTABLE;
COMMENT ON FUNCTION grid_br.xyS_draw_anycell(int,int,int)
IS 'Draws a square-cell using the requested ref-point, with the requested side size.'
;
Construtor do identificador cbits
Existem várias formas de obter o identificador baseado em curva de morton. O identificador IJ da célula, como vimos, é fácil de ser obtido no GGeohash.
drop function if exists vbit_to_int(varbit);
drop function if exists vbit_to_int(varbit[]);
drop function if exists vbit_interleave;
drop function if exists ints_to_interleavedbits;
drop function if exists vbit_DeInterleave;
DROP AGGREGATE if exists bitcat_agg(varbit);
CREATE AGGREGATE bitcat_agg(varbit) (
STYPE = varbit
,SFUNC = bitcat -- \do+ "||"
);
CREATE FUNCTION vbit_to_int(x varbit) returns int as $f$
SELECT ( substring(0::bit(32),bit_length(x)+1) || x )::bit(32)::int
$f$ LANGUAGE SQL IMMUTABLE;
CREATE FUNCTION vbit_to_bigint(x varbit) returns int as $f$
SELECT ( substring(0::bit(64),bit_length(x)+1) || x )::bit(64)::bigint
$f$ LANGUAGE SQL IMMUTABLE;
CREATE FUNCTION vbit_to_int(x varbit[]) returns int[] as $f$
SELECT array_agg( vbit_to_int(x_i) ORDER BY i) FROM unnest(x) WITH ORDINALITY t(x_i,i)
$f$ LANGUAGE SQL IMMUTABLE;
CREATE FUNCTION vbit_interleave(x varbit, y varbit) returns varbit as $f$
SELECT bitcat_agg( substring(x,i,1) || substring(y,i,1) )
FROM generate_series(1,bit_length(x)) t(i)
$f$ LANGUAGE SQL IMMUTABLE;
CREATE FUNCTION ints_to_interleavedbits(x int, y int, len int default 32) returns varbit as $f$
SELECT vbit_interleave( substring(x::bit(32),33-len), substring(y::bit(32),33-len) )
$f$ LANGUAGE SQL IMMUTABLE;
-- SELECT ints_to_interleavedbits(47,19,6); -- 100110101111
CREATE FUNCTION vbit_DeInterleave(x varbit) returns varbit[] as $f$
SELECT array[ string_agg(substring(x,i,1)::text,'')::varbit,
string_agg(substring(x,i+1,1)::text,'')::varbit
]
FROM generate_series(1,bit_length(x),2) t(i)
$f$ LANGUAGE SQL IMMUTABLE;
-- SELECT vbit_to_int(vbit_DeInterleave(b'100110101111')); -- 47, 19
CREATE FUNCTION vbit_DeInterleave_to_int(x varbit) returns int[] as $wrap$
SELECT vbit_to_int( vbit_DeInterleave(x) )
$wrap$ LANGUAGE SQL IMMUTABLE;
Os indexadores IJ são a forma mais compacta (normalizada) de representação da grade, e independente da projeção ou métrica usadas. A função ints_to_interleavedbits(i,j)
nos fornece exatamente o indexador desejado.
Cada célula da cobertura L0 já tem seu código de 4 bits (dígito hexadecimal), em seguida qualquer outra célula de qualquer nível hierárquico terá um código cbits do nível L com L*2 bits obtidos das coordenadas IJ através da função ints_to_interleavedbits(i,j)
. O tamanho de lado da célula do nível L, no caso do Brasil, será metros. A célula de 1 m é justamente a célula do nível L20, portanto terá bits de comprimento, concatenados aos 4 bits de L0 e ao prefixo do país.
A degeneração geométrica, de quadrado para retângulo, é relativa ao segundo argumento de vbit_interleave(x,y)
. Como a função é sempre chamada com a mesma ordem dos argumentos, sempre teremos ou só retangulos orizontais (XY) ou só verticais (YX).
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 na projeção do Brasil, ou seja, com coordenadas planas x e y. O primeiro passo é obter as coordenadas ij0 do ponto na cobertura.
O tamanho de célula s0 da cobertura nacional é o default em ij0=grid_br.xyS_collapseTo_ijS(x,y)
. Alternativamente uma cadeia de IFs (árvore de decisão BBOX) pode nos fornecer rapidamente o valor. O que é mais rápido, depende de benchmark em cada linguagem (C ou SQL).
Como ij0 é um caso muito especial de ij podemos obter o seu valor num formato mais conveniente, por exemplo usando 3 bits (0=000, 1=001, 2=010, 3=011, 4=100, 5=101) e concatenando os valores de i e j: 40=b'100000'=32 41=b'100001'=33 42=b'100010'=34 ... 13=b'001011'=11 02=b'000010'=2 44=b'100100'=36 24=b'010100'=20. Resulta numa array esparsa de 40 posições. Teremos:
- coordenadas iniciais x0 e y0 dadas por arrays
x0_from_ij0
ex0_from_ij0
. - índice cbits0 dado por array
cbits0_from_ij0
.
Por fim as 'arrays podem ser "hardcoded" nas funções usuárias, otimizando a obtenção dos valores desejado.
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 é YX!
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
As bordas, ou seja, a linha divisória entre duas células, não pode ser uma "linha sem dono", precisa pertencer a uma das células de um dos lados da borda. Por convenção DNGS as linhas afetadas pelo ponto de referência pertencem à célula referida.
Na ilsutração mostramos que dois lados da célula permanecem com o rótulo da referência, e outros dois lados pertencem às respectivas células contíguas. Por exemplo a célula "6" da ilustração possui dois lados próprios, o inferior e o esquerdo, e dois lados impróprios, das células "2" e "7".
No caso especial de fronteira nacional, como as células "0" ou "3", onde não existem "células contíguas", o dono é o país estrangeiro, portanto, a rigor, a linha não existe na grade.
Geração do sistema de grades
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.
Algorítmo genérico, para gerar grades nacionais ou municipais. O algoritmo de grid_br.parents_to_children_cuting
é simples, pode ser expresso como recorê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$;
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(yxl int[]) RETURNS int[] AS $wrap$
SELECT grid_br.xyS_collapseTo_ijS(yxl[1], yxl[2], yxl[3])
$wrap$ LANGUAGE SQL IMMUTABLE
;
-----
drop FUNCTION if exists grid_br.ijS_to_xySref(int,int,int,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(yxl int[], is_half boolean default false) RETURNS int[] AS $wrap$
SELECT grid_br.ijS_to_xySref(yxl[1], yxl[2], yxl[3],is_half)
$wrap$ LANGUAGE SQL IMMUTABLE;
Testes da lib
Primeiro demonstrações (de que as funções funcionam), testes para validar hipóteses e requisitos básicos; depois testes de regressão.
Validação primária
Testes para validação, usando funções ou tabelas já homologadas:
select code_b16h, ijs[1] i, ijs[2] j, count(*) n from (
select a.code_b16h, grid_br.xylevel_collapseTo_ijS(st_y(pt.geom)::int,st_x(pt.geom)::int,hlevel::int) as ijs
from grid_br.all_levels a, lateral ST_DumpPoints( ST_GeneratePoints(a.geom, 200) ) as pt
where hlevel=0
) t group by 1,2,3 order by 1;
Grade da Colômbia, de zero a "e":
grid_l0_cell_ij: 40 41 42 43 30 31 32 33 34 21 22 23 11 12 13 grid_l0_cell_b16: 0 1 2 3 4 5 6 7 8 9 a b c d e code_b16h | i | j | n -----------+---+---+----- 0 | 4 | 0 | 200 1 | 4 | 1 | 200 2 | 4 | 2 | 200 3 | 4 | 3 | 200 4 | 3 | 0 | 200 5 | 3 | 1 | 200 6 | 3 | 2 | 200 7 | 3 | 3 | 200 8 | 3 | 4 | 200 9 | 2 | 1 | 200 a | 2 | 2 | 200 b | 2 | 3 | 200 c | 1 | 1 | 200 d | 1 | 2 | 200 e | 1 | 3 | 200
Conferindo a localização do ponto xyRef:
drop view test_l0_refpoint;
create view test_l0_refpoint as
select *, ST_SetSRID(ST_POINT(yxs[2], yxs[1]), 952019) as geom
from (
select code_b16h, ijs[1] i, ijs[2] j, grid_br.ijS_to_xySref(ijs) as yxs, count(*) n
from (
select a.code_b16h, grid_br.xylevel_collapseTo_ijS(st_y(pt.geom)::int,st_x(pt.geom)::int,hlevel::int) as ijs
from grid_br.all_levels a, lateral ST_DumpPoints( ST_GeneratePoints(a.geom, 20) ) as pt
where hlevel=0 and code_b16h !~ '^f'
) t1 group by 1,2,3,4 order by 1
) t2;
Portanto grid_br.ijS_to_xySref(grid_br.xylevel_collapseTo_ijS(ponto,nivel))
fornece automaticamente a célula L0 de referência onde será aplicado algoritmo GGeohash.