Documentação

osmc:Metodologia/Algoritmo SQL/Lib: mudanças entre as edições

osmc Discussão

osmc:Metodologia/Algoritmo SQL/Lib (ver código-fonte)

Edição das 08h42min de 2 de julho de 2024

5 884 bytes adicionados ,  Terça-feira às 08h42min
m
→‎Algoritmo e funções finais de resolução
(10 revisões intermediárias pelo mesmo usuário não estão sendo mostradas)
Linha 1: Linha 1:
== Lib de Grades ==
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.


Biblioteca para simplificar a geração de grades. Exemplificando implementação dentro do ''schema'' <code>grid_br</code> das grades do Brasil. Convencionou-se em GGeohash usar YX no lugar de XY.
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''.  
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''.  


===  Core ===
== Sistemas de coordenadas ==
'''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 &mdash; ''grid hierarchical level'' do sistema de grades &mdash; 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''.  
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 &mdash; ''grid hierarchical level'' do sistema de grades &mdash; 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>.
: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_").
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_").
Linha 14: Linha 26:
<syntaxhighlight lang="sql" style="font-size: 80%;">
<syntaxhighlight lang="sql" style="font-size: 80%;">
-- Core functions:
-- Core functions:
drop FUNCTION if exists grid_br.xyS_collapseTo_ijS(int,int,int,boolean);
drop FUNCTION if exists grid_br.yxS_collapseTo_ijS(int,int,int,boolean);
drop FUNCTION if exists grid_br.xyL_collapseTo_ijL(int,int,int,boolean);
drop FUNCTION if exists grid_br.yxL_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$
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 ]
   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;
$f$ language SQL IMMUTABLE;


CREATE FUNCTION grid_br.xyL_collapseTo_ijL(y int, x int, intlevel int default 0) RETURNS int[] AS $f$
CREATE FUNCTION grid_br.yxL_collapseTo_ijL(yxL int[]) RETURNS int[] AS $f$
   SELECT array[ $ijS[1], $ijS[2], intlevel ]
   SELECT array[ $ijS[1], $ijS[2], intlevel ]
   FROM ( SELECT grid_br.xyS_collapseTo_ijS($1, $2, (2^(20-intlevel/10.0))::int, (intlevel%2)=1) ) t(ijS)
   FROM (  
$f$ language SQL IMMUTABLE
    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>


'''Construtores de geometria''':
'''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.
Linha 34: Linha 48:
DROP  FUNCTION if exists grid_br.xyS_draw_anycell(int,int,int)
DROP  FUNCTION if exists grid_br.xyS_draw_anycell(int,int,int)
;
;
DROP  FUNCTION if exists grid_br.xyS_draw_anycell(int[])
CREATE FUNCTION grid_br.yxS_draw_anycell(
;
CREATE FUNCTION grid_br.xyS_draw_anycell(
   y int,  -- Yref
   y int,  -- Yref
   x int,  -- Xref  (XY=canto inferior esquerdo)
   x int,  -- Xref  (XY=canto inferior esquerdo)
Linha 49: 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.'
;
;
CREATE FUNCTION grid_br.xyS_draw_anycell(int[]) RETURNS geometry AS $wrap$
  SELECT grid_br.xyS_draw_anycell($1[1],$1[2],$1[3])
$wrap$ LANGUAGE SQL IMMUTABLE;
</syntaxhighlight>
</syntaxhighlight>


Linha 110: Linha 119:


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).
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 e funções finais de resolução ===
Algoritmo principal:
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>


# ''pt'' na projeção do Brasil, ou seja, com coordenadas planas ''x'' e ''y''.
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>ij0=grid_br.xyS_collapseTo_ijS(x,y);  xy0=grid_br.ijS_to_xySref(ij0);</code>
:<code>SE lenght(cbits0)>4 e level_desejado<1.5 THEN NULL;  ELSE recalcula xy0 dentro da célula especial.</code>
# <code>cbits0 = grid_br.IJ0_to_L0( ij0, false )</code>
# <code>cbits = cbits0 || ints_to_interleavedbits(x-x0, y-y0, 20-level_desejado)</code>
 
É um pouco mais complexo, para contemplar os casos de pontos sobre cobertura fantasma. Existe uma condição de validade e um ajuste do ponto ao nível:
:<code>SE lenght(cbits0)>4 e level_desejado<1.5 THEN NULL;  ELSE recalcula xy0 dentro da célula.</code>


<syntaxhighlight lang="sql" style="font-size: 80%;">
<syntaxhighlight lang="sql" style="font-size: 80%;">
drop FUNCTION if exists grid_br.xyS_to_cbits(int,int,int,boolean)
drop FUNCTION if exists grid_br.xyS_to_cbits(int,int,int,boolean)
;
;
-- cuidado é YX! ver subtração
-- revisar se é YX!
CREATE FUNCTION grid_br.xyS_to_cbits(
CREATE FUNCTION grid_br.xyS_to_cbits(
   y int,
   y int,
Linha 134: Linha 170:
-- falta conferir se lenght(cbits0)>4 AND s<3
-- falta conferir se lenght(cbits0)>4 AND s<3
  SELECT CASE WHEN use_country_prefix THEN b'0001001100' ELSE b'' END  
  SELECT CASE WHEN use_country_prefix THEN b'0001001100' ELSE b'' END  
         || cbits0 || ints_to_interleavedbits(x-xy0[1], y-xy0[2], 20-s)
         || cbits0 || ints_to_interleavedbits( grid_br.xyS_collapseTo_ijL(x-xy0[1], y-xy0[2], S) )
  FROM (
  FROM (
   SELECT grid_br.ijS_to_xySref(ij0) as xy0, grid_br.IJ0_to_L0(ij0,false) as cbits0
   SELECT grid_br.ijS_to_xySref(ij0) as xy0, grid_br.IJ0_to_L0(ij0,false) as cbits0
Linha 142: Linha 178:
</syntaxhighlight>
</syntaxhighlight>


==Demais convenções==
===Demais convenções===
...
...


Linha 154: Linha 190:


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%;">
-- 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>


== Helper lib ==
== Helper lib ==
Peter
2 402

edições

Disponível em “https://wiki.addressforall.org/doc/Especial:MobileDiff/4051...4199”