Thursday 17 October 2019

Rotating a vector layer in QGIS with qgsAffine (or other method)?


I would like to rotate a set of vector points in QGIS an arbitrary number of degrees around a central point (or arbitrary point).


This is similar to a recent question about creating a regular grid; it was suggested there to use the "Affine Transformation" tool (which I assume meant the plugin) to rotate or shift a grid of points an arbitrary angle or distance. I suspect I'm not understanding how it works, and haven't been able to get it to work.


I create a regular grid of points in QGIS, and ensure that the UTM zone is set correctly for both layer and project, enable editing for the layer, then open the plugin dialog (qgsAffine): Affine Transformation dialog


I select 'whole layer' and then, wanting to rotate the entire field of points by 15°, put 15 in both 'rotation' boxes (which may be where things are going wrong). The operation results in rotating the points somewhere off-planet!


Is this the right tool for the job? I would like to rotate a set of points about their common centre, ideally.


Update: qgsAffine is just a thought; if we can do this in any QGIS tool I'll be happy!


Update 2: qgsAffine is usable IF you know the right numbers to plug in (see answer below, thanks Mike!). Spreadsheet / calculator works fine, or here's the R function to get the numbers directly:


## Compute correct affine numbers for qgsAffine plugin
affine <- function(originX, originY, rotAngle) {

A <- rotAngle * pi / 180
scaleX <- scaleY <- cos(A)
rotX <- sin(A)
rotY <- -sin(A)
transX <- originX - cos(A) * originX + sin(A) * originY
transY <- originY - sin(A) * originX - cos(A) * originY
aff <- data.frame(scaleX, scaleY, rotX, rotY, transX, transY)
return(aff)
}


So, to rotate a grid of points in northern Uganda (UTM 36N), affine(578988, 419210, 30) gives:


     scaleX    scaleY rotX rotY   transX    transY
1 0.8660254 0.8660254 0.5 -0.5 287174.7 -233330.5

... which, entered into the qgsAffine dialog, correctly rotates the points.



Answer



You can do this in PostGIS using ST_Affine. The functionality to rotate around an arbitrary point was added to ST_Rotate for PostGIS 2.0.


If you have an earlier version (like PostGIS 1.5, or even earlier), you can add these functions:


CREATE OR REPLACE FUNCTION st_rotate(geometry, double precision, geometry)
RETURNS geometry AS

'SELECT ST_Affine($1, cos($2), -sin($2), 0, sin($2), cos($2), 0, 0, 0, 1, ST_X($3) - cos($2) * ST_X($3) + sin($2) * ST_Y($3), ST_Y($3) - sin($2) * ST_X($3) - cos($2) * ST_Y($3), 0)'
LANGUAGE sql IMMUTABLE STRICT
COST 100;
COMMENT ON FUNCTION st_rotate(geometry, double precision, geometry) IS 'args: geomA, rotRadians, pointOrigin - Rotate a geometry rotRadians counter-clockwise about an origin.';

CREATE OR REPLACE FUNCTION st_rotate(geometry, double precision, double precision, double precision)
RETURNS geometry AS
'SELECT ST_Affine($1, cos($2), -sin($2), 0, sin($2), cos($2), 0, 0, 0, 1, $3 - cos($2) * $3 + sin($2) * $4, $4 - sin($2) * $3 - cos($2) * $4, 0)'
LANGUAGE sql IMMUTABLE STRICT
COST 100;

COMMENT ON FUNCTION st_rotate(geometry, double precision, double precision, double precision) IS 'args: geomA, rotRadians, x0, y0 - Rotate a geometry rotRadians counter-clockwise about an origin.';

See examples at ST_Rotate to get an idea on how to use it to rotate a geometry around an x, y point, including the centroid (common centre).


Because we all like math, the transformation matrix from the above functions is represented as:


[ cos(θ)  | -sin(θ)  ||  x0 - cos(θ) * x0 + sin(θ) * y0 ]
[ sin(θ) | cos(θ) || y0 - sin(θ) * x0 - cos(θ) * y0 ]

Where θ is the counter-clockwise rotation about an origin, x0 is the Easting/Longitude of the origin point, and y0 is the Northing/Latitude. This math could possibly be adapted to any affine transformation tool.




To use the qgsAffine tool, you need to understand where the values of the matrix flow to. A good spreadsheet template is also required to do pre-calculations. The qgsAffine dialog looks something like this:



              X   Y
+---+---+
Scale | a | e |
+---+---+
Rotation | d | b |
+---+---+
Translation | c | f |
+---+---+

where:




  • a : cos(θ)

  • b : -sin(θ)

  • c : x0 - cos(θ) * x0 + sin(θ) * y0

  • d : sin(θ)

  • e : cos(θ)

  • f : y0 - sin(θ) * x0 - cos(θ) * y0


For example, if you want to rotate a polygon 30° clockwise around 42°S, 174°E, here are your inputs to your spreadsheet:




  • x0 = 174

  • y0 = -42

  • θ=-30 degrees or -0.523598776 radians


Then, copy/paste the results from a spreadsheet to the right box. Using the tab order in the from the dialog:



  • a : 0.866025404

  • d : -0.5

  • c : 44.31157974

  • e : 0.866025404


  • b : 0.5

  • f : 81.37306696


qgsAffine


The same example from PostGIS would look something like:


SELECT ST_Rotate(geom, -30*pi()/180, 174.0, -42.0)

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...