Message posté par : Sylvain M.
----------------------------------------
C'est bon, j'ai compris le problème, et je l'ai résolu.
Il venait de la fonction hexagoncoordinates().
Comme présenté en commentaire, celle-ci :
-----------------
Citation :
Étant donné un carré, trouve toutes les cellules hexagonales d'un pavage hexagonal
(déterminé par la taille des côtés) qui pourraient recouvrir ce carré (légèrement
surdéterminé).
-----------------
Et elle n'était donc pas adaptée pour prendre en entrée une géométrie de point.
J'ai donc créé une nouvelle fonction qui répond à mes besoins (j'avoue, je me suis
fait aider par une iA) :
-----------------
Code :
CREATE OR REPLACE FUNCTION hexagoncoordinates_frompoint(
pt geometry,
edge float8,
OUT i integer,
OUT j integer
)
RETURNS record
AS $$
DECLARE
h float8 := edge * cos(pi()/6);
x float8 := ST_X(pt);
y float8 := ST_Y(pt);
srid integer := ST_SRID(pt);
i0 integer := floor(x / (1.5*edge));
j0 integer;
ci integer;
cj integer;
cx float8;
cy float8;
hex geometry;
BEGIN
-- estimation brute du j en tenant compte du décalage éventuel
IF (i0 % 2) = 0 THEN
j0 := floor((y) / (2*h));
ELSE
j0 := floor((y - h) / (2*h));
END IF;
-- tester la cellule et les voisines
FOR ci, cj IN
SELECT i0 + dx, j0 + dy
FROM (VALUES
(0,0),(1,0),(-1,0),(0,1),(0,-1),(1,-1),(-1,1),
(1,1),(-1,-1) -- un peu plus sûr
) AS v(dx,dy)
LOOP
-- centre de la cellule (ci,cj)
cx := (1.5 * edge) * ci;
cy := (2 * h) * cj + (CASE WHEN (ci % 2) = 0 THEN 0 ELSE h END);
-- polygone hexagonal autour du centre
hex := ST_SetSRID(
ST_MakePolygon(
ST_MakeLine(ARRAY[
ST_MakePoint(cx + edge*cos(radians(0)), cy +
edge*sin(radians(0))),
ST_MakePoint(cx + edge*cos(radians(60)), cy +
edge*sin(radians(60))),
ST_MakePoint(cx + edge*cos(radians(120)), cy +
edge*sin(radians(120))),
ST_MakePoint(cx + edge*cos(radians(180)), cy +
edge*sin(radians(180))),
ST_MakePoint(cx + edge*cos(radians(240)), cy +
edge*sin(radians(240))),
ST_MakePoint(cx + edge*cos(radians(300)), cy +
edge*sin(radians(300))),
ST_MakePoint(cx + edge*cos(radians(0)), cy + edge*sin(radians(0)))
])
),
srid
);
IF ST_Contains(hex, pt) THEN
i := ci;
j := cj;
RETURN;
END IF;
END LOOP;
RAISE EXCEPTION
'Le point % ne correspond à aucun hexagone pour edge % (vérifiez la
grille)',
ST_AsText(pt), edge;
END;
$$ LANGUAGE plpgsql IMMUTABLE STRICT PARALLEL SAFE;
-----------------
Et ça marche très bien !
Et pour ceux que ça intéresse, voici donc une requête qui créé, à partir d'une table
d'observations (ponctuelles), une table de mailles hexagonales de 0.01 degrés, avec
des statistiques (nb_obs, et nb_taxons) pour chaque maille :
-----------------
Code :
CREATE TABLE obs_hex_grid_001 AS (
SELECT ROW_NUMBER() OVER () AS id,
count(id) AS nb_obs,
count(distinct taxon_id) AS nb_taxons,
ST_SetSRID(hexagon((hexagoncoordinates_frompoint(geom, 0.01)).i,
(hexagoncoordinates_frompoint(geom, 0.01)).j, 0.01),4326) AS hex_geom
FROM observations
GROUP BY hex_geom);
-----------------
----------------------------------------
Le message est situé
https://georezo.net/forum/viewtopic.php?pid=375315#p375315
Pour vous désabonner connectez-vous sur le forum puis Profil / Abonnement
--
Association GeoRezo - le portail géomatique
https://georezo.net