Crash when reading NTV2 grid files containing subgrids not referenced to top level

[bpayne]

Attempting to read the attached ntv2 datum shift file (sk27-98.gsb) crashes proj4. This is because one of the subgrids is attached to a subgrid that is already a child of the top level grid, instead of the top level grid itself. (Well, actually it is the attempt to recover from this condition that causes the crash, but even if it did not crash, the intended logic did not handle this condition properly). This is a legal interpretation of the format specification.

Also, once the header data structure is read in correctly, there is an additional issue finding and subgrids of subgrids. This would be consistent with the original code not handling this general case. Code for this fix is also provided in this bug report.

To test positioning error after the initial crash, tranform point Lat/Lon? : 50.754, -103.680 To UTM zone 13 north, using the datum shift in sk27-98.gsb. The correct result, from subgrid SKfqsub is: 593150.462, 5623077.968 Proj 4 will give the result from subgrid SKfortQ: 593150.682, 5623078.040

I have also attached the following bits of code for your convenience.

The following snippet of code will update proj4 to allow the datum file to be read in correctly. This code will prevent the crash and create the datum structure accurately, but it can still be read incorrectly without the second piece.

///////////////////////////////////////////////////////////////////////////////// Code snippet from pj_gridinfo.c, function pj_gridinfo_init_ntv2, approx start line 513.

/* -------------------------------------------------------------------- */
/*      Attach to the correct list or sublist.                          */
/* -------------------------------------------------------------------- */
        if( strncmp((const char *)header+24,"NONE",4) == 0 )
        {
            if( gi != gilist )
            {
                PJ_GRIDINFO *lnk;

                for( lnk = gilist; lnk->next != NULL; lnk = lnk->next ) {}
                lnk->next = gi;
            }
        }

        else
        {
            PJ_GRIDINFO *lnk;
            PJ_GRIDINFO *gpChild;
            PJ_GRIDINFO *gp;
            
            for (gpChild = gilist; gpChild != NULL; gpChild = gpChild->child )
            {
                gp = gpChild;
                while( gp != NULL 
                       && strncmp(gp->ct->id,(const char*)header+24,8) != 0 )
                    gp = gp->next;

                if( gp == NULL )
                {
                    // Grid header not found, search the next level down
                    continue;
                }
                else if( gp->child == NULL )
                {
                    gp->child = gi;
                    break;
                }
                else
                {
                    for( lnk = gp->child; lnk->next != NULL; lnk = lnk->next ) {}
                    lnk->next = gi;
                    break;
                }
            }

            if (gpChild == NULL)
            {
                // Didn't find an attachment point for this grid, stick it at the end of the list
                pj_log( ctx, PJ_LOG_ERROR,
                    "pj_gridinfo_init_ntv2(): "
                    "failed to find parent %8.8s for %s.\n", 
                    (const char *) header+24, gi->ct->id );

                for( lnk = gilist; lnk->next != NULL; lnk = lnk->next ) {}
                lnk->next = gi;
            }
        }

/* -------------------------------------------------------------------- */
/*      Seek past the data.                                             */
/* -------------------------------------------------------------------- */


The following code will allow correct use of nested subgrids in the grid shift algrorithm.  While the example data we have only contains a single nested subgrid, this code should be usable for an arbitrary depth of nesting.  Previously, we would select the grid from the incorrect subgrid, and come up with an inaccurate shift for points contained in the subgrid.


/////////////////////////////////////////////////////////////////////////////////
Code snippet from pj_apply_gridshift.c, function pj_apply_gridshift_3, approx line 131.

    for( i = 0; i < point_count; i++ )
    {
        long io = i * point_offset;
        LP   input, output;
        int  itable;
        int foundSubgrid = 1;
            
        input.phi = y[io];
        input.lam = x[io];
        output.phi = HUGE_VAL;
        output.lam = HUGE_VAL;

        /* 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 != NULL && foundSubgrid == 1)
            {
                PJ_GRIDINFO *child;
                foundSubgrid = 0;
                for( child = gi->child; child != NULL; child = child->next )
                {
                    struct CTABLE *ct1 = child->ct;
                    double 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;

                    // Use this child instead, it is better
                    gi = child;
                    ct = child->ct;
                    foundSubgrid=1;  // Force us to search for a child of this subgrid
                    break;
                }
            }

            /* load the grid shift info if we don't have it. */
            if( ct->cvs == NULL && !pj_gridinfo_load( ctx, gi ) )
            {
                pj_ctx_set_errno( ctx, -38 );
                return -38;
            }
            
            output = nad_cvt( input, inverse, ct );
            if( output.lam != HUGE_VAL )
            {
                if( debug_count++ < 20 )
                    pj_log( ctx, PJ_LOG_DEBUG_MINOR,
                            "pj_apply_gridshift(): used %s", ct->id );
                break;
            }
        }

        if( output.lam == HUGE_VAL )

[bpayne]

Code snippets that fix bug

[bpayne]

Whoops! The datum file is too large to be attached, 1.7MB, even 1MB zipped. If you need a copy of the datum file for testing, feel free to contact me, I can email or ftp the file to you directly.