Comment créer une grille de points régulière à l'intérieur d'un polygone dans Postgis?

Comment créer, à l'intérieur d'un polygone, une grille régulière de points espacés x, y en postgis? Comme l'exemple:

Vous faites cela avec generate_series.

Si vous ne voulez pas écrire manuellement où la grille doit commencer et s'arrêter, le plus simple est de créer une fonction.

Je n'ai pas testé correctement ce qui suit, mais je pense que cela devrait fonctionner:

CREATE OR REPLACE FUNCTION makegrid(geometry, integer)
RETURNS geometry AS
'SELECT ST_Collect(ST_POINT(x,y)) FROM 
generate_series(floor(st_xmin($1))::int, ceiling(st_xmax($1)-st_xmin($1))::int, $2) as x
,generate_series(floor(st_ymin($1))::int, ceiling(st_ymax($1)-st_ymin($1))::int,$2) as y 
where st_intersects($1,ST_POINT(x,y))'
LANGUAGE sql

Pour l'utiliser, vous pouvez faire:

SELECT makegrid(the_geom, 1000) from mytable;

où le premier argument est le polygone dans lequel vous souhaitez insérer la grille et le deuxième argument est la distance entre les points de la grille.

Si vous voulez un point par ligne, vous utilisez simplement ST_Dump comme:

SELECT (ST_Dump(makegrid(the_geom, 1000))).geom as the_geom from mytable;

HTH

Nicklas

J'ai récupéré le code de la fonction makegrid de Nicklas Avén et je l'ai rendu un peu plus générique en lisant et en utilisant le srid de la géométrie du polygone. Sinon, l'utilisation d'un polygone avec un srid défini donnerait une erreur.

La fonction:

CREATE OR REPLACE FUNCTION makegrid(geometry, integer)
RETURNS geometry AS
'SELECT ST_Collect(ST_SetSRID(ST_POINT(x,y),ST_SRID($1))) FROM 
generate_series(floor(st_xmin($1))::int, ceiling(st_xmax($1)-st_xmin($1))::int, $2) as x
,generate_series(floor(st_ymin($1))::int, ceiling(st_ymax($1)-st_ymin($1))::int,$2) as y 
where st_intersects($1,ST_SetSRID(ST_POINT(x,y),ST_SRID($1)))'
LANGUAGE sql

Pour utiliser la fonction se fait exactement comme Nicklas Avén l'a écrit:

SELECT makegrid(the_geom, 1000) from mytable;

ou si vous voulez un point par ligne:

SELECT (ST_Dump(makegrid(the_geom, 1000))).geom as the_geom from mytable;

J'espère que cela sera utile pour quelqu'un.

Alex

Les personnes utilisant une géométrie wgs84 auront probablement des problèmes avec cette fonction car

generate_series(floor(st_xmin($1))::int, ceiling(st_xmax($1))::int,$2) as x
,generate_series(floor(st_ymin($1))::int, ceiling(st_ymax($1))::int,$2) as y 

ne renvoient que des entiers. Sauf pour les très grandes géométries telles que les pays (qui reposent sur plusieurs degrés lat, lng), cela ne fera que collecter 1 point qui la plupart du temps ne coupe même pas la géométrie elle-même ... => résultat vide!

Mon problème était que je n'arrive pas à utiliser generate_series () avec une distance décimale sur des nombres flottants tels que ceux WSG84 ... C'est pourquoi j'ai modifié la fonction pour la faire fonctionner de toute façon:

SELECT ST_Collect(st_setsrid(ST_POINT(x/1000000::float,y/1000000::float),st_srid($1))) FROM 
  generate_series(floor(st_xmin($1)*1000000)::int, ceiling(st_xmax($1)*1000000)::int,$2) as x ,
  generate_series(floor(st_ymin($1)*1000000)::int, ceiling(st_ymax($1)*1000000)::int,$2) as y 
