do_proj.c

Go to the documentation of this file.

/**
   \file do_proj.c
 
   \brief GProj library - Functions for re-projecting point data
 
   \author Original Author unknown, probably Soil Conservation Service
   Eric Miller, Paul Kelly, Markus Metz
 
   (C) 2003-2008,2018 by the GRASS Development Team
 
   This program is free software under the GNU General Public
   License (>=v2). Read the file COPYING that comes with GRASS
   for details.
**/
 
#include <stdio.h>
#include <string.h>
#include <ctype.h>
#include <math.h>
 
#include <grass/gis.h>
#include <grass/gprojects.h>
#include <grass/glocale.h>
 
/* a couple defines to simplify reading the function */
#define MULTIPLY_LOOP(x, y, c, m) \
    do {                          \
        for (i = 0; i < c; ++i) { \
            x[i] *= m;            \
            y[i] *= m;            \
        }                         \
    } while (0)
 
#define DIVIDE_LOOP(x, y, c, m)   \
    do {                          \
        for (i = 0; i < c; ++i) { \
            x[i] /= m;            \
            y[i] /= m;            \
        }                         \
    } while (0)
 
static double METERS_in = 1.0, METERS_out = 1.0;
 
int get_pj_area(const struct pj_info *iproj, double *xmin, double *xmax,
                double *ymin, double *ymax)
{
    struct Cell_head window;
 
    /* modules must set the current window, do not unset this window here */
    /* G_unset_window(); */
    G_get_set_window(&window);
    *xmin = window.west;
    *xmax = window.east;
    *ymin = window.south;
    *ymax = window.north;
 
    if (window.proj != PROJECTION_LL) {
        /* transform to ll equivalent */
        double estep, nstep;
        double x[85], y[85];
        int i;
        const char *projstr = NULL;
        char *indef = NULL;
        struct pj_info oproj, tproj; /* proj parameters  */
 
        oproj.pj = NULL;
        oproj.proj[0] = '\0';
        tproj.def = NULL;
 
        if (proj_get_type(iproj->pj) == PJ_TYPE_BOUND_CRS) {
            PJ *source_crs;
 
            source_crs = proj_get_source_crs(NULL, iproj->pj);
            if (source_crs) {
                projstr =
                    proj_as_proj_string(NULL, source_crs, PJ_PROJ_5, NULL);
                if (projstr) {
                    indef = G_store(projstr);
                }
                proj_destroy(source_crs);
            }
        }
        else {
            projstr = proj_as_proj_string(NULL, iproj->pj, PJ_PROJ_5, NULL);
            if (projstr) {
                indef = G_store(projstr);
            }
        }
 
        if (indef == NULL)
            indef = G_store(iproj->def);
 
        /* needs +over to properly cross the anti-meridian
         * the +over switch can be used to disable the default wrapping
         * of output longitudes to the range -180 to 180 */
        G_asprintf(&tproj.def, "+proj=pipeline +step +inv %s +over", indef);
        G_debug(1, "get_pj_area() tproj.def: %s", tproj.def);
        tproj.pj = proj_create(PJ_DEFAULT_CTX, tproj.def);
 
        if (tproj.pj == NULL) {
            G_warning(_("proj_create() failed for '%s'"), tproj.def);
            G_free(indef);
            G_free(tproj.def);
            proj_destroy(tproj.pj);
 
            return 0;
        }
        projstr = proj_as_proj_string(NULL, tproj.pj, PJ_PROJ_5, NULL);
        if (projstr == NULL) {
            G_warning(_("proj_create() failed for '%s'"), tproj.def);
            G_free(indef);
            G_free(tproj.def);
            proj_destroy(tproj.pj);
 
            return 0;
        }
        else {
            G_debug(1, "proj_create() projstr '%s'", projstr);
        }
        G_free(indef);
 
        /* inpired by gdal/ogr/ogrct.cpp OGRProjCT::ListCoordinateOperations()
         */
        estep = (window.east - window.west) / 21.;
        nstep = (window.north - window.south) / 21.;
        for (i = 0; i < 20; i++) {
            x[i] = window.west + estep * (i + 1);
            y[i] = window.north;
 
            x[i + 20] = window.west + estep * (i + 1);
            y[i + 20] = window.south;
 
            x[i + 40] = window.west;
            y[i + 40] = window.south + nstep * (i + 1);
 
            x[i + 60] = window.east;
            y[i + 60] = window.south + nstep * (i + 1);
        }
        x[80] = window.west;
        y[80] = window.north;
        x[81] = window.west;
        y[81] = window.south;
        x[82] = window.east;
        y[82] = window.north;
        x[83] = window.east;
        y[83] = window.south;
        x[84] = (window.west + window.east) / 2.;
        y[84] = (window.north + window.south) / 2.;
 
        GPJ_transform_array(iproj, &oproj, &tproj, PJ_FWD, x, y, NULL, 85);
 
        proj_destroy(tproj.pj);
        G_free(tproj.def);
 
        *xmin = *xmax = x[84];
        *ymin = *ymax = y[84];
        for (i = 0; i < 84; i++) {
            if (*xmin > x[i])
                *xmin = x[i];
            if (*xmax < x[i])
                *xmax = x[i];
            if (*ymin > y[i])
                *ymin = y[i];
            if (*ymax < y[i])
                *ymax = y[i];
        }
 
        /* The west longitude is generally lower than the east longitude,
         * except for areas of interest that go across the anti-meridian.
         * do not reduce global coverage to a small north-south strip
         */
        if (*xmin < -180 && *xmax < 180 && *xmin + 360 > *xmax) {
            /* must be crossing the anti-meridian at 180W */
            *xmin += 360;
        }
        else if (*xmax > 180 && *xmin > -180 && *xmax - 360 < *xmin) {
            /* must be crossing the anti-meridian at 180E */
            *xmax -= 360;
        }
 
        G_debug(1, "input window north: %.8f", window.north);
        G_debug(1, "input window south: %.8f", window.south);
        G_debug(1, "input window east: %.8f", window.east);
        G_debug(1, "input window west: %.8f", window.west);
 
        G_debug(1, "transformed xmin: %.8f", *xmin);
        G_debug(1, "transformed xmax: %.8f", *xmax);
        G_debug(1, "transformed ymin: %.8f", *ymin);
        G_debug(1, "transformed ymax: %.8f", *ymax);
 
        /* test valid values, as in
         * gdal/ogr/ogrct.cpp
         * OGRCoordinateTransformationOptions::SetAreaOfInterest()
         */
        if (fabs(*xmin) > 180) {
            G_warning(_("Invalid west longitude %g"), *xmin);
            return 0;
        }
        if (fabs(*xmax) > 180) {
            G_warning(_("Invalid east longitude %g"), *xmax);
            return 0;
        }
        if (fabs(*ymin) > 90) {
            G_warning(_("Invalid south latitude %g"), *ymin);
            return 0;
        }
        if (fabs(*ymax) > 90) {
            G_warning(_("Invalid north latitude %g"), *ymax);
            return 0;
        }
        if (*ymin > *ymax) {
            G_warning(_("South %g is larger than north %g"), *ymin, *ymax);
            return 0;
        }
    }
    G_debug(1, "get_pj_area(): xmin %g, xmax %g, ymin %g, ymax %g", *xmin,
            *xmax, *ymin, *ymax);
 
    return 1;
}
 
