3

I have a set of data in a Postgres database. Originally the data was point data that I converted to polygons with makegrid_2d from this post. Now I want to be able to rotate the polygons (around each polygons centroid). I tried to do it with several methods hence the result is the same as if I just rotate the features (separately) in QGIS like the ones down to the left:

First try to rotate

I end up with new polygons overlapping and not covering the entire area.

Any suggestions on how to improve it?

EDIT1: I release that I wasn't clear on what I wanted, I want to skew the polygons. Meaning that I want to change the polygons instead of having borders pointing (0-180deg resp 90-270deg) I want to "tilt" the polygon so its borders are pointing 20-200 resp 110-290deg. Does it make any sense?

I used st_astext(ST_Rotate(polygon,radians(20),ST_Centroid(ST_GeomFromText(ST_AsText(polygon))))) to rotate the polygons in the image.

  • The rotation must use the center of the whole grid! – Andreas Müller Jun 26 '18 at 21:30
  • wait...those nine polygons seem to be quite neatly rotated? what ecactly do you want? did you forget a 'not' somewhere ,) ? and what methods didn't work how? – ThingumaBob Jun 26 '18 at 21:35
  • What do u mean with that the rotation must use the center of the whole grid @AndreasMüller? – axel_ande Jun 27 '18 at 4:25
  • Can you share the function you used and the steps you took to get to the current state. – Cushen Jun 27 '18 at 6:38
3

As the comments suggested, you need to rotate around the point of origin, which in your case is the centroid of all the polygons. Essentially, you need to collect all the polygon geometries, define their centroid, and rotate around that. The following works for me (including the code I used to generate a basic grid with makegrid_2d):

--Defines the grid (1km rectangles inside the bounds)
WITH grid AS (
  SELECT 
    ROW_NUMBER() OVER () AS grid_id,
    m.geom 
  FROM (
    SELECT (
      ST_Dump(
        MAKEGRID_2D(
          ST_GeomFromText('POLYGON((-72.0268017 46.5390344,-72.0268017 47.0837218,-70.8129371 47.0837218,-70.8129371 46.5390344,-72.0268017 46.5390344))',4326), 
          1000, 
          32187
        )
      )
    ).geom
  ) m
),
--Defines the centroid of the whole grid
centroid AS (
  SELECT ST_Centroid(ST_Collect(grid.geom)) AS geometry FROM grid
)
--Rotates around the defined centroid
SELECT ST_Rotate(grid.geom,radians(20),(SELECT geometry FROM centroid)) FROM grid

Result :

enter image description here

1

Looking at the answer you are referring to, try editing these portions of code.

NextX := ST_X(ST_Project(CPoint, $2, radians(90))::geometry);
NextY := ST_Y(ST_Project(CPoint, $3, radians(0))::geometry);

to be:
(instead of 0,90 values are changed to 20, 110)

NextX := ST_X(ST_Project(CPoint, $2, radians(110))::geometry);
NextY := ST_Y(ST_Project(CPoint, $3, radians(20))::geometry);

The full function, when edited, would be:

CREATE OR REPLACE FUNCTION public.makegrid_2d (
  bound_polygon public.geometry,
  width_step integer,
  height_step integer
)
RETURNS public.geometry AS
$body$
DECLARE
  Xmin DOUBLE PRECISION;
  Xmax DOUBLE PRECISION;
  Ymax DOUBLE PRECISION;
  X DOUBLE PRECISION;
  Y DOUBLE PRECISION;
  NextX DOUBLE PRECISION;
  NextY DOUBLE PRECISION;
  CPoint public.geometry;
  sectors public.geometry[];
  i INTEGER;
  SRID INTEGER;
BEGIN
  Xmin := ST_XMin(bound_polygon);
  Xmax := ST_XMax(bound_polygon);
  Ymax := ST_YMax(bound_polygon);
  SRID := ST_SRID(bound_polygon);

  Y := ST_YMin(bound_polygon); --current sector's corner coordinate
  i := -1;
  <<yloop>>
  LOOP
    IF (Y > Ymax) THEN  
        EXIT;
    END IF;

    X := Xmin;
    <<xloop>>
    LOOP
      IF (X > Xmax) THEN
          EXIT;
      END IF;

      CPoint := ST_SetSRID(ST_MakePoint(X, Y), SRID);
      --orientation of the grid cells:
      NextX := ST_X(ST_Project(CPoint, $2, radians(110))::geometry);
      NextY := ST_Y(ST_Project(CPoint, $3, radians(20))::geometry);

      i := i + 1;
      sectors[i] := ST_MakeEnvelope(X, Y, NextX, NextY, SRID);

      X := NextX;
    END LOOP xloop;
    CPoint := ST_SetSRID(ST_MakePoint(X, Y), SRID);
    NextY := ST_Y(ST_Project(CPoint, $3, radians(20))::geometry);
    Y := NextY;
  END LOOP yloop;

  RETURN ST_Collect(sectors);
END;
$body$
LANGUAGE 'plpgsql';
0

Thanks for the help! The "problem" with the solution of changing the make_grid function is that I want to have the rotation angle as a variable and the "problem" with with just rotating the data is that it wont cover the original spot.

So, inspired by Jufaua's solution, I ended up with a few more with blocks creating a "too large" grid, then rotating that, checking which of the rotated data that covered the original data and finally "joining in" other parameters. The two first steps: (from green to purple)Step1

And the last two steps:(From purple (step1) to orange to blue)Step2

And here is the code (with pos as geometry column and field_info as a double precision column in "schema.table"):

    WITH grid AS (
      SELECT 
        ROW_NUMBER() OVER () AS grid_id,
        m.geom 
      FROM (
        SELECT (
          ST_Dump(
            MAKEGRID_2D(
              ST_SetSRID(st_buffer(ST_Extent(pos),
                                                       GREATEST(((select max(st_x(pos)) from schema.table) - 
                                                            (select min(st_x(pos)) from schema.table)),
                                                           ((select max(st_y(pos)) from schema.table) - 
                                                            (select min(st_y(pos)) from schema.table)))/4
                                                      ),4326),50,50))
             ).geom  from schema.table
      ) m
    ),
    --Defines the centroid of the whole grid
    centroid AS (
      SELECT ST_Centroid(ST_Collect(grid.geom)) AS geometry FROM grid
    ), 
    --Rotates around the defined centroid
    rotated as(SELECT ST_Rotate(grid.geom,radians(20),(SELECT geometry FROM centroid)) as polys 
               FROM grid
              ),
    --Selectes the polygons that are intersecting the orignal data
    select_data as (select polys
                   from rotated
                   where st_intersects(ST_SetSRID((select ST_Extent(pos) from schema.table), 4326),
                                       polys))
    --Do the final selections and joining in some average data
    select row_number () over() as r_num, 
    polys,
    avg(2*field_info) 
    from select_data
    left join schema.table on st_intersects(polys, pos)
    group by polys;

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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