1

I am looking to create an irregular grid of points via PostGIS, using the four provided corner (point) coordinates, and number of point-rows, and point-columns.

Four corners may or may not be perpendicular to each other (see image, 7 rows, 5 columns). Terrain is considered flat. What might be the easiest approach to solve this, within PostGIS ?

enter image description here

6
  • 4
    Can you provide an example image/sketch of what you're trying to achieve?
    – Erik
    Aug 8, 2022 at 11:38
  • 1
    Automation of regular grids is pretty basic, bur, by definition, an irregular grid would seem automation-resistant.
    – Vince
    Aug 8, 2022 at 11:42
  • Please include relevant images in your question.
    – Erik
    Aug 8, 2022 at 13:43
  • My approach was first creating four sides (ST_InterpolatePoints), then filling the inside, row by row. But not sure if this can be applicable in SQL/PostGIS limits.
    – Alper ALT
    Aug 9, 2022 at 6:34
  • It is still not clear from your question what you would like to have in the output, for example: 1) only geometry from points as in the picture, or 2) geometry from points and numbered columns and rows calculated from one corner (1 or 2 or 3 or 4)...This is important... Aug 10, 2022 at 19:31

2 Answers 2

1

So, there are many ways to solve your question and this approach is one of them.

Create a fun function called ST_RegularPointsGridOfCornerPoints

DROP FUNCTION ST_RegularPointsGridOfCornerPoints

CREATE OR REPLACE FUNCTION ST_RegularPointsGridOfCornerPoints(
    geom GEOMETRY,
    r bigint,
    c bigint)
RETURNS GEOMETRY AS  
$BODY$
WITH 
    tbla AS (SELECT ST_Boundary(ST_Union(geom)) geom FROM (SELECT ((ST_DelaunayTriangles(ST_Collect(geom)))) geom) foo),            
    tblb AS (SELECT row_number() over() AS id,
             ST_MakeLine(pt1, pt2) geom FROM (SELECT ST_PointN(geom, generate_series(1, ST_NPoints(geom)-1)) pt1,
     ST_PointN(geom, generate_series(2, ST_NPoints(geom))) pt2 FROM tbla) AS geom),
    tblc AS (SELECT generate_series (0,r-1) as steps),
    tbld AS (SELECT steps AS stp1, ST_LineInterpolatePoint(geom, steps/(SELECT count(steps)::float-1 FROM tblc)) geom1 FROM tblc, tblb WHERE tblb.id 
IN (2) GROUP BY tblc.steps, geom),
    tble AS (SELECT steps AS stp2, ST_LineInterpolatePoint(ST_Reverse(geom), steps/(SELECT count(steps)::float-1 FROM tblc)) geom2 FROM tblc, tblb 
WHERE tblb.id IN (4) GROUP BY tblc.steps, geom),
    tblf AS (SELECT row_number() over() AS id, ST_MakeLine(geom1, geom2) geom FROM tbld JOIN tble ON true AND stp1=stp2),
    tblg AS (SELECT generate_series (0,c-1) as steps)
      (SELECT ST_LineInterpolatePoint(geom, steps/(SELECT count(steps-1)::float-1 FROM tblg)) geom FROM tblg, tblf geom);
$BODY$
LANGUAGE SQL

Run

SELECT ST_RegularPointsGridOfCornerPoints(ST_Union(geom), 7, 5) geom FROM <name_table>

See the result - Unfortunately something went wrong and it works not on all versions of PostgreSQL builds (For example, for PostgreSQL 14.0, compiled by Visual C++ build 1914, 64-bit and higher should work :-))... Remember my comment, its future hasn't come yet :-(...

As a consequence, for now, run the body of the function as a CTE and set the required values of columns and rows, for example, as specified in your question for your example. The architecture of the SQL-code is shown below:

create table <name_table> AS
WITH 
    tbla AS (SELECT ST_Boundary(ST_Union(geom)) geom FROM (SELECT ((ST_DelaunayTriangles(ST_Collect(geom)))) geom FROM layer_1) foo),           
    tblb AS (SELECT row_number() over() AS id,
             ST_MakeLine(pt1, pt2) geom FROM (SELECT ST_PointN(geom, generate_series(1, ST_NPoints(geom)-1)) pt1,
     ST_PointN(geom, generate_series(2, ST_NPoints(geom))) pt2 FROM tbla) AS geom),
    tblc AS (SELECT generate_series (0,4) as steps),
    tbld AS (SELECT steps AS stp1, ST_LineInterpolatePoint(geom, steps/(SELECT count(steps)::float-1 FROM tblc)) geom1 FROM tblc, tblb WHERE tblb.id 
IN (2) GROUP BY tblc.steps, geom),
    tble AS (SELECT steps AS stp2, ST_LineInterpolatePoint(ST_Reverse(geom), steps/(SELECT count(steps)::float-1 FROM tblc)) geom2 FROM tblc, tblb 
WHERE tblb.id IN (4) GROUP BY tblc.steps, geom),
    tblf AS (SELECT row_number() over() AS id, ST_MakeLine(geom1, geom2) geom FROM tbld JOIN tble ON true AND stp1=stp2),
    tblg AS (SELECT generate_series (0,6) as steps)
      (SELECT ST_LineInterpolatePoint(geom, steps/(SELECT count(steps-1)::float-1 FROM tblg)) geom FROM tblg, tblf);

The figure below shows the result, you should get the same one for yourself...

enter image description here

The figure

Unfortunately, I only fancy fun and customizable functions and they are not always simple :-(...

P.S. In the following questions, try to present the SQL-code and an explanation of what prevented you from getting the expected result...

Original GeoSpatial Solutions...

Translated with www.DeepL.com/Translator (free version)

2
  • 1
    This is working exceptionally well, thank you very much !
    – Alper ALT
    Aug 13, 2022 at 11:17
  • I'm glad I helped you and others and thank you for the question... 8-)... Aug 15, 2022 at 17:24
0

An elegant way to do this is to transform a square grid of points into the required quadrilateral. Since the transformation does not preserve parallel lines, a projective transformation is normally required. This is complex to derive. But this answer describes a clever and simple transformation of the unit square to an arbitrary quadrilateral. The diagram below shows how it uses a combination of three vectors derived from the vertices of the quadrilateral.

Quadrilateral
transformation

Here's SQL implementing the transformation, with an example quadrilateral:

WITH quad AS (SELECT 
  5 AS LLx, 5 AS LLy,
  10 AS ULx, 30 AS ULy,
  25 AS URx, 25 AS URy,
  30 AS LRx, 0 AS LRy
),
vec AS (
  SELECT  LLx AS ox, LLy AS oy,
          ULx - LLx AS ux, ULy - LLy AS uy, 
          LRx - LLx AS vx, LRy - LLy As vy, 
          URx - LLx - ((ULX - LLx) + (LRx - LLx)) AS wx, 
          URy - LLy - ((ULy - LLy) + (LRy - LLy)) AS wy 
  FROM quad
),
grid AS (SELECT x / 10.0 AS x, y / 10.0 AS y 
  FROM        generate_series(0, 10) AS sx(x)
  CROSS JOIN  generate_series(0, 10) AS sy(y)
)
SELECT ST_Point(ox + ux * x + vx * y + wx * x * y, 
                oy + uy * x + vy * y + wy * x * y) AS geom
  FROM vec CROSS JOIN grid;

The result is:

enter image description here

1
  • 1
    Thank you for this nice solution, it appears both answers are solving my problem perfectly and equally well.
    – Alper ALT
    Sep 10, 2022 at 6:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.