char *get_pj_type_string(PJ *pj)
{
    char *pj_type = NULL;
 
    switch (proj_get_type(pj)) {
    case PJ_TYPE_UNKNOWN:
        G_asprintf(&pj_type, "unknown");
        break;
    case PJ_TYPE_ELLIPSOID:
        G_asprintf(&pj_type, "ellipsoid");
        break;
    case PJ_TYPE_PRIME_MERIDIAN:
        G_asprintf(&pj_type, "prime meridian");
        break;
    case PJ_TYPE_GEODETIC_REFERENCE_FRAME:
        G_asprintf(&pj_type, "geodetic reference frame");
        break;
    case PJ_TYPE_DYNAMIC_GEODETIC_REFERENCE_FRAME:
        G_asprintf(&pj_type, "dynamic geodetic reference frame");
        break;
    case PJ_TYPE_VERTICAL_REFERENCE_FRAME:
        G_asprintf(&pj_type, "vertical reference frame");
        break;
    case PJ_TYPE_DYNAMIC_VERTICAL_REFERENCE_FRAME:
        G_asprintf(&pj_type, "dynamic vertical reference frame");
        break;
    case PJ_TYPE_DATUM_ENSEMBLE:
        G_asprintf(&pj_type, "datum ensemble");
        break;
 
    /** Abstract type, not returned by proj_get_type() */
    case PJ_TYPE_CRS:
        G_asprintf(&pj_type, "crs");
        break;
    case PJ_TYPE_GEODETIC_CRS:
        G_asprintf(&pj_type, "geodetic crs");
        break;
    case PJ_TYPE_GEOCENTRIC_CRS:
        G_asprintf(&pj_type, "geocentric crs");
        break;
 
    /** proj_get_type() will never return that type, but
     * PJ_TYPE_GEOGRAPHIC_2D_CRS or PJ_TYPE_GEOGRAPHIC_3D_CRS. */
    case PJ_TYPE_GEOGRAPHIC_CRS:
        G_asprintf(&pj_type, "geographic crs");
        break;
    case PJ_TYPE_GEOGRAPHIC_2D_CRS:
        G_asprintf(&pj_type, "geographic 2D crs");
        break;
    case PJ_TYPE_GEOGRAPHIC_3D_CRS:
        G_asprintf(&pj_type, "geographic 3D crs");
        break;
    case PJ_TYPE_VERTICAL_CRS:
        G_asprintf(&pj_type, "vertical crs");
        break;
    case PJ_TYPE_PROJECTED_CRS:
        G_asprintf(&pj_type, "projected crs");
        break;
    case PJ_TYPE_COMPOUND_CRS:
        G_asprintf(&pj_type, "compound crs");
        break;
    case PJ_TYPE_TEMPORAL_CRS:
        G_asprintf(&pj_type, "temporal crs");
        break;
    case PJ_TYPE_ENGINEERING_CRS:
        G_asprintf(&pj_type, "engineering crs");
        break;
    case PJ_TYPE_BOUND_CRS:
        G_asprintf(&pj_type, "bound crs");
        break;
    case PJ_TYPE_OTHER_CRS:
        G_asprintf(&pj_type, "other crs");
        break;
    case PJ_TYPE_CONVERSION:
        G_asprintf(&pj_type, "conversion");
        break;
    case PJ_TYPE_TRANSFORMATION:
        G_asprintf(&pj_type, "transformation");
        break;
    case PJ_TYPE_CONCATENATED_OPERATION:
        G_asprintf(&pj_type, "concatenated operation");
        break;
    case PJ_TYPE_OTHER_COORDINATE_OPERATION:
        G_asprintf(&pj_type, "other coordinate operation");
        break;
    default:
        G_asprintf(&pj_type, "unknown");
        break;
    }
 
    return pj_type;
}
 
PJ *get_pj_object(const struct pj_info *in_gpj, char **in_defstr)
{
    PJ *in_pj = NULL;
 
    *in_defstr = NULL;
 
    /* 1. SRID, 2. WKT, 3. standard pj from pj_get_kv */
    if (in_gpj->srid) {
        G_debug(1, "Trying SRID '%s' ...", in_gpj->srid);
        in_pj = proj_create(PJ_DEFAULT_CTX, in_gpj->srid);
        if (!in_pj) {
            G_warning(_("Unrecognized SRID '%s'"), in_gpj->srid);
        }
        else {
            *in_defstr = G_store(in_gpj->srid);
            /* PROJ will do the unit conversion if set up from srid
             * -> disable unit conversion for GPJ_transform */
            /* ugly hack */
            ((struct pj_info *)in_gpj)->meters = 1;
        }
    }
    if (!in_pj && in_gpj->wkt) {
        G_debug(1, "Trying WKT '%s' ...", in_gpj->wkt);
        in_pj = proj_create(PJ_DEFAULT_CTX, in_gpj->wkt);
        if (!in_pj) {
            G_warning(_("Unrecognized WKT '%s'"), in_gpj->wkt);
        }
        else {
            *in_defstr = G_store(in_gpj->wkt);
            /* PROJ will do the unit conversion if set up from wkt
             * -> disable unit conversion for GPJ_transform */
            /* ugly hack */
            ((struct pj_info *)in_gpj)->meters = 1;
        }
    }
    if (!in_pj && in_gpj->pj) {
        in_pj = proj_clone(PJ_DEFAULT_CTX, in_gpj->pj);
        *in_defstr = G_store(proj_as_wkt(NULL, in_pj, PJ_WKT2_LATEST, NULL));
        if (*in_defstr && !**in_defstr)
            *in_defstr = NULL;
    }
 
    if (!in_pj) {
        G_warning(_("Unable to create PROJ object"));
 
        return NULL;
    }
 
    /* Even Rouault:
     * if info_in->def contains a +towgs84/+nadgrids clause,
     * this pipeline would apply it, whereas you probably only want
     * the reverse projection, and no datum shift.
     * The easiest would probably to mess up with the PROJ string.
     * Otherwise with the PROJ API, you could
     * instantiate a PJ object from the string,
     * check if it is a BoundCRS with proj_get_source_crs(),
     * and in that case, take the source CRS with proj_get_source_crs(),
     * and do the inverse transform on it */
 
    if (proj_get_type(in_pj) == PJ_TYPE_BOUND_CRS) {
        PJ *source_crs;
 
        G_debug(1, "found bound crs");
        source_crs = proj_get_source_crs(NULL, in_pj);
        if (source_crs) {
            *in_defstr =
                G_store(proj_as_wkt(NULL, source_crs, PJ_WKT2_LATEST, NULL));
            if (*in_defstr && !**in_defstr)
                *in_defstr = NULL;
            in_pj = source_crs;
        }
    }
 
    return in_pj;
}
 
