I need to plot 3d PostGIS data into vertical planes for symbolization and drawing vertical cross-section interpretations for mineral exploration models. I have written a short plpgsql function that takes the PostGIS 3D geometry, moves an xy-coordinate to (0,0,0) then rotates it on the z-axis based upon a given rotation(degrees) and then rotates it vertically along the x-axis. This new geometry is then filtered/cropped by the new/rotated z-coordinates so just the right plane of data can be viewed. It works and functions well with QGIS but the downside is that it is not a very update able view for any points/lines/polygons that I draw for any particular cross-section interpretation.
Is it possible to write a proj4/proj6 project that rotates the data on the x-axis? It would have to be version of the Oblique Mercator to provide the horizontal/z-axis rotation?
I know I am asking for a lot but vertical cross-section interpretation is crucial in the mineral/petroleum industry.
CREATE OR REPLACE FUNCTION public.xsec_vert_rotate(geometry, double precision, double precision, double precision) RETURNS geometry AS $BODY$ DECLARE orig_geom alias for $1; x1 alias for $2; y1 alias for $3; rot_deg alias for $4; rot_geom geometry; BEGIN rot_geom := ST_RotateX(ST_RotateZ(ST_Translate(orig_geom,(0.0::double precision - x1),(0.0::double precision - y1), >0),radians(rot_deg)),-pi()/2); RETURN rot_geom; END; $BODY$ LANGUAGE plpgsql VOLATILE COST 100; COMMENT ON FUNCTION public.xsec_vert_rotate(geometry, double precision, double precision, double precision) IS 'Cross-Section (Vertical) Rotate: geometry3d, x1/y1 Point Coordinates to move to (0,0,0), Rotation_degress along z-axis to orient data to xy grid system, sets SRID to Zero'; SELECT l.sampleid, l.holeid, l.sampfrom, l.sampto, st_z(xsec_vert_rotate(st_transform(st_startpoint(l.bhline_geom), 3158), 375800, 6307000, -65)) AS z_xsec, xsec_vert_rotate(st_transform(l.bhline_geom, 3158), 375800, 6307000, 65) AS bhline_geom FROM public.lab WHERE st_z(xsec_vert_rotate(st_transform(st_startpoint(l.bhline_geom), 3158), 375800, 6307000, -65)) between -50 and 50;