WHERE st_intersects($1,ST_SetSRID(ST_POINT(x/1000000::float,y/1000000::float),ST_SRID($1)))

Fondamentalement, exactement la même chose. Il suffit de multiplier et de diviser par 1000000 pour obtenir les décimales dans le jeu quand j'en ai besoin.

Il y a sûrement une meilleure solution pour y parvenir. ++

Cet algorithme devrait être bien:

createGridInPolygon(polygon, resolution) {
    for(x=polygon.xmin; x<polygon.xmax; x+=resolution) {
       for(y=polygon.ymin; y<polygon.ymax; y+=resolution) {
          if(polygon.contains(x,y)) createPoint(x,y);
       }
    }
}

où «polygone» est le polygone et «résolution» est la résolution de grille requise.

Pour l'implémenter dans PostGIS, les fonctions suivantes peuvent être nécessaires:

• ST_XMin , ST_XMax , ST_YMin et ST_YMax pour obtenir les coordonnées minimales et maximales du polygone,
• ST_Contient pour tester si le polygone contient le point,
• et ST_Point pour créer un point.

Bonne chance!

Trois algorithmes utilisant différentes méthodes.

Lien de dépôt Github

1. L'approche simple et la meilleure, en utilisant la distance réelle de la terre des coordonnées de la direction x et y. L'algorithme fonctionne avec n'importe quel SRID, en interne, il fonctionne avec WGS 1984 (EPSG: 4326) et la transformation de résultat en entrée SRID.

Fonction ================================================= ==================

CREATE OR REPLACE FUNCTION public.I_Grid_Point_Distance(geom public.geometry, x_side decimal, y_side decimal)
RETURNS public.geometry AS $BODY$
DECLARE
x_min decimal;
x_max decimal;
y_max decimal;
x decimal;
y decimal;
returnGeom public.geometry[];
i integer := -1;
srid integer := 4326;
input_srid integer;
BEGIN
CASE st_srid(geom) WHEN 0 THEN
    geom := ST_SetSRID(geom, srid);
        ----RAISE NOTICE 'No SRID Found.';
    ELSE
        ----RAISE NOTICE 'SRID Found.';
END CASE;
    input_srid:=st_srid(geom);
    geom := st_transform(geom, srid);
    x_min := ST_XMin(geom);
    x_max := ST_XMax(geom);
    y_max := ST_YMax(geom);
    y := ST_YMin(geom);
    x := x_min;
    i := i + 1;
    returnGeom[i] := st_setsrid(ST_MakePoint(x, y), srid);
<<yloop>>
LOOP
IF (y > y_max) THEN
    EXIT;
END IF;

CASE i WHEN 0 THEN 
    y := ST_Y(returnGeom[0]);
ELSE 
    y := ST_Y(ST_Project(st_setsrid(ST_MakePoint(x, y), srid), y_side, radians(0))::geometry);
END CASE;

x := x_min;
<<xloop>>
LOOP
  IF (x > x_max) THEN
      EXIT;
  END IF;
    i := i + 1;
    returnGeom[i] := st_setsrid(ST_MakePoint(x, y), srid);
    x := ST_X(ST_Project(st_setsrid(ST_MakePoint(x, y), srid), x_side, radians(90))::geometry);
END LOOP xloop;
END LOOP yloop;
RETURN
ST_CollectionExtract(st_transform(ST_Intersection(st_collect(returnGeom), geom), input_srid), 1);
END;
$BODY$ LANGUAGE plpgsql IMMUTABLE;

Utilisez la fonction avec une requête simple, la géométrie doit être valide et de type polygone, multi-polygones ou enveloppe

SELECT I_Grid_Point_Distance(geom, 50, 61) from polygons limit 1;

Résultat ================================================= =====================

