EVOLUTION-MANAGER

Edit File: apply_gridshift.cpp

/******************************************************************************
 * Project:  PROJ.4
 * Purpose:  Apply datum shifts based on grid shift files (normally NAD27 to
 *           NAD83 or the reverse).  This module is responsible for keeping
 *           a list of loaded grids, and calling with each one that is
 *           allowed for a given datum (expressed as the nadgrids= parameter).
 * Author:   Frank Warmerdam, warmerdam@pobox.com
 *
 ******************************************************************************
 * Copyright (c) 2000, Frank Warmerdam <warmerdam@pobox.com>
 *
 * Permission is hereby granted, free of charge, to any person obtaining a
 * copy of this software and associated documentation files (the "Software"),
 * to deal in the Software without restriction, including without limitation
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 * and/or sell copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included
 * in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 * DEALINGS IN THE SOFTWARE.
 *****************************************************************************/

#define PJ_LIB__

#include <math.h>
#include <stdio.h>
#include <string.h>

#include "proj.h"
#include "proj_internal.h"

/************************************************************************/
/*                         pj_apply_gridshift()                         */
/*                                                                      */
/*      This is the externally callable interface - part of the         */
/*      public API - though it is not used internally any more and I    */
/*      doubt it is used by any other applications.  But we preserve    */
/*      it to honour our public api.                                    */
/************************************************************************/

int pj_apply_gridshift( projCtx ctx, const char *nadgrids, int inverse,
                        long point_count, int point_offset,
                        double *x, double *y, double *z )

{
    PJ_GRIDINFO **gridlist;
    int           grid_count;
    int           ret;

gridlist = pj_gridlist_from_nadgrids( ctx, nadgrids, &grid_count );

if( gridlist == nullptr || grid_count == 0 )
    {
        pj_dalloc( gridlist );
        return ctx->last_errno;
    }

ret = pj_apply_gridshift_3( ctx, gridlist, grid_count, inverse,
                                point_count, point_offset, x, y, z );

/*
    ** Note this frees the array of grid list pointers, but not the grids
    ** which is as intended.  The grids themselves live on.
    */
    pj_dalloc( gridlist );

return ret;
}

/************************************************************************/
/*                        pj_apply_gridshift_2()                        */
/*                                                                      */
/*      This implementation uses the gridlist from a coordinate         */
/*      system definition.  If the gridlist has not yet been            */
/*      populated in the coordinate system definition we set it up      */
/*      now.                                                            */
/************************************************************************/

int pj_apply_gridshift_2( PJ *defn, int inverse,
                          long point_count, int point_offset,
                          double *x, double *y, double *z )

{
    if( defn->catalog_name != nullptr )
        return pj_gc_apply_gridshift( defn, inverse, point_count, point_offset,
                                      x, y, z );

if( defn->gridlist == nullptr )
    {
        defn->gridlist =
            pj_gridlist_from_nadgrids( pj_get_ctx( defn ),
                                       pj_param(defn->ctx, defn->params,"snadgrids").s,
                                       &(defn->gridlist_count) );

if( defn->gridlist == nullptr || defn->gridlist_count == 0 )
            return defn->ctx->last_errno;
    }

return pj_apply_gridshift_3( pj_get_ctx( defn ),
                                 defn->gridlist, defn->gridlist_count, inverse,
                                 point_count, point_offset, x, y, z );
}

/************************************************************************/
/*                             find_ctable()                            */
/*                                                                      */
/*    Determine which grid is the correct given an input coordinate.    */
/************************************************************************/

