Le portail francophone de la géomatique

Sylvain M. · Sat 15 November 2025 00:59

Bonsoir à tous,

J'ai plusieurs projets persos où je souhaite faire une synthèse par mailles hexagonales de jeux de données ponctuels volumineux (~200 000 entités) et d'emprise mondiale (pour ceux qui sont curieux : photothèque ou observations naturalistes, mais le principe est le même).
Voici un exemple de carte interactive utilisant ce principe de mailles hexagonales : https://api.gbif.org/v2/map/demo2.html

Plutôt que de générer des couches de mailles mondiales à chaque échelle, ce qui seraient extrêmement volumineux vu que je souhaite descendre jusqu'à des mailles de quelques centaines de mètres, j'essaie de générer les géométries de mailles à la volée depuis les points source.

J'ai trouvé pour cela les fonctions hexagon() et hexagoncoordinates() sur cette page github :
https://github.com/CrunchyData/pg_tiles … in-example

J'essaie de les adapter/comprendre, mais je constate un comportement étrange sur mes données : pour un point en entrée, les fonctions retournent plusieurs mailles hexagonales.

Voici un exemple qui illustre mon problème (1 seul point, et 1 maille de 0.1 degrés) :

Code:

WITH    point_test AS (SELECT st_geomfromtext('POINT (0.47003 47.65044)',4326) AS geom)
SELECT    geom,
        hexagoncoordinates(geom,0.1) AS hexagon_coordinates,
        ST_SetSRID(hexagon((hexagoncoordinates(geom, 0.1)).i, (hexagoncoordinates(geom, 0.1)).j, 0.1),4326) AS hexagon_geom
FROM    point_test

Cette requête depuis un point fictif ne devrait sortir qu'un polygone, mais il en sort 4 :

Code:

geom                    hexagon_coordinates
POINT (0.47003 47.65044)    (3,275)
POINT (0.47003 47.65044)    (3,276)
POINT (0.47003 47.65044)    (4,275)
POINT (0.47003 47.65044)    (4,276)

Si jamais un expert PostGis passe par là, et si le challenge l'intéresse, je serai très reconnaissant s'il me permet d'avancer !

Dernière modification par Sylvain M. (Sat 15 November 2025 01:12)

Sylvain M. · Hier 13:53

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 :

É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);

Communauté	Ressources	L'emploi	Le marché	L'association	Fils RSS
Forum Blogs Wiki Planet Sites amis	Agenda Biblio Fiches AMO-CNIG Paris PCGIS Téléchargements	Offres d'emploi Banque de CV Formations Démarche Métiers	Geo-Entreprises Geo-Communiqué Appels d'offres	L'association Qui sommes-nous ? Nous aider Bienfaiteurs Partenaires Historique Statuts	Tout le site Job Géomatique Services GeoBlogs ...	Faire un don Mentions légales

Le portail francophone de la géomatique

Annonce

#1 Sat 15 November 2025 00:59

[PostGis] Grille Hexagonale d'analyse

Code:

Code:

#2 Hier 13:53

Re: [PostGis] Grille Hexagonale d'analyse

Code:

Code:

Pied de page des forums