/**
 * \brief Create a PROJ transformation object to transform coordinates
 *        from an input SRS to an output SRS
 *
 * After the transformation has been initialized with this function,
 * coordinates can be transformed from input SRS to output SRS with
 * GPJ_transform() and direction = PJ_FWD, and back from output SRS to
 * input SRS with direction = OJ_INV.
 * If coordinates should be transformed between the input SRS and its
 * latlong equivalent, an uninitialized info_out with
 * info_out->pj = NULL can be passed to the function. In this case,
 * coordinates will be transformed between the input SRS and its
 * latlong equivalent, and for PROJ 5+, the transformation object is
 * created accordingly
 *
 * PROJ 5+:
 *   info_in->pj must not be null
 *   if info_out->pj is null, assume info_out to be the ll equivalent
 *   of info_in
 *   create info_trans as conversion from info_in to its ll equivalent
 *   NOTE: this is the inverse of the logic of PROJ 5 which by default
 *         converts from ll to a given SRS, not from a given SRS to ll
 *         thus PROJ 5+ itself uses an inverse transformation in the
 *         first step of the pipeline for proj_create_crs_to_crs()
 *   if info_trans->def is not NULL, this pipeline definition will be
 *   used to create a transformation object
 *
 * \param info_in pointer to pj_info struct for input co-ordinate system
 * \param info_out pointer to pj_info struct for output co-ordinate system
 * \param info_trans pointer to pj_info struct for a transformation object (PROJ
 *5+)
 *
 * \return 1 on success, -1 on failure
 **/
int GPJ_init_transform(const struct pj_info *info_in,
                       const struct pj_info *info_out,
                       struct pj_info *info_trans)
{
    if (info_in->pj == NULL)
        G_fatal_error(_("Input coordinate system is NULL"));
 
    if (info_in->def == NULL)
        G_fatal_error(_("Input coordinate system definition is NULL"));
 
    /* PROJ6+: enforce axis order easting, northing
     * +axis=enu (works with proj-4.8+) */
 
    info_trans->pj = NULL;
    info_trans->meters = 1.;
    info_trans->zone = 0;
    snprintf(info_trans->proj, sizeof(info_trans->proj), "pipeline");
 