struct CTABLE* find_ctable(projCtx ctx, PJ_LP input, int grid_count, PJ_GRIDINFO **tables) {
    int itable;

/* keep trying till we find a table that works */
    for( itable = 0; itable < grid_count; itable++ )
    {

PJ_GRIDINFO *gi = tables[itable];
        struct CTABLE *ct = gi->ct;
        double epsilon = (fabs(ct->del.phi)+fabs(ct->del.lam))/10000.0;
        /* skip tables that don't match our point at all.  */
        if ( ct->ll.phi - epsilon > input.phi
                || ct->ll.lam - epsilon > input.lam
                || (ct->ll.phi + (ct->lim.phi-1) * ct->del.phi + epsilon < input.phi)
                || (ct->ll.lam + (ct->lim.lam-1) * ct->del.lam + epsilon < input.lam) ) {
            continue;
        }

/* If we have child nodes, check to see if any of them apply. */
        while( gi->child )
        {
            PJ_GRIDINFO *child;

for( child = gi->child; child != nullptr; child = child->next )
            {
                struct CTABLE *ct1 = child->ct;
                epsilon = (fabs(ct1->del.phi)+fabs(ct1->del.lam))/10000.0;

if( ct1->ll.phi - epsilon > input.phi
                    || ct1->ll.lam - epsilon > input.lam
                    || (ct1->ll.phi+(ct1->lim.phi-1)*ct1->del.phi + epsilon < input.phi)
                    || (ct1->ll.lam+(ct1->lim.lam-1)*ct1->del.lam + epsilon < input.lam) ) {
                    continue;
                }
                break;
            }

/* If we didn't find a child then nothing more to do */
            if( child == nullptr ) break;

/* Otherwise use the child, first checking it's children */
            gi = child;
            ct = child->ct;
        }
        /* load the grid shift info if we don't have it. */
        if( ct->cvs == nullptr) {
            if (!pj_gridinfo_load( ctx, gi ) ) {
                pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
                return nullptr;
            }
        }
        /* if we get this far we have found a suitable grid */
        return ct;
    }

return nullptr;
}

/************************************************************************/
/*                        pj_apply_gridshift_3()                        */
/*                                                                      */
/*      This is the real workhorse, given a gridlist.                   */
/************************************************************************/

int pj_apply_gridshift_3( projCtx ctx, PJ_GRIDINFO **gridlist, int gridlist_count,
                          int inverse, long point_count, int point_offset,
                          double *x, double *y, double *z )
{
    int  i;
    struct CTABLE *ct;
    static int debug_count = 0;
    (void) z;

if( gridlist== nullptr || gridlist_count == 0 )
    {
        pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
        return PJD_ERR_FAILED_TO_LOAD_GRID;
    }

ctx->last_errno = 0;

for( i = 0; i < point_count; i++ )
    {
        long io = i * point_offset;
        PJ_LP   input, output;
        int  itable;

input.phi = y[io];
        input.lam = x[io];
        output.phi = HUGE_VAL;
        output.lam = HUGE_VAL;

ct = find_ctable(ctx, input, gridlist_count, gridlist);
        if( ct != nullptr )
        {
            output = nad_cvt( ctx, input, inverse, ct, gridlist_count, gridlist);

if ( output.lam != HUGE_VAL && debug_count++ < 20 )
                pj_log( ctx, PJ_LOG_DEBUG_MINOR, "pj_apply_gridshift(): used %s", ct->id );
        }

if ( output.lam == HUGE_VAL )
        {
            if( ctx->debug_level >= PJ_LOG_DEBUG_MAJOR )
            {
                pj_log( ctx, PJ_LOG_DEBUG_MAJOR,
                    "pj_apply_gridshift(): failed to find a grid shift table for\n"
                    "                      location (%.7fdW,%.7fdN)",
                    x[io] * RAD_TO_DEG,
                    y[io] * RAD_TO_DEG );
                for( itable = 0; itable < gridlist_count; itable++ )
                {
                    PJ_GRIDINFO *gi = gridlist[itable];
                    if( itable == 0 )
                        pj_log( ctx, PJ_LOG_DEBUG_MAJOR, "   tried: %s", gi->gridname );
                    else
                        pj_log( ctx, PJ_LOG_DEBUG_MAJOR, ",%s", gi->gridname );
                }
            }

/*
             * We don't actually have any machinery currently to set the
             * following macro, so this is mostly kept here to make it clear
             * how we ought to operate if we wanted to make it super clear
             * that an error has occurred when points are outside our available
             * datum shift areas.  But if this is on, we will find that "low
             * value" points on the fringes of some datasets will completely
             * fail causing lots of problems when it is more or less ok to
             * just not apply a datum shift.  So rather than deal with
             * that we just fallback to no shift. (see also bug #45).
             */
#ifdef ERR_GRID_AREA_TRANSIENT_SEVERE
            y[io] = HUGE_VAL;
            x[io] = HUGE_VAL;
#else
            /* leave x/y unshifted. */
#endif
        }
        else
        {
            y[io] = output.phi;
            x[io] = output.lam;
        }
    }

return 0;
}

