Migrating to proj API 6+ without forcing crs input

Question

I'm posting this question here instead of gis.stackexchange.com because I think this is more of a programming question than a GIS question.

The proj migration guide explains how to convert code to proj 6 and later by turning code that looked like this (error handling omitted for brevity):

projPJ pj_longlat = pj_init_plus(source);
projPJ pj_merc = pj_init_plus(target);
pj_transform(pj_longlat, pj_merc, ...)

into code that looks like this:

PJ *P = proj_create_crs_to_crs(PJ_DEFAULT_CTX, source, target, NULL);
proj_trans(P, PJ_FWD, ...);

This is what was done in several projects including mapnik but there is one downside: with the new code, both source and target must be a CRS or otherwise proj_create_crs_to_crs will fail with:

proj_create_operations: source_crs is not a CRS

This means that some applications (for example I successfully used mapnik rendering with a +proj=pipeline projection, which is not a CRS) are not possible anymore.

In a proj issue report a developer points out that it should be possible to avoid proj_create_crs_to_crs by using proj_create manually. But if I understand it correctly, then to create a PJ object that transforms from source to target (and vice-versa) I need a pipeline but if one of my projections contains a pipeline, then this will not work because I cannot have nested pipelines.

So how would I convert code from the PROJ 4 to PROJ 6 API but without removing the ability to input arbitrary transformations, even those that are not a CRS?

Edit

I asked the same question on the PROJ mailing list and got this very helpful reply from Even Rouault. I'll use that to try and implement a solution.

Answer 1

With help from Even Rouault I managed to submit a pull request to mapnik which will use an alternative code path if the input transforms are not a CRS:

https://github.com/mapnik/mapnik/pull/4309

Essentially, this is doing the following (error checking omitted for brevity):

ctx_ = proj_context_create();
std::string crs_source = pj_add_type_crs_if_needed(source);
PJ *src = proj_create(ctx_, crs_source.c_str());
std::string crs_dest = pj_add_type_crs_if_needed(dest);
PJ *dst = proj_create(ctx_, crs_dest.c_str());
if (proj_is_crs(src) && proj_is_crs(dst)) {
    /* insert proj_create_crs_to_crs and friends from
     * the official migration guide here */
} else {
    /* either src or dst is not a crs, so we need to set
     * up the pipeline ourselves */
    auto formatter = osgeo::proj::io::PROJStringFormatter::create();
    formatter->startInversion();
    formatter->ingestPROJString(source);
    formatter->stopInversion();
    formatter->ingestPROJString(dest);
    transform_ = proj_create(ctx_, formatter->toString().c_str());
}
proj_destroy(src);
proj_destroy(dst);
return transform_;

Most notably, the above alternative codepath does not call proj_normalize_for_visualization, so the user is on their own in case they provide a non-CRS transform that would need normalization.

The function pj_add_type_crs_if_needed comes from the PROJ codebase in src/4D_api.cpp and adds a +type=crs suffix to a PROJ string (if it is a PROJ string) and can be written like this:

std::string pj_add_type_crs_if_needed(const std::string &str) {
    std::string ret(str);
    if ((str.rfind("proj=", 0) == 0 || str.rfind("+proj=", 0) == 0 ||
         str.rfind("+init=", 0) == 0 || str.rfind("+title=", 0) == 0) &&
        str.find("type=crs") == std::string::npos) {
        ret += " +type=crs";
    }
    return ret;
}