2. Deuxième fonction basée sur l' algorithme de Nicklas Avén . Je l'ai amélioré pour gérer n'importe quel SRID.

   avoir mis à niveau les modifications suivantes dans l'algorithme.

   1. Variable séparée pour les directions x et y pour la taille des pixels,
   2. Nouvelle variable pour calculer la distance en sphéroïde ou ellipsoïde.
   3. Saisissez n'importe quel SRID, transformez la fonction Geom dans l'environnement de travail de Spheroid ou Ellipsoid Datum, puis appliquez la distance de chaque côté, obtenez le résultat et transformez-le en SRID d'entrée.

Fonction ================================================= ==================

CREATE OR REPLACE FUNCTION I_Grid_Point(geom geometry, x_side decimal, y_side decimal, spheroid boolean default false)
RETURNS SETOF geometry AS $BODY$ 
DECLARE
x_max decimal; 
y_max decimal;
x_min decimal;
y_min decimal;
srid integer := 4326;
input_srid integer; 
BEGIN
CASE st_srid(geom) WHEN 0 THEN
  geom := ST_SetSRID(geom, srid);
  RAISE NOTICE 'SRID Not Found.';
    ELSE
        RAISE NOTICE 'SRID Found.';
    END CASE;

    CASE spheroid WHEN false THEN
        RAISE NOTICE 'Spheroid False';
        srid := 4326;
        x_side := x_side / 100000;
        y_side := y_side / 100000;
    else
        srid := 900913;
        RAISE NOTICE 'Spheroid True';
    END CASE;
    input_srid:=st_srid(geom);
    geom := st_transform(geom, srid);
    x_max := ST_XMax(geom);
    y_max := ST_YMax(geom);
    x_min := ST_XMin(geom);
    y_min := ST_YMin(geom);
RETURN QUERY
WITH res as (SELECT ST_SetSRID(ST_MakePoint(x, y), srid) point FROM
generate_series(x_min, x_max, x_side) as x,
generate_series(y_min, y_max, y_side) as y
WHERE st_intersects(geom, ST_SetSRID(ST_MakePoint(x, y), srid))
) select ST_TRANSFORM(ST_COLLECT(point), input_srid) from res;
END;
$BODY$ LANGUAGE plpgsql IMMUTABLE STRICT;

Utilisez-le avec une simple requête.

SELECT I_Grid_Point(geom, 22, 15, false) from polygons;

Résultat ================================================= ==================

3. Fonction basée sur un générateur en série.

Fonction ================================================= =================

CREATE OR REPLACE FUNCTION I_Grid_Point_Series(geom geometry, x_side decimal, y_side decimal, spheroid boolean default false)
RETURNS SETOF geometry AS $BODY$
DECLARE
x_max decimal;
y_max decimal;
x_min decimal;
y_min decimal;
srid integer := 4326;
input_srid integer;
x_series DECIMAL;
y_series DECIMAL;
BEGIN
CASE st_srid(geom) WHEN 0 THEN
  geom := ST_SetSRID(geom, srid);
  RAISE NOTICE 'SRID Not Found.';
    ELSE
        RAISE NOTICE 'SRID Found.';
    END CASE;

    CASE spheroid WHEN false THEN
        RAISE NOTICE 'Spheroid False';
    else
        srid := 900913;
        RAISE NOTICE 'Spheroid True';
    END CASE;
    input_srid:=st_srid(geom);
    geom := st_transform(geom, srid);
    x_max := ST_XMax(geom);
    y_max := ST_YMax(geom);
    x_min := ST_XMin(geom);
    y_min := ST_YMin(geom);

    x_series := CEIL ( @( x_max - x_min ) / x_side);
    y_series := CEIL ( @( y_max - y_min ) / y_side );
RETURN QUERY
SELECT st_collect(st_setsrid(ST_MakePoint(x * x_side + x_min, y*y_side + y_min), srid)) FROM
generate_series(0, x_series) as x,
generate_series(0, y_series) as y
WHERE st_intersects(st_setsrid(ST_MakePoint(x*x_side + x_min, y*y_side + y_min), srid), geom);
END;
$BODY$ LANGUAGE plpgsql IMMUTABLE STRICT;

