[Proj] vertical projection

[Even Rouault]

Le mercredi 05 octobre 2016 10:21:31, Vincent Mora a écrit :
> Hi,
> 
> I need to draw terrain sections for several projects (hydrology, geology).
> 
> The vertical section is defined by a linestring in the (x,y) plane (e.g.
> +proj=lcc). Data in the vicinity of the section are ortho-projected on
> the folded plane with coordinates (s,z), where s is the point location
> on the line (meters or feet).
> 
> For the moment the CRS part of the QGis plugin we are developping is a
> bit of a hack, and I'd like to introduce that kind of projection in
> proj4 since I believe it belongs there.
> 
> A projection on a folded plane is rather common, but in this case:
> - the folds are sharp,
> - the projection is not well defined for points that lie on the bisector
> of two line-segments
> - the projection is meaninigless for points that are not in the plane
> vicinity.
> 
> The proj4 string should contain:
> - the line definition (e.g. +wkt=LINESTRING(...))
> - the classical definition of the CRS of the line,
> - and a thickness (vicinity of the plane)
> 
> Questions:
> - has such a CRS it's place in proj4 ?
> - do you see obvious errors/misconceptions ?

Vincent,

The fact that the line definition can be defined in any CRS handled by proj.4
would make it quite a odd object for proj.4 if we used current mechanisms.
You'd need to have a definition like :
+proj=vertplane +wkt=LINESTRING(...) +wkt_proj=lcc +wkt_proj_lat_1= ...

Probably a better fit would be to rely on Thomas Knudsen's ongoing work on
transformation pipelines (
http://osgeo-org.1560.x6.nabble.com/Transformation-pipelines-your-opinion-td5269960.html ,
http://thomasknudsen.net/tran_users_guide.html ).
I'm not sure this has yet landed into master. Perhaps Thomas can confirm ?

With that framework, I think you could define a 2 stage pipeline :
- first converting from geographical coordinates to LCC projected ones
- then doing the projection onto the vertical plane

+proj=pipeline +step=lcc +lat_1=...  +step +proj=vertplane +wkt=LINESTRING(...) +thickness=

Regarding QGIS integration, this would require some changes, due to pipeline support
requiring a new dedicated API in proj.4 "TRIPLET pj_apply_projection (TRIPLET point, int direction, PJ *P)" since you cannot use the
existing int pj_transform( projPJ srcdefn, projPJ dstdefn, long point_count, int point_offset, double *x, double *y, double *z )

By the way, Thomas, wouldn't that make sense to have a pipeline transformation
API that would accept several points at once ? The default implementation could
just loop over them, but that would be a provision for potential later optimizations
where you could potentially vectorize computations.

Even

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com