3

I have a polygon dataset containing hundreds of rectangle map 'frames' that represent map extents across the CONUS. The dataset containing the frames is in EPSG:3857, but the frames themselves represent printed map extents in various UTM zones across the country. Each polygon 'frame' feature has an attribute containing the proj4 string for the printed frame, for example EPSG:32614 is:

+proj=utm +zone=14 +datum=WGS84 +units=m +no_defs

I want to calculate the rotation angle of each polygon frame, and I found this excellent post that does just that. But, of course this:

(180/pi())* atan2($y_at(1)-$y_at(0), $x_at(1)-$x_at(0))

calculates the rotation angle of the feature based on the EPSG:3857 SRS. For example, one of the output angles was -55.5 degrees, but when I extract this feature and run the calculation with the dataset defined as 32614, the angle is calculated at -53.9 degrees.

Is it possible to use this proj4 string attribute to calculate the rotation angle of the feature with the EPSG:32614 instead of the EPSG:3857 of the table itself?

3
  • if you had the EPSG code in a field, it would be easy with the transform function.
    – csk
    Nov 13 '19 at 17:43
  • I looked at the transform function, but can't quite sort it out and it's coming up invalid in the field calc dialogue: (transform( 180/pi())* atan2($y_at(1)-$y_at(0), $x_at(1)-$x_at(0)), 'EPSG:3857', 'EPSG:32614')
    – jamierob
    Nov 13 '19 at 17:52
  • Normally you would do transform($geometry, 'EPSG:3857', 'EPSG:32614'), and substitute that into your equation wherever you used the $geometry term. But your expression uses $x_at and $y_at, which always apply to the current geometry, so there's nowhere to substitute the transformed geometry. This will require some creative problem solving.
    – csk
    Nov 13 '19 at 18:02
4

Normally, in order to do a calculation using a geometry in a different CRS than the layer it's in, you would transform the geometry like this: transform($geometry, 'EPSG:3857', 'EPSG:32614'). Then you would substitute the transformed geometry into your equation wherever you used the $geometry term. But your expression uses $x_at and $y_at, which always apply to the current geometry, so there's nowhere to substitute the transformed geometry.

Instead, let's store the two vertices in the attribute table. In the process, we'll transform each of these vertices into the target CRS.

  1. This expression creates a point geometry from the first vertex of the geometry: make_point($x_at(0),$y_at(0)

  2. Use the transform function to re-project the point into the target CRS. I'm assuming your layer has a field called "EPSG" which contains the EPSG code for the target CRS. transform(make_point($x_at(0),$y_at(0)),'EPSG:3857',"EPSG")

  3. Finally, convert the point geometry to a string so it can be stored in the attribute table, using the geom_to_wkt function.

    • Make a field called "Vertex_0" with this expression: geom_to_wkt(transform(make_point($x_at(0),$y_at(0)),'EPSG:3857',"EPSG")))
    • Make a field called "Vertex_1" with this expression: geom_to_wkt(transform(make_point($x_at(1),$y_at(1)),'EPSG:3857',"EPSG")))

To convert the WKT string back into a geometry, use geom_from_wkt. Use the functions x(geometry) and y(geometry` to extract the x and y coordinates.

Make the following substitutions into the original equation:

  • substitute x(geom_from_wkt("Vertex_0")) where the original equation says $x_at(0)
  • substitute y(geom_from_wkt("Vertex_0")) where the original equation says $y_at(0)
  • substitute x(geom_from_wkt("Vertex_1")) where the original equation says $x_at(1)
  • substitute y(geom_from_wkt("Vertex_1")) where the original equation says $y_at(1)

New equation:

(180/pi())* atan2(y(geom_from_wkt("Vertex_1"))-y(geom_from_wkt("Vertex_0")), x(geom_from_wkt("Vertex_1"))-x(geom_from_wkt("Vertex_0")))

Note: Of course you could do this all in one step if you substitute x(transform(make_point($x_at(0),$y_at(0)),'EPSG:3857',"EPSG")) for $x_at(0) (and so on) in the original equation. Storing the vertices as fields just makes it easier to troubleshoot.

1
  • 1
    That worked beautifully. I continue to be astounded by the knowledge and generosity of people on the gis stackexchange.
    – jamierob
    Nov 13 '19 at 20:39

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.