Utilisez-le avec une simple requête.

SELECT I_Grid_Point_Series(geom, 22, 15, false) from polygons; Résultat ================================================= =========================

Donc ma version fixe:

CREATE OR REPLACE FUNCTION makegrid(geometry, integer, integer)
RETURNS geometry AS
'SELECT ST_Collect(st_setsrid(ST_POINT(x,y),$3)) FROM 
  generate_series(floor(st_xmin($1))::int, ceiling(st_xmax($1))::int,$2) as x
  ,generate_series(floor(st_ymin($1))::int, ceiling(st_ymax($1))::int,$2) as y 
where st_intersects($1,st_setsrid(ST_POINT(x,y),$3))'
LANGUAGE sql

Usage:

SELECT (ST_Dump(makegrid(the_geom, 1000, 3857))).geom as the_geom from my_polygon_table

Voici une autre approche qui est certainement plus rapide et plus facile à comprendre.

Par exemple pour une grille de 1000 m sur 1000 m:

SELECT (ST_PixelAsCentroids(ST_AsRaster(the_geom,1000.0,1000.0))).geom 
FROM the_polygon

Le SRID d'origine est également conservé.

Cet extrait de code convertit la géométrie du polygone en un raster vide, puis convertit chaque pixel en un point. Avantage: il n'est pas nécessaire de vérifier à nouveau si le polygone d'origine coupe les points.

Optionnel:

Vous pouvez également ajouter l'alignement de la grille avec les paramètres gridx et gridy. Mais puisque nous utilisons le centroïde de chaque pixel (et non un coin), nous devons utiliser un modulo pour calculer la bonne valeur:

SELECT (ST_PixelAsCentroids(ST_AsRaster(the_geom,1000.0,1000.0,mod(1000/2,100),mod(1000/2,100)))).geom 
FROM the_polygon

Avec mod(grid_size::integer/2,grid_precision)

Voici la fonction postgres:

CREATE OR REPLACE FUNCTION st_makegrid(geometry, float, integer)
RETURNS SETOF geometry AS
'SELECT (ST_PixelAsCentroids(ST_AsRaster($1,$2::float,$2::float,mod($2::int/2,$3),mod($2::int/2,$3)))).geom'
LANGUAGE sql;

Peut être utilisé avec:

SELECT makegrid(the_geom,1000.0,100) as geom from the_polygon  
-- makegrid(the_geom,grid_size,alignement)

Une petite mise à jour potentielle des réponses précédentes - troisième argument comme échelle pour wgs84 (ou utiliser 1 pour les normales), et également arrondi à l'intérieur du code afin que les points mis à l'échelle sur plusieurs formes soient alignés.

J'espère que cela vous aide, Martin

CREATE OR REPLACE FUNCTION makegrid(geometry, integer, integer)
RETURNS geometry AS



/*geometry column , integer: distance between points, integer: scale factor for distance (useful for wgs84, e.g. use there 50000 as distance and 1000000 as scale factor*/

'
SELECT ST_Collect(st_setsrid(ST_POINT(x/$3::float,y/$3::float),st_srid($1))) FROM 
  generate_series(
                (round(floor(st_xmin($1)*$3)::int/$2)*$2)::int, 
                (round(ceiling(st_xmax($1)*$3)::int/$2)*$2)::int,
                $2) as x ,
  generate_series(
                (round(floor(st_ymin($1)*$3)::int/$2)*$2)::int, 
                (round(ceiling(st_ymax($1)*$3)::int/$2)*$2)::int,
                $2) as y 
WHERE st_intersects($1,ST_SetSRID(ST_POINT(x/$3::float,y/$3::float),ST_SRID($1)))
'

LANGUAGE sql