[PROJ] Migrating to proj API 6+ without forcing crs input
Johannes Schauer Marin Rodrigues
josch at mister-muffin.de
Mon Apr 18 06:43:16 PDT 2022
Hi Even,
Quoting Even Rouault (2022-04-14 19:31:07)
> Le 14/04/2022 à 18:55, Johannes Schauer Marin Rodrigues a écrit :
> > thanks for your reply! So you suggest that projects that are upgrading to proj
> > API 6+ and want to continue supporting non-CRS user-supplied strings should
> > resort to implementing their own proj string parser and transformation code?
>
> Well, tinshift support was added after PROJ 6, so there's no real
> backward compatibility issue here.
ah I didn't mean it as "I want to run this on old PROJ". I'm happy to run the
most recent version. I just wanted to keep the stuff that worked with software
(mapnik) that used to use the old proj API also working with the new API.
> Most of the past users of pj_transform() used it with PROJ strings that were
> CRS definitions, so we're a bit in the category of unanticipated / marginal
> use cases.
Yes, I realize that. I'm thus even more grateful that you took the time and
replied to my marginal use-case query! :)
> The reason for PROJ 6 to differenciate CRS from
> transformations/pipelines is that in past version PROJ strings that were
> meant to be CRS definitions used PROJ operation methods that can also be
> used to do transformations, hence it was a bit like a CRS was a
> transformation from something (some 'base' CRS not always well defined)
> to the CRS of interest. Stopping here my confusing digging into past
> design decisions :-), but basically there are good reasons of separating
> CRS and transformation/pipelines and not lightly mixing them together.
Got it. I'm not trying to challenge that design decision. Especially since I've
got it working now with your help. :)
> > Is there not something that can be provided by proj itself to make this
> > easier? If proj would add some code to support this, then this would also
> > prevent a number of differently buggy proj parsing and transformation
> > codebases to be implemented by every project that cares about non-CRS
> > input. It seems odd to me that with upgrading to a new API version, the
> > use-case of using proj for non-CRS input is now so hard where before it
> > just worked out-of-the box.
>
> Obviously it is messy to do that outside of PROJ . It's not immediate to
> my mind which reasonable and well-defined C API PROJ could offer for
> what you need.
Lucky for me, the project I'm interested in (mapnik) is C++. :D
> But there's definitely in the PROJ C++ code things to ingest a PROJ
> string, invert it (including reverting order and direction of pipeline
> steps), concatenate it, in a hopefully not-too-buggy way, as that's used
> internally by it.
>
> Particularly the PROJStringFormatter class:
> https://github.com/OSGeo/PROJ/blob/master/include/proj/io.hpp#L376
>
> with the ingestPROJString(), startInversion(), stopInversion(),
> addStep(), addParam() and toString() methods
>
> In pseudo code, your below proposal would be something like
>
> osgeo::proj::io::PROJStringFormatter formatter;
>
> formatter.startInversion();
>
> formatter.ingestPROJString(source_pipeline);
>
> formatter.stopInversion();
>
> formatter.ingestPROJString(target_pipeline)
>
> P = proj_create(ctx, formatter.toString().c_str() )
Wow, thanks a ton! This indeed gave me all the tools to get this to work. :D
> > What to do if source or target is not a CRS and also not a pipeline?
>
> What kind of objects are you thinking too ? If it is a single
> transformation step, then it is semantically a pipeline with one step.
>
> You can't do anything useful with other types of objects (datums ,
> ellipsoids, etc.)
In the mapnik test cases I used "+proj=axisswap +order=2,1" as a source that is
not a CRS and it seems to work fine.
> Mixing a CRS with a pipeline/transformation can also be a bit tricky
> because of unit and axis order adjustments that you may need to insert.
I think those users who supply custom non-CRS projection should accept that
they are on their own and should not expect any magic to happen when they don't
provide something with the right unit or axis adjustments.
In any case, for anybody who might be interested, here is the code that I came
up with and which seems to work for mapnik:
https://github.com/mapnik/mapnik/pull/4309
If you have a stackoverflow account and care about the bounty, feel free to
copypaste the following code into an answer so that I can award the bounty:
https://stackoverflow.com/questions/71823105/migrating-to-proj-api-6-without-forcing-crs-input
Anyways, here is the code. I'm happy about any comments you can give and
otherwise leave the code here for posterity in case somebody else is facing the
same problem as I. The return value of source.params() and dest.params() is the
proj string. The function throw_proj_exception calls proj_context_errno_string
and is omitted for brevity. The resulting projection is stored in transform_ in
both codepaths.
PJ* transform_ = nullptr;
ctx_ = proj_context_create();
proj_log_level(ctx_, PJ_LOG_ERROR);
// we replicate what proj_create_crs_to_crs() does and then call
// proj_create_crs_to_crs_from_pj because we need a PJ object for
// proj_is_crs and we want to avoid running proj_create twice as
// as much (once ourselves and once in proj_create_crs_to_crs)
std::string crs_source = pj_add_type_crs_if_needed(source.params());
std::string crs_dest = pj_add_type_crs_if_needed(dest.params());
PJ* src;
PJ* dst;
try {
src = proj_create(ctx_, crs_source.c_str());
} catch( const std::exception& ) {
src = nullptr;
}
if (!src) {
throw_proj_exception(ctx_,
std::string("Cannot instantiate source projection for '") + source.params() + "'");
}
try {
dst = proj_create(ctx_, crs_dest.c_str());
} catch( const std::exception& ) {
dst = nullptr;
}
if (!dst) {
throw_proj_exception(ctx_,
std::string("Cannot instantiate dest projection for '") + dest.params() + "'");
}
if (proj_is_crs(src) && proj_is_crs(dst)) {
// if both source and destination are a CRS, we can use
// proj_create_crs_to_crs
transform_ = proj_create_crs_to_crs_from_pj(ctx_, src, dst, nullptr, nullptr);
if (transform_ == nullptr)
{
throw_proj_exception(ctx_,
std::string("Cannot initialize proj_transform (crs_to_crs) for given projections: '") +
source.params() + "'->'" + dest.params() + "'");
}
PJ* transform_gis = proj_normalize_for_visualization(ctx_, transform_);
if (transform_gis == nullptr)
{
throw_proj_exception(ctx_,
std::string("Cannot initialize proj_transform (normalize) for given projections: '") +
source.params() + "'->'" + dest.params() + "'");
}
proj_destroy(transform_);
transform_ = transform_gis;
} else {
// if one of the projections is not a CRS, then this was
// probably a user-supplied transform and we have to create
// the transform pipeline ourselves
auto formatter = osgeo::proj::io::PROJStringFormatter::create();
formatter->startInversion();
try {
formatter->ingestPROJString(source.params());
} catch (const osgeo::proj::io::ParsingException &e) {
throw_proj_exception(ctx_, std::string("Failed to ingest source: '") + source.params() + "'");
}
formatter->stopInversion();
try {
formatter->ingestPROJString(dest.params());
} catch (const osgeo::proj::io::ParsingException &e) {
throw_proj_exception(ctx_, std::string("Failed to ingest dest: '") + dest.params() + "'");
}
try {
transform_ = proj_create(ctx_, formatter->toString().c_str());
} catch( const std::exception& ) {
transform_ = nullptr;
}
if (!transform_) {
throw_proj_exception(ctx_,
std::string("Cannot instantiate custom projection pipeline: '") + formatter->toString() + "'");
}
}
proj_destroy(src);
proj_destroy(dst);
return transform_;
Since I need to call proj_is_crs to find out if source and dest are a CRS or
not, and since I want to avoid calling proj_create more often than necessary,
the code recreates parts of whatproj_create_crs_to_crs does and uses
proj_create_crs_to_crs_from_pj instead if both source and dest are a CRS. For
this to work, I added my own version of pj_add_type_crs_if_needed because that
function wasn't exposed in proj.h.
Thanks a ton for your super helpful insights!
cheers, josch
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 833 bytes
Desc: signature
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20220418/6d0ac2f8/attachment.sig>
More information about the PROJ
mailing list