    /* user-provided pipeline */
    if (info_trans->def) {
        const char *projstr;
 
        /* info_in->pj, info_in->proj, info_out->pj, info_out->proj
         * must be set */
        if (!info_in->pj || !info_in->proj[0] || !info_out->pj ||
            !info_out->proj[0]) {
            G_warning(_(
                "A custom pipeline requires input and output projection info"));
 
            return -1;
        }
 
        /* create a pj from user-defined transformation pipeline */
        info_trans->pj = proj_create(PJ_DEFAULT_CTX, info_trans->def);
        if (info_trans->pj == NULL) {
            G_warning(_("proj_create() failed for '%s'"), info_trans->def);
 
            return -1;
        }
        projstr = proj_as_proj_string(NULL, info_trans->pj, PJ_PROJ_5, NULL);
        if (projstr == NULL) {
            G_warning(_("proj_create() failed for '%s'"), info_trans->def);
 
            return -1;
        }
        else {
            /* make sure axis order is easting, northing
             * proj_normalize_for_visualization() does not work here
             * because source and target CRS are unknown to PROJ
             * remove any "+step +proj=axisswap +order=2,1" ?
             *  */
            info_trans->def = G_store(projstr);
 
            if (strstr(info_trans->def, "axisswap")) {
                G_warning(
                    _("The transformation pipeline contains an '%s' step. "
                      "Remove this step if easting and northing are swapped in "
                      "the output."),
                    "axisswap");
            }
 
            G_debug(1, "proj_create() pipeline: %s", info_trans->def);
 
            /* the user-provided PROJ pipeline is supposed to do
             * all the needed unit conversions */
            /* ugly hack */
            ((struct pj_info *)info_in)->meters = 1;
            ((struct pj_info *)info_out)->meters = 1;
        }
    }
    /* if no output CRS is defined,
     * assume info_out to be ll equivalent of info_in */
    else if (info_out->pj == NULL) {
        const char *projstr = NULL;
        char *indef = NULL;
 
        /* get PROJ-style definition */
        indef = G_store(info_in->def);
        G_debug(1, "ll equivalent definition: %s", indef);
 
        /* what about axis order?
         * is it always enu?
         * probably yes, as long as there is no +proj=axisswap step */
        G_asprintf(&(info_trans->def), "+proj=pipeline +step +inv %s", indef);
        info_trans->pj = proj_create(PJ_DEFAULT_CTX, info_trans->def);
        if (info_trans->pj == NULL) {
            G_warning(_("proj_create() failed for '%s'"), info_trans->def);
            G_free(indef);
 
            return -1;
        }
        projstr = proj_as_proj_string(NULL, info_trans->pj, PJ_PROJ_5, NULL);
        if (projstr == NULL) {
            G_warning(_("proj_create() failed for '%s'"), info_trans->def);
            G_free(indef);
 
            return -1;
        }
        G_free(indef);
    }
    /* input and output CRS are available */
    else if (info_in->def && info_out->pj && info_out->def) {
        char *indef = NULL, *outdef = NULL;
        char *insrid = NULL, *outsrid = NULL;
        PJ *in_pj, *out_pj;
        PJ_OBJ_LIST *op_list;
        PJ_OPERATION_FACTORY_CONTEXT *operation_ctx;
        PJ_AREA *pj_area = NULL;
        double xmin, xmax, ymin, ymax;
        int op_count = 0, op_count_area = 0;
 
        /* get pj_area */
        /* do it here because get_pj_area() will use
         * the PROJ definition for simple transformation to the
         * ll equivalent and we need to do unit conversion */
        if (get_pj_area(info_in, &xmin, &xmax, &ymin, &ymax)) {
            pj_area = proj_area_create();
            proj_area_set_bbox(pj_area, xmin, ymin, xmax, ymax);
        }
        else {
            G_warning(_("Unable to determine area of interest for '%s'"),
                      info_in->def);
 
            return -1;
        }
 
        G_debug(1, "source proj string: %s", info_in->def);
        G_debug(1, "source type: %s", get_pj_type_string(info_in->pj));
 
        /* PROJ6+: EPSG must be uppercase EPSG */
        if (info_in->srid) {
            if (strncmp(info_in->srid, "epsg", 4) == 0) {
                insrid = G_store_upper(info_in->srid);
                G_free(info_in->srid);
                ((struct pj_info *)info_in)->srid = insrid;
                insrid = NULL;
            }
        }
 
        in_pj = get_pj_object(info_in, &indef);
        if (in_pj == NULL || indef == NULL) {
            G_warning(_("Input CRS not available for '%s'"), info_in->def);
 
            return -1;
        }
        G_debug(1, "Input CRS definition: %s", indef);
 
        G_debug(1, "target proj string: %s", info_out->def);
        G_debug(1, "target type: %s", get_pj_type_string(info_out->pj));
 
        /* PROJ6+: EPSG must be uppercase EPSG */
        if (info_out->srid) {
            if (strncmp(info_out->srid, "epsg", 4) == 0) {
                outsrid = G_store_upper(info_out->srid);
                G_free(info_out->srid);
                ((struct pj_info *)info_out)->srid = outsrid;
                outsrid = NULL;
            }
        }
 
        out_pj = get_pj_object(info_out, &outdef);
        if (out_pj == NULL || outdef == NULL) {
            G_warning(_("Output CRS not available for '%s'"), info_out->def);
 
            return -1;
        }
        G_debug(1, "Output CRS definition: %s", outdef);
 
        /* check number of operations */
 
        operation_ctx =
            proj_create_operation_factory_context(PJ_DEFAULT_CTX, NULL);
        /* proj_create_operations() works only if both source_crs
         * and target_crs are found in the proj db
         * if any is not found, proj can not get a list of operations
         * and we have to take care of datumshift manually */
        /* list all operations irrespecitve of area and
         * grid availability */
        op_list = proj_create_operations(PJ_DEFAULT_CTX, in_pj, out_pj,
                                         operation_ctx);
        proj_operation_factory_context_destroy(operation_ctx);
 
        op_count = 0;
        if (op_list)
            op_count = proj_list_get_count(op_list);
        if (op_count > 1) {
            int i;
 
            G_important_message(_("Found %d possible transformations"),
                                op_count);
            for (i = 0; i < op_count; i++) {
                const char *area_of_use, *projstr;
                double e, w, s, n;
                PJ_PROJ_INFO pj_info;
                PJ *op, *op_norm;
 
                op = proj_list_get(PJ_DEFAULT_CTX, op_list, i);
                op_norm = proj_normalize_for_visualization(PJ_DEFAULT_CTX, op);
 
                if (!op_norm) {
                    G_warning(_("proj_normalize_for_visualization() failed for "
                                "operation %d"),
                              i + 1);
                }
                else {
                    proj_destroy(op);
                    op = op_norm;
                }
 
                pj_info = proj_pj_info(op);
                proj_get_area_of_use(NULL, op, &w, &s, &e, &n, &area_of_use);
                G_important_message("************************");
                G_important_message(_("Operation %d:"), i + 1);
                if (pj_info.description) {
                    G_important_message(_("Description: %s"),
                                        pj_info.description);
                }
                if (area_of_use) {
                    G_important_message(" ");
                    G_important_message(_("Area of use: %s"), area_of_use);
                }
                if (pj_info.accuracy > 0) {
                    G_important_message(" ");
                    G_important_message(_("Accuracy within area of use: %g m"),
                                        pj_info.accuracy);
                }
                const char *str = proj_get_remarks(op);
 
                if (str && *str) {
                    G_important_message(" ");
                    G_important_message(_("Remarks: %s"), str);
                }
                str = proj_get_scope(op);
                if (str && *str) {
                    G_important_message(" ");
                    G_important_message(_("Scope: %s"), str);
                }
 
                projstr = proj_as_proj_string(NULL, op, PJ_PROJ_5, NULL);
                if (projstr) {
                    G_important_message(" ");
                    G_important_message(_("PROJ string:"));
                    G_important_message("%s", projstr);
                }
                proj_destroy(op);
            }
            G_important_message("************************");
 
            G_important_message(_("See also output of:"));
            G_important_message("projinfo -o PROJ -s \"%s\" -t \"%s\"", indef,
                                outdef);
            G_important_message(_("Please provide the appropriate PROJ string "
                                  "with the %s option"),
                                "pipeline");
            G_important_message("************************");
        }
 
        if (op_list)
            proj_list_destroy(op_list);
 
        /* following code copied from proj_create_crs_to_crs_from_pj()
         * in proj src/4D_api.cpp
         * using PROJ_SPATIAL_CRITERION_PARTIAL_INTERSECTION
         * this can cause problems and artefacts
         * switch to PROJ_SPATIAL_CRITERION_STRICT_CONTAINMENT
         * in case of problems
         * but results can be different from gdalwarp:
         * shifted geolocation in some areas
         * in these cases there is no right or wrong,
         * different pipelines are all regarded as valid by PROJ
         * depending on the area of interest
         *
         * see also:
         * OGRProjCT::ListCoordinateOperations() in GDAL ogr/ogrct.cpp
         * create_operation_to_geog_crs() in PROJ src/4D_api.cpp
         * proj_create_crs_to_crs_from_pj() in PROJ src/4D_api.cpp
         * proj_operation_factory_context_set_spatial_criterion() in PROJ
         * src/iso19111/c_api.cpp
         *  */
 
        /* now use the current region as area of interest */
        operation_ctx =
            proj_create_operation_factory_context(PJ_DEFAULT_CTX, NULL);
        proj_operation_factory_context_set_area_of_interest(
            PJ_DEFAULT_CTX, operation_ctx, xmin, ymin, xmax, ymax);
        proj_operation_factory_context_set_spatial_criterion(
            PJ_DEFAULT_CTX, operation_ctx,
            PROJ_SPATIAL_CRITERION_PARTIAL_INTERSECTION);
        /* from GDAL OGRProjCT::ListCoordinateOperations() */
        proj_operation_factory_context_set_grid_availability_use(
            PJ_DEFAULT_CTX, operation_ctx,
            proj_context_is_network_enabled(PJ_DEFAULT_CTX)
                ? PROJ_GRID_AVAILABILITY_KNOWN_AVAILABLE
                : PROJ_GRID_AVAILABILITY_DISCARD_OPERATION_IF_MISSING_GRID);
 
        /* The operations are sorted with the most relevant ones first:
         * by descending area (intersection of the transformation area
         * with the area of interest, or intersection of the
         * transformation with the area of use of the CRS),
         * and by increasing accuracy.
         * Operations with unknown accuracy are sorted last,
         * whatever their area.
         */
        op_list = proj_create_operations(PJ_DEFAULT_CTX, in_pj, out_pj,
                                         operation_ctx);
        proj_operation_factory_context_destroy(operation_ctx);
        op_count_area = 0;
        if (op_list)
            op_count_area = proj_list_get_count(op_list);
        if (op_count_area == 0) {
            /* no operations */
            info_trans->pj = NULL;
        }
        else if (op_count_area == 1) {
            info_trans->pj = proj_list_get(PJ_DEFAULT_CTX, op_list, 0);
        }
        else { /* op_count_area > 1 */
            /* can't use pj_create_prepared_operations()
             * this is a PROJ-internal function
             * trust the sorting of PROJ and use the first one */
            info_trans->pj = proj_list_get(PJ_DEFAULT_CTX, op_list, 0);
        }
        if (op_list)
            proj_list_destroy(op_list);
 
        /* try proj_create_crs_to_crs() */
        /*
           G_debug(1, "trying %s to %s", indef, outdef);
         */
 
        /* proj_create_crs_to_crs() does not work because it calls
         * proj_create_crs_to_crs_from_pj() which calls
         * proj_operation_factory_context_set_spatial_criterion()
         * with PROJ_SPATIAL_CRITERION_PARTIAL_INTERSECTION
         * instead of
         * PROJ_SPATIAL_CRITERION_STRICT_CONTAINMENT
         *
         * fixed in PROJ master, probably available with PROJ 7.3.x */
 
        /*
           info_trans->pj = proj_create_crs_to_crs(PJ_DEFAULT_CTX,
           indef,
           outdef,
           pj_area);
         */
 
        if (in_pj)
            proj_destroy(in_pj);
        if (out_pj)
            proj_destroy(out_pj);
 
        if (info_trans->pj) {
            const char *projstr;
            PJ *pj_norm = NULL;
 
            G_debug(1, "proj_create_crs_to_crs() succeeded with PROJ%d",
                    PROJ_VERSION_MAJOR);
 
            projstr =
                proj_as_proj_string(NULL, info_trans->pj, PJ_PROJ_5, NULL);
 
            info_trans->def = G_store(projstr);
 
            if (projstr) {
                /* make sure axis order is easting, northing
                 * proj_normalize_for_visualization() requires
                 * source and target CRS
                 * -> does not work with ll equivalent of input:
                 * no target CRS in +proj=pipeline +step +inv %s */
                pj_norm = proj_normalize_for_visualization(PJ_DEFAULT_CTX,
                                                           info_trans->pj);
 
                if (!pj_norm) {
                    G_warning(
                        _("proj_normalize_for_visualization() failed for '%s'"),
                        info_trans->def);
                }
                else {
                    projstr =
                        proj_as_proj_string(NULL, pj_norm, PJ_PROJ_5, NULL);
                    if (projstr && *projstr) {
                        proj_destroy(info_trans->pj);
                        info_trans->pj = pj_norm;
                        info_trans->def = G_store(projstr);
                    }
                    else {
                        proj_destroy(pj_norm);
                        G_warning(_("No PROJ definition for normalized version "
                                    "of '%s'"),
                                  info_trans->def);
                    }
                }
                G_important_message(_("Selected PROJ pipeline:"));
                G_important_message(_("%s"), info_trans->def);
                G_important_message("************************");
            }
            else {
                proj_destroy(info_trans->pj);
                info_trans->pj = NULL;
            }
        }
 
        if (pj_area)
            proj_area_destroy(pj_area);
 
        if (insrid)
            G_free(insrid);
        if (outsrid)
            G_free(outsrid);
        G_free(indef);
        G_free(outdef);
    }
    if (info_trans->pj == NULL) {
        G_warning(_("proj_create() failed for '%s'"), info_trans->def);
 
        return -1;
    }
 
