[PROJ] Proj transformation for arbitrary (quaternion/Euler) spatial rotation

[DAVEN P QUINN]

Hello all,

I’m a geologist who is working with paleogeographic reconstructions (see example here<https://davenquinn.com/viz/corelle-demo-pbdb/?time=295>). These require composite reproductions where different parts of a feature dataset (continents) are rotated with different axis-angle transformations. I’ve had good luck doing these rotations with quaternion math in the browser/Python environments, but I am now trying to use Proj transformations in order to apply rotations directly within PostGIS queries.

This linked image<https://pbs.twimg.com/media/Fgpx4UGXEAQmPOt?format=jpg&name=4096x4096> shows the desired result, a plate reconstruction to 250 Ma with a different rotation applied to each plate. This was produced by applying the desired quaternions through pl/pgsql math. Moving this math to Proj internals would result in a ~50-100x speedup.

The tool that seems most fitting is the `ob_trans` family of projections, but I have been unable to define a rotation that can handle my preferred representation (a pole defined in Lon-lat coordinates and an associated angle of rotation around it). From my reading of the docs, I believe the `o_lon_c`, `o_lat_c` and `o_alpha` parameters should do this, but I cannot get them to work reliably. In fact, I can’t even reliably define a ’no-op’ transformation that leaves coordinates unchanged. Perhaps I have the angular coordinate system wrong, or the rotation is being done in Cartesian space even though `proj=lonlat` is used.

Is it possible to define an arbitrary spherical spatial rotation in Proj transformations? Maybe I need to use pipelines instead? I’d appreciate any guidance! More details below the fold...

Regards,

Daven P. Quinn
Research scientist II · U of Wisconsin Madison
PhD · structural geology · Caltech ‘18
https://davenquinn.com
+1 704 920 8487

-------------

Here is an example of the pl/pgsql I am currently using to assemble a projection (I have gone through many iterations testing different offsets and the different ways to specify the transformation):
```
RETURN '+proj=ob_tran +o_proj=longlat +o_alpha=' || pi()/2+angle || 'r +o_lon_c=' || pi()/2+lon || 'r +o_lat_c=' || lat || 'r' proj
```
where (lat, lon, angle) defines a rotation pole.

This results in some geometries (with a fortuitous set of poles, I guess) attaining ballpark-correct transformations while other features are shifted far outside of their origin tiles.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20221109/982c84e3/attachment.htm>