/**********************************************/
int proj_hgrid_init(PJ* P, const char *grids) {
/**********************************************

Initizalize and populate list of horizontal
  grids.

Takes a PJ-object and the plus-parameter
    name that is used in the proj-string to
    specify the grids to load, e.g. "+grids".
    The + should be left out here.

Returns the number of loaded grids.

***********************************************/

/* prepend "s" to the "grids" string to allow usage with pj_param */
    char *sgrids = (char *) pj_malloc( (strlen(grids)+1+1) *sizeof(char) );
    sprintf(sgrids, "%s%s", "s", grids);

if (P->gridlist == nullptr) {
        P->gridlist = pj_gridlist_from_nadgrids(
            P->ctx,
            pj_param(P->ctx, P->params, sgrids).s,
            &(P->gridlist_count)
        );

if( P->gridlist == nullptr || P->gridlist_count == 0 ) {
            pj_dealloc(sgrids);
            return 0;
        }
    }

if (P->gridlist_count == 0) {
        proj_errno_set(P, PJD_ERR_FAILED_TO_LOAD_GRID);
    }

pj_dealloc(sgrids);
    return P->gridlist_count;
}

/********************************************/
/*           proj_hgrid_value()             */
/*                                          */
/*    Return coordinate offset in grid      */
/********************************************/
PJ_LP proj_hgrid_value(PJ *P, PJ_LP lp) {
    struct CTABLE *ct;
    PJ_LP out = proj_coord_error().lp;

ct = find_ctable(P->ctx, lp, P->gridlist_count, P->gridlist);
    if (ct == nullptr) {
        pj_ctx_set_errno( P->ctx, PJD_ERR_GRID_AREA);
        return out;
    }

/* normalize input to ll origin */
    lp.lam -= ct->ll.lam;
    lp.phi -= ct->ll.phi;

lp.lam = adjlon(lp.lam - M_PI) + M_PI;

out = nad_intr(lp, ct);

if (out.lam == HUGE_VAL || out.phi == HUGE_VAL) {
        pj_ctx_set_errno(P->ctx, PJD_ERR_GRID_AREA);
    }

return out;
}

PJ_LP proj_hgrid_apply(PJ *P, PJ_LP lp, PJ_DIRECTION direction) {
    struct CTABLE *ct;
    int inverse;
    PJ_LP out;

out.lam = HUGE_VAL; out.phi = HUGE_VAL;

ct = find_ctable(P->ctx, lp, P->gridlist_count, P->gridlist);

if (ct == nullptr || ct->cvs == nullptr) {
        if( P->gridlist_count == 1 &&
            strcmp(P->gridlist[0]->gridname, "null") == 0) {
          // TODO: remove this particular case that is put there just to be
          // able to handle longitudes outside of -180,180
          out = lp;
        } else {
          pj_ctx_set_errno( P->ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
        }
        return out;
    }

inverse = direction == PJ_FWD ? 0 : 1;
    out = nad_cvt(P->ctx, lp, inverse, ct, P->gridlist_count, P->gridlist);

if (out.lam == HUGE_VAL || out.phi == HUGE_VAL)
        pj_ctx_set_errno(P->ctx, PJD_ERR_GRID_AREA);

return out;

}