    return 1;
}
 
/* TODO: rename pj_ to GPJ_ to avoid symbol clash with PROJ lib */
 
/**
 * \brief Re-project a point between two co-ordinate systems using a
 *        transformation object prepared with GPJ_prepare_pj()
 *
 * This function takes pointers to three pj_info structures as arguments,
 * and projects a point between the input and output co-ordinate system.
 * The pj_info structure info_trans must have been initialized with
 * GPJ_init_transform().
 * The direction determines if a point is projected from input CRS to
 * output CRS (PJ_FWD) or from output CRS to input CRS (PJ_INV).
 * The easting, northing, and height of the point are contained in the
 * pointers passed to the function; these will be overwritten by the
 * coordinates of the transformed point.
 *
 * \param info_in pointer to pj_info struct for input co-ordinate system
 * \param info_out pointer to pj_info struct for output co-ordinate system
 * \param info_trans pointer to pj_info struct for a transformation object (PROJ
 *5+) \param dir direction of the transformation (PJ_FWD or PJ_INV) \param x
 *Pointer to a double containing easting or longitude \param y Pointer to a
 *double containing northing or latitude \param z Pointer to a double containing
 *height, or NULL
 *
 * \return Return value from PROJ proj_trans() function
 **/
 
int GPJ_transform(const struct pj_info *info_in, const struct pj_info *info_out,
                  const struct pj_info *info_trans, int dir, double *x,
                  double *y, double *z)
{
    int ok = 0;
 
    int in_is_ll, out_is_ll, in_deg2rad, out_rad2deg;
    PJ_COORD c;
 
    if (info_in->pj == NULL)
        G_fatal_error(_("No input projection"));
 
    if (info_trans->pj == NULL)
        G_fatal_error(_("No transformation object"));
 
    in_deg2rad = out_rad2deg = 1;
    if (dir == PJ_FWD) {
        /* info_in -> info_out */
        METERS_in = info_in->meters;
        in_is_ll = !strncmp(info_in->proj, "ll", 2);
        /* PROJ 6+: conversion to radians is not always needed:
         * if proj_angular_input(info_trans->pj, dir) == 1
         * -> convert from degrees to radians */
        if (in_is_ll && proj_angular_input(info_trans->pj, dir) == 0) {
            in_deg2rad = 0;
        }
        if (info_out->pj) {
            METERS_out = info_out->meters;
            out_is_ll = !strncmp(info_out->proj, "ll", 2);
            /* PROJ 6+: conversion to radians is not always needed:
             * if proj_angular_input(info_trans->pj, dir) == 1
             * -> convert from degrees to radians */
            if (out_is_ll && proj_angular_output(info_trans->pj, dir) == 0) {
                out_rad2deg = 0;
            }
        }
        else {
            METERS_out = 1.0;
            out_is_ll = 1;
        }
    }
    else {
        /* info_out -> info_in */
        METERS_out = info_in->meters;
        out_is_ll = !strncmp(info_in->proj, "ll", 2);
        /* PROJ 6+: conversion to radians is not always needed:
         * if proj_angular_input(info_trans->pj, dir) == 1
         * -> convert from degrees to radians */
        if (out_is_ll && proj_angular_output(info_trans->pj, dir) == 0) {
            out_rad2deg = 0;
        }
        if (info_out->pj) {
            METERS_in = info_out->meters;
            in_is_ll = !strncmp(info_out->proj, "ll", 2);
            /* PROJ 6+: conversion to radians is not always needed:
             * if proj_angular_input(info_trans->pj, dir) == 1
             * -> convert from degrees to radians */
            if (in_is_ll && proj_angular_input(info_trans->pj, dir) == 0) {
                in_deg2rad = 0;
            }
        }
        else {
            METERS_in = 1.0;
            in_is_ll = 1;
        }
    }
 
    /* prepare */
    if (in_is_ll) {
        if (in_deg2rad) {
            /* convert degrees to radians */
            c.lpzt.lam = (*x) / RAD_TO_DEG;
            c.lpzt.phi = (*y) / RAD_TO_DEG;
        }
        else {
            c.lpzt.lam = (*x);
            c.lpzt.phi = (*y);
        }
        c.lpzt.z = 0;
        if (z)
            c.lpzt.z = *z;
        c.lpzt.t = 0;
    }
    else {
        /* convert to meters */
        c.xyzt.x = *x * METERS_in;
        c.xyzt.y = *y * METERS_in;
        c.xyzt.z = 0;
        if (z)
            c.xyzt.z = *z;
        c.xyzt.t = 0;
    }
 
    G_debug(1, "c.xyzt.x: %g", c.xyzt.x);
    G_debug(1, "c.xyzt.y: %g", c.xyzt.y);
    G_debug(1, "c.xyzt.z: %g", c.xyzt.z);
 
    /* transform */
    c = proj_trans(info_trans->pj, dir, c);
    ok = proj_errno(info_trans->pj);
 
    if (ok < 0) {
        G_warning(_("proj_trans() failed: %s"), proj_errno_string(ok));
        return ok;
    }
 
    /* output */
    if (out_is_ll) {
        /* convert to degrees */
        if (out_rad2deg) {
            /* convert radians to degrees */
            *x = c.lp.lam * RAD_TO_DEG;
            *y = c.lp.phi * RAD_TO_DEG;
        }
        else {
            *x = c.lp.lam;
            *y = c.lp.phi;
        }
        if (z)
            *z = c.lpzt.z;
    }
    else {
        /* convert to map units */
        *x = c.xyzt.x / METERS_out;
        *y = c.xyzt.y / METERS_out;
        if (z)
            *z = c.xyzt.z;
    }
 
    return ok;
}
 
/**
 * \brief Re-project an array of points between two co-ordinate systems
 *        using a transformation object prepared with GPJ_prepare_pj()
 *
 * This function takes pointers to three pj_info structures as arguments,
 * and projects an array of pointd between the input and output
 * co-ordinate system. The pj_info structure info_trans must have been
 * initialized with GPJ_init_transform().
 * The direction determines if a point is projected from input CRS to
 * output CRS (PJ_FWD) or from output CRS to input CRS (PJ_INV).
 * The easting, northing, and height of the point are contained in the
 * pointers passed to the function; these will be overwritten by the
 * coordinates of the transformed point.
 *
 * \param info_in pointer to pj_info struct for input co-ordinate system
 * \param info_out pointer to pj_info struct for output co-ordinate system
 * \param info_trans pointer to pj_info struct for a transformation object (PROJ
 *5+) \param dir direction of the transformation (PJ_FWD or PJ_INV) \param x
 *pointer to an array of type double containing easting or longitude \param y
 *pointer to an array of type double containing northing or latitude \param z
 *pointer to an array of type double containing height, or NULL \param n number
 *of points in the arrays to be transformed
 *
 * \return Return value from PROJ proj_trans() function
 **/
 
int GPJ_transform_array(const struct pj_info *info_in,
                        const struct pj_info *info_out,
                        const struct pj_info *info_trans, int dir, double *x,
                        double *y, double *z, int n)
{
    int ok;
    int i;
    int has_z = 1;
 
    /* PROJ 5+ variant */
    int in_is_ll, out_is_ll, in_deg2rad, out_rad2deg;
    PJ_COORD c;
 
    if (info_trans->pj == NULL)
        G_fatal_error(_("No transformation object"));
 
    in_deg2rad = out_rad2deg = 1;
    if (dir == PJ_FWD) {
        /* info_in -> info_out */
        METERS_in = info_in->meters;
        in_is_ll = !strncmp(info_in->proj, "ll", 2);
        /* PROJ 6+: conversion to radians is not always needed:
         * if proj_angular_input(info_trans->pj, dir) == 1
         * -> convert from degrees to radians */
        if (in_is_ll && proj_angular_input(info_trans->pj, dir) == 0) {
            in_deg2rad = 0;
        }
        if (info_out->pj) {
            METERS_out = info_out->meters;
            out_is_ll = !strncmp(info_out->proj, "ll", 2);
            /* PROJ 6+: conversion to radians is not always needed:
             * if proj_angular_input(info_trans->pj, dir) == 1
             * -> convert from degrees to radians */
            if (out_is_ll && proj_angular_output(info_trans->pj, dir) == 0) {
                out_rad2deg = 0;
            }
        }
        else {
            METERS_out = 1.0;
            out_is_ll = 1;
        }
    }
    else {
        /* info_out -> info_in */
        METERS_out = info_in->meters;
        out_is_ll = !strncmp(info_in->proj, "ll", 2);
        /* PROJ 6+: conversion to radians is not always needed:
         * if proj_angular_input(info_trans->pj, dir) == 1
         * -> convert from degrees to radians */
        if (out_is_ll && proj_angular_output(info_trans->pj, dir) == 0) {
            out_rad2deg = 0;
        }
        if (info_out->pj) {
            METERS_in = info_out->meters;
            in_is_ll = !strncmp(info_out->proj, "ll", 2);
            /* PROJ 6+: conversion to degrees is not always needed:
             * if proj_angular_output(info_trans->pj, dir) == 1
             * -> convert from degrees to radians */
            if (in_is_ll && proj_angular_input(info_trans->pj, dir) == 0) {
                in_deg2rad = 0;
            }
        }
        else {
            METERS_in = 1.0;
            in_is_ll = 1;
        }
    }
 
    if (z == NULL) {
        z = G_malloc(sizeof(double) * n);
        /* they say memset is only guaranteed for chars ;-( */
        for (i = 0; i < n; i++)
            z[i] = 0.0;
        has_z = 0;
    }
    ok = 0;
    if (in_is_ll) {
        c.lpzt.t = 0;
        if (out_is_ll) {
            /* what is more costly ?
             * calling proj_trans for each point
             * or having three loops over all points ?
             * proj_trans_array() itself calls proj_trans() in a loop
             * -> one loop over all points is better than
             *    three loops over all points
             */
            for (i = 0; i < n; i++) {
                if (in_deg2rad) {
                    /* convert degrees to radians */
                    c.lpzt.lam = (*x) / RAD_TO_DEG;
                    c.lpzt.phi = (*y) / RAD_TO_DEG;
                }
                else {
                    c.lpzt.lam = (*x);
                    c.lpzt.phi = (*y);
                }
                c.lpzt.z = z[i];
                c = proj_trans(info_trans->pj, dir, c);
                if ((ok = proj_errno(info_trans->pj)) < 0)
                    break;
                if (out_rad2deg) {
                    /* convert radians to degrees */
                    x[i] = c.lp.lam * RAD_TO_DEG;
                    y[i] = c.lp.phi * RAD_TO_DEG;
                }
                else {
                    x[i] = c.lp.lam;
                    y[i] = c.lp.phi;
                }
            }
        }
        else {
            for (i = 0; i < n; i++) {
                if (in_deg2rad) {
                    /* convert degrees to radians */
                    c.lpzt.lam = (*x) / RAD_TO_DEG;
                    c.lpzt.phi = (*y) / RAD_TO_DEG;
                }
                else {
                    c.lpzt.lam = (*x);
                    c.lpzt.phi = (*y);
                }
                c.lpzt.z = z[i];
                c = proj_trans(info_trans->pj, dir, c);
                if ((ok = proj_errno(info_trans->pj)) < 0)
                    break;
 
                /* convert to map units */
                x[i] = c.xy.x / METERS_out;
                y[i] = c.xy.y / METERS_out;
            }
        }
    }
    else {
        c.xyzt.t = 0;
        if (out_is_ll) {
            for (i = 0; i < n; i++) {
                /* convert to meters */
                c.xyzt.x = x[i] * METERS_in;
                c.xyzt.y = y[i] * METERS_in;
                c.xyzt.z = z[i];
                c = proj_trans(info_trans->pj, dir, c);
                if ((ok = proj_errno(info_trans->pj)) < 0)
                    break;
                if (out_rad2deg) {
                    /* convert radians to degrees */
                    x[i] = c.lp.lam * RAD_TO_DEG;
                    y[i] = c.lp.phi * RAD_TO_DEG;
                }
                else {
                    x[i] = c.lp.lam;
                    y[i] = c.lp.phi;
                }
            }
        }
        else {
            for (i = 0; i < n; i++) {
                /* convert to meters */
                c.xyzt.x = x[i] * METERS_in;
                c.xyzt.y = y[i] * METERS_in;
                c.xyzt.z = z[i];
                c = proj_trans(info_trans->pj, dir, c);
                if ((ok = proj_errno(info_trans->pj)) < 0)
                    break;
                /* convert to map units */
                x[i] = c.xy.x / METERS_out;
                y[i] = c.xy.y / METERS_out;
            }
        }
    }
    if (!has_z)
        G_free(z);
 
    if (ok < 0) {
        G_warning(_("proj_trans() failed: %s"), proj_errno_string(ok));
    }
 
    return ok;
}
 
/*
 * old API, to be deleted
 */
 
/**
 * \brief Re-project a point between two co-ordinate systems
 *
 * This function takes pointers to two pj_info structures as arguments,
 * and projects a point between the co-ordinate systems represented by them.
 * The easting and northing of the point are contained in two pointers passed
 * to the function; these will be overwritten by the co-ordinates of the
 * re-projected point.
 *
 * \param x Pointer to a double containing easting or longitude
 * \param y Pointer to a double containing northing or latitude
 * \param info_in pointer to pj_info struct for input co-ordinate system
 * \param info_out pointer to pj_info struct for output co-ordinate system
 *
 * \return Return value from PROJ proj_trans() function
 **/
 
int pj_do_proj(double *x, double *y, const struct pj_info *info_in,
               const struct pj_info *info_out)
{
    int ok;
 
    struct pj_info info_trans;
    PJ_COORD c;
 
    if (GPJ_init_transform(info_in, info_out, &info_trans) < 0) {
        return -1;
    }
 
    METERS_in = info_in->meters;
    METERS_out = info_out->meters;
 
    if (strncmp(info_in->proj, "ll", 2) == 0) {
        /* convert to radians */
        c.lpzt.lam = (*x) / RAD_TO_DEG;
        c.lpzt.phi = (*y) / RAD_TO_DEG;
        c.lpzt.z = 0;
        c.lpzt.t = 0;
        c = proj_trans(info_trans.pj, PJ_FWD, c);
        ok = proj_errno(info_trans.pj);
 
        if (strncmp(info_out->proj, "ll", 2) == 0) {
            /* convert to degrees */
            *x = c.lp.lam * RAD_TO_DEG;
            *y = c.lp.phi * RAD_TO_DEG;
        }
        else {
            /* convert to map units */
            *x = c.xy.x / METERS_out;
            *y = c.xy.y / METERS_out;
        }
    }
    else {
        /* convert to meters */
        c.xyzt.x = *x * METERS_in;
        c.xyzt.y = *y * METERS_in;
        c.xyzt.z = 0;
        c.xyzt.t = 0;
        c = proj_trans(info_trans.pj, PJ_FWD, c);
        ok = proj_errno(info_trans.pj);
 
        if (strncmp(info_out->proj, "ll", 2) == 0) {
            /* convert to degrees */
            *x = c.lp.lam * RAD_TO_DEG;
            *y = c.lp.phi * RAD_TO_DEG;
        }
        else {
            /* convert to map units */
            *x = c.xy.x / METERS_out;
            *y = c.xy.y / METERS_out;
        }
    }
    proj_destroy(info_trans.pj);
 
    if (ok < 0) {
        G_warning(_("proj_trans() failed: %d"), ok);
    }
    return ok;
}
 
/**
 * \brief Re-project an array of points between two co-ordinate systems with
 *        optional ellipsoidal height conversion
 *
 * This function takes pointers to two pj_info structures as arguments,
 * and projects an array of points between the co-ordinate systems
 * represented by them. Pointers to the three arrays of easting, northing,
 * and ellipsoidal height of the point (this one may be NULL) are passed
 * to the function; these will be overwritten by the co-ordinates of the
 * re-projected points.
 *
 * \param count Number of points in the arrays to be transformed
 * \param x Pointer to an array of type double containing easting or longitude
 * \param y Pointer to an array of type double containing northing or latitude
 * \param h Pointer to an array of type double containing ellipsoidal height.
 *          May be null in which case a two-dimensional re-projection will be
 *          done
 * \param info_in pointer to pj_info struct for input co-ordinate system
 * \param info_out pointer to pj_info struct for output co-ordinate system
 *
 * \return Return value from PROJ proj_trans() function
 **/
 
int pj_do_transform(int count, double *x, double *y, double *h,
                    const struct pj_info *info_in,
                    const struct pj_info *info_out)
{
    int ok;
    int i;
    int has_h = 1;
 
    struct pj_info info_trans;
    PJ_COORD c;
 
    if (GPJ_init_transform(info_in, info_out, &info_trans) < 0) {
        return -1;
    }
 
    METERS_in = info_in->meters;
    METERS_out = info_out->meters;
 
    if (h == NULL) {
        h = G_malloc(sizeof *h * count);
        /* they say memset is only guaranteed for chars ;-( */
        for (i = 0; i < count; ++i)
            h[i] = 0.0;
        has_h = 0;
    }
    ok = 0;
    if (strncmp(info_in->proj, "ll", 2) == 0) {
        c.lpzt.t = 0;
        if (strncmp(info_out->proj, "ll", 2) == 0) {
            for (i = 0; i < count; i++) {
                /* convert to radians */
                c.lpzt.lam = x[i] / RAD_TO_DEG;
                c.lpzt.phi = y[i] / RAD_TO_DEG;
                c.lpzt.z = h[i];
                c = proj_trans(info_trans.pj, PJ_FWD, c);
                if ((ok = proj_errno(info_trans.pj)) < 0)
                    break;
                /* convert to degrees */
                x[i] = c.lp.lam * RAD_TO_DEG;
                y[i] = c.lp.phi * RAD_TO_DEG;
            }
        }
        else {
            for (i = 0; i < count; i++) {
                /* convert to radians */
                c.lpzt.lam = x[i] / RAD_TO_DEG;
                c.lpzt.phi = y[i] / RAD_TO_DEG;
                c.lpzt.z = h[i];
                c = proj_trans(info_trans.pj, PJ_FWD, c);
                if ((ok = proj_errno(info_trans.pj)) < 0)
                    break;
                /* convert to map units */
                x[i] = c.xy.x / METERS_out;
                y[i] = c.xy.y / METERS_out;
            }
        }
    }
    else {
        c.xyzt.t = 0;
        if (strncmp(info_out->proj, "ll", 2) == 0) {
            for (i = 0; i < count; i++) {
                /* convert to meters */
                c.xyzt.x = x[i] * METERS_in;
                c.xyzt.y = y[i] * METERS_in;
                c.xyzt.z = h[i];
                c = proj_trans(info_trans.pj, PJ_FWD, c);
                if ((ok = proj_errno(info_trans.pj)) < 0)
                    break;
                /* convert to degrees */
                x[i] = c.lp.lam * RAD_TO_DEG;
                y[i] = c.lp.phi * RAD_TO_DEG;
            }
        }
        else {
            for (i = 0; i < count; i++) {
                /* convert to meters */
                c.xyzt.x = x[i] * METERS_in;
                c.xyzt.y = y[i] * METERS_in;
                c.xyzt.z = h[i];
                c = proj_trans(info_trans.pj, PJ_FWD, c);
                if ((ok = proj_errno(info_trans.pj)) < 0)
                    break;
                /* convert to map units */
                x[i] = c.xy.x / METERS_out;
                y[i] = c.xy.y / METERS_out;
            }
        }
    }
    if (!has_h)
        G_free(h);
    proj_destroy(info_trans.pj);
 
    if (ok < 0) {
        G_warning(_("proj_trans() failed: %d"), ok);
    }
    return ok;
}