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

31

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:

texte alternatif

Pablo
la source
J'ai essayé de créer des polygones de découpage en fusionnant ce code avec le code de découpage "postGIS en action" mais un seul polygone est créé ... Ai-je oublié quelque chose? CRÉER OU REMPLACER LA FONCTION makegrid (géométrie, entier, entier) RETOURNE la géométrie AS 'SELECT st_intersection (g1.geom1, g2.geom2) AS geom FROM (SELECT $ 1 AS geom1) AS g1 INNER JOIN (Sélectionnez st_setsrid (CAST (ST_MakeBox2d (st_setsrid (( ST_Point (x, y), $ 3), st_setsrid (ST_Point (x + $ 2, y + $ 2), $ 3)) comme géométrie), $ 3) comme geom2 FROM generate_series (étage (st_xmin ($ 1)) :: int, plafond ( st_xmax ($ 1)) :: int, $ 2) as x, generate_series (floor (st_ymin ($ 1)) :: int, plafond (st_ymax (
aurel_nc
regardez ma réponse détaillée.
Muhammad Imran Siddique

Réponses:

30

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

Nicklas Avén
la source
2
Vous devrez peut-être ajouter st_setSRID () aux fonctions st_point, sinon st_intersects ne fonctionnera pas.
JaakL
Ajout de ma version testée comme réponse séparée.
JaakL
12

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

Alexandre Neto
la source
La réponse acceptée ne fonctionne pas avec mes données en raison d'erreurs SRID. Cette modification fonctionne mieux.
Vitaly Isaev
Vous voudrez peut-être ajouter quelque chose lorsque le polygone traverse l'antiméridien? Je peux imaginer que cela causerait un problème avec xmin / xmax.
Thomas
2
Ça n'a pas marché pour moi. Utilisation de Postgres 9.6 et PostGIS 2.3.3. À l'intérieur de l'appel de generate_series, j'ai dû mettre ceci comme deuxième paramètre "plafond (st_xmax ($ 1)) :: int" au lieu de "plafond (st_xmax ($ 1) -st_xmin ($ 1)) :: int" et "plafond ( st_ymax ($ 1)) :: int "au lieu de" plafond (st_ymax ($ 1) -st_ymin ($ 1)) :: int "
Vitor Sapucaia
J'approuve le commentaire précédent; la limite supérieure de generate_series doit être le plafond du max, et non le plafond de la différence (max - min).
R. Bourgeon
10

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. ++

Julien Garcia
la source
C'est une solution de contournement intelligente. Avez-vous vérifié les résultats? Sont-ils cohérents?
Pablo
Salut. Oui pablo. Je suis satisfait des résultats jusqu'à présent. J'en avais besoin pour construire un polygone avec une altitude relative au-dessus du sol. (J'utilise SRTM pour calculer l'altitude que je veux pour chaque point de grille). Il me manque maintenant un moyen d'inclure également les points qui se trouvent sur le périmètre du polygone. Actuellement, la forme rendue est quelque peu tronquée au bord.
Julien Garcia
travaillé, quand toutes les autres solutions ont échoué, merci!
Jordan Arseno
7

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:

Bonne chance!

julien
la source
1
Notez que si vous avez de grandes zones de polygones complexes (par exemple, j'ai un tampon de côte), cette approche n'est pas tout à fait optimale.
JaakL
Alors, que proposez-vous à la place?
julien
4

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 ================================================= =====================

entrez la description de l'image ici

  1. 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 ================================================= ==================entrez la description de l'image ici

  1. 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 ================================================= =========================

entrez la description de l'image ici

Muhammad Imran Siddique
la source
3

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
JaakL
la source
1
Salut, j'obtiens un résultat vide avec la fonction makegrid. Le fichier de formes a été importé dans PostGIS à l'aide de shp2pgsql. Vous n'avez aucune idée de ce qui pourrait causer des problèmes, le srs est réglé sur wgs84.
Michal Zimmermann
3

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)
obchardon
la source
1

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
Martin
la source
Transformer la géométrie en un SRID spécifique (comme EPSG: 3857) ne serait-il pas préférable de simplement multiplier par un facteur d'échelle?
Nikolaus Krismer