diff options
| -rw-r--r-- | src/4D_api.cpp | 62 | ||||
| -rw-r--r-- | src/Makefile.am | 3 | ||||
| -rw-r--r-- | src/apply_gridshift.cpp | 502 | ||||
| -rw-r--r-- | src/gridinfo.cpp | 834 | ||||
| -rw-r--r-- | src/gridlist.cpp | 220 | ||||
| -rw-r--r-- | src/grids.cpp | 578 | ||||
| -rw-r--r-- | src/grids.hpp | 51 | ||||
| -rw-r--r-- | src/init.cpp | 3 | ||||
| -rw-r--r-- | src/lib_proj.cmake | 5 | ||||
| -rw-r--r-- | src/malloc.cpp | 3 | ||||
| -rw-r--r-- | src/nad_cvt.cpp | 77 | ||||
| -rw-r--r-- | src/nad_init.cpp | 308 | ||||
| -rw-r--r-- | src/nad_intr.cpp | 67 | ||||
| -rw-r--r-- | src/proj_internal.h | 58 | ||||
| -rw-r--r-- | src/transform.cpp | 57 | ||||
| -rw-r--r-- | src/transformations/hgridshift.cpp | 4 |
16 files changed, 982 insertions, 1850 deletions
diff --git a/src/4D_api.cpp b/src/4D_api.cpp index ed6a62d6..c994df63 100644 --- a/src/4D_api.cpp +++ b/src/4D_api.cpp @@ -1553,12 +1553,9 @@ PJ_GRID_INFO proj_grid_info(const char *gridname) { /*PJ_CONTEXT *ctx = proj_context_create(); */ PJ_CONTEXT *ctx = pj_get_default_ctx(); - PJ_GRIDINFO *gridinfo = pj_gridinfo_init(ctx, gridname); memset(&grinfo, 0, sizeof(PJ_GRID_INFO)); - /* in case the grid wasn't found */ - if (gridinfo->filename == nullptr || gridinfo->ct == nullptr) { - + { const auto gridSet = NS_PROJ::VerticalShiftGridSet::open(ctx, gridname); if( gridSet ) { @@ -1591,44 +1588,49 @@ PJ_GRID_INFO proj_grid_info(const char *gridname) { grinfo.upperright.lam = extent.eastLon; grinfo.upperright.phi = extent.northLat; - pj_gridinfo_free(ctx, gridinfo); - return grinfo; } } - - pj_gridinfo_free(ctx, gridinfo); - strcpy(grinfo.format, "missing"); - return grinfo; } - /* The string copies below are automatically null-terminated due to */ - /* the memset above, so strncpy is safe */ - - /* name of grid */ - strncpy (grinfo.gridname, gridname, sizeof(grinfo.gridname) - 1); + { + const auto gridSet = NS_PROJ::HorizontalShiftGridSet::open(ctx, gridname); + if( gridSet ) + { + const auto& grids = gridSet->grids(); + if( !grids.empty() ) + { + const auto& grid = grids.front(); + const auto& extent = grid->extentAndRes(); - /* full path of grid */ - pj_find_file(ctx, gridname, grinfo.filename, sizeof(grinfo.filename) - 1); + /* name of grid */ + strncpy (grinfo.gridname, gridname, sizeof(grinfo.gridname) - 1); - /* grid format */ - strncpy (grinfo.format, gridinfo->format, sizeof(grinfo.format) - 1); + /* full path of grid */ + pj_find_file(ctx, gridname, grinfo.filename, sizeof(grinfo.filename) - 1); - /* grid size */ - grinfo.n_lon = gridinfo->ct->lim.lam; - grinfo.n_lat = gridinfo->ct->lim.phi; + /* grid format */ + strncpy (grinfo.format, gridSet->format().c_str(), sizeof(grinfo.format) - 1); - /* cell size */ - grinfo.cs_lon = gridinfo->ct->del.lam; - grinfo.cs_lat = gridinfo->ct->del.phi; + /* grid size */ + grinfo.n_lon = grid->width(); + grinfo.n_lat = grid->height(); - /* bounds of grid */ - grinfo.lowerleft = gridinfo->ct->ll; - grinfo.upperright.lam = grinfo.lowerleft.lam + (grinfo.n_lon-1)*grinfo.cs_lon; - grinfo.upperright.phi = grinfo.lowerleft.phi + (grinfo.n_lat-1)*grinfo.cs_lat; + /* cell size */ + grinfo.cs_lon = extent.resLon; + grinfo.cs_lat = extent.resLat; - pj_gridinfo_free(ctx, gridinfo); + /* bounds of grid */ + grinfo.lowerleft.lam = extent.westLon; + grinfo.lowerleft.phi = extent.southLat; + grinfo.upperright.lam = extent.eastLon; + grinfo.upperright.phi = extent.northLat; + return grinfo; + } + } + } + strcpy(grinfo.format, "missing"); return grinfo; } diff --git a/src/Makefile.am b/src/Makefile.am index 3d0868a7..a12de4e1 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -195,9 +195,8 @@ libproj_la_SOURCES = \ release.cpp gauss.cpp \ fileapi.cpp \ \ - nad_cvt.cpp nad_init.cpp nad_intr.cpp \ apply_gridshift.cpp datums.cpp datum_set.cpp transform.cpp \ - geocent.cpp geocent.h utils.cpp gridinfo.cpp gridlist.cpp \ + geocent.cpp geocent.h utils.cpp \ jniproj.cpp mutex.cpp initcache.cpp apply_vgridshift.cpp geodesic.c \ strtod.cpp \ \ diff --git a/src/apply_gridshift.cpp b/src/apply_gridshift.cpp index 2f46ca11..5a5404da 100644 --- a/src/apply_gridshift.cpp +++ b/src/apply_gridshift.cpp @@ -28,7 +28,9 @@ * DEALINGS IN THE SOFTWARE. *****************************************************************************/ -#define PJ_LIB__ +#ifndef FROM_PROJ_CPP +#define FROM_PROJ_CPP +#endif #include <math.h> #include <stdio.h> @@ -36,231 +38,62 @@ #include "proj.h" #include "proj_internal.h" +#include "proj/internal/internal.hpp" -/************************************************************************/ -/* 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. */ -/************************************************************************/ +using namespace NS_PROJ; -int pj_apply_gridshift( projCtx ctx, const char *nadgrids, int inverse, - long point_count, int point_offset, - double *x, double *y, double *z ) +// --------------------------------------------------------------------------- +static const HorizontalShiftGrid* findGrid( + const std::vector<std::unique_ptr<HorizontalShiftGridSet>>& grids, + PJ_LP input) { - 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->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. */ -/************************************************************************/ - -static 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++ ) + for( const auto& gridset: grids ) { - - 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; + auto grid = gridset->gridAt(input.lam, input.phi); + if( grid ) + return grid; } - 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 ) +static std::vector<std::unique_ptr<NS_PROJ::HorizontalShiftGridSet>> + getListOfGridSets(PJ_CONTEXT* ctx, + const char* grids) { - int i; - struct CTABLE *ct; - static int debug_count = 0; - (void) z; - - if( gridlist== nullptr || gridlist_count == 0 ) + std::vector<std::unique_ptr<NS_PROJ::HorizontalShiftGridSet>> list; + auto listOfGrids = internal::split(std::string(grids), ','); + for( const auto& grid: listOfGrids ) { - 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 ) + const char* gridname = grid.c_str(); + bool canFail = false; + if( gridname[0] == '@' ) { - output = nad_cvt( input, inverse, ct ); - - if ( output.lam != HUGE_VAL && debug_count++ < 20 ) - pj_log( ctx, PJ_LOG_DEBUG_MINOR, "pj_apply_gridshift(): used %s", ct->id ); + canFail = true; + gridname ++; } - - if ( output.lam == HUGE_VAL ) + auto gridSet = HorizontalShiftGridSet::open(ctx, gridname); + if( !gridSet ) { - if( ctx->debug_level >= PJ_LOG_DEBUG_MAJOR ) + if( !canFail ) { - 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 ); - } + pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); + return {}; } - - /* - * 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; + list.emplace_back(std::move(gridSet)); } } - - return 0; + return list; } + /**********************************************/ -int proj_hgrid_init(PJ* P, const char *grids) { +int proj_hgrid_init(PJ* P, const char *gridkey) { /********************************************** Initizalize and populate list of horizontal @@ -275,29 +108,158 @@ int proj_hgrid_init(PJ* P, const char *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); + std::string key("s"); + key += gridkey; + const char* grids = pj_param(P->ctx, P->params, key.c_str()).s; + if( grids == nullptr ) + return 0; - if (P->gridlist == nullptr) { - P->gridlist = pj_gridlist_from_nadgrids( - P->ctx, - pj_param(P->ctx, P->params, sgrids).s, - &(P->gridlist_count) - ); + P->hgrids = getListOfGridSets(P->ctx, grids); + return static_cast<int>(P->hgrids.size()); +} - if( P->gridlist == nullptr || P->gridlist_count == 0 ) { - pj_dealloc(sgrids); - return 0; +typedef struct { pj_int32 lam, phi; } ILP; + +static PJ_LP nad_intr(PJ_LP t, const HorizontalShiftGrid* grid) { + PJ_LP val, frct; + ILP indx; + int in; + + const auto& extent = grid->extentAndRes(); + t.lam /= extent.resLon; + indx.lam = isnan(t.lam) ? 0 : (pj_int32)lround(floor(t.lam)); + t.phi /= extent.resLat; + indx.phi = isnan(t.phi) ? 0 : (pj_int32)lround(floor(t.phi)); + + frct.lam = t.lam - indx.lam; + frct.phi = t.phi - indx.phi; + val.lam = val.phi = HUGE_VAL; + if (indx.lam < 0) { + if (indx.lam == -1 && frct.lam > 0.99999999999) { + ++indx.lam; + frct.lam = 0.; + } else + return val; + } else if ((in = indx.lam + 1) >= grid->width()) { + if (in == grid->width() && frct.lam < 1e-11) { + --indx.lam; + frct.lam = 1.; + } else + return val; + } + if (indx.phi < 0) { + if (indx.phi == -1 && frct.phi > 0.99999999999) { + ++indx.phi; + frct.phi = 0.; + } else + return val; + } else if ((in = indx.phi + 1) >= grid->height()) { + if (in == grid->height() && frct.phi < 1e-11) { + --indx.phi; + frct.phi = 1.; + } else + return val; + } + + float f00Lon = 0, f00Lat = 0; + float f10Lon = 0, f10Lat = 0; + float f01Lon = 0, f01Lat = 0; + float f11Lon = 0, f11Lat = 0; + if( !grid->valueAt(indx.lam, indx.phi, f00Lon, f00Lat) || + !grid->valueAt(indx.lam + 1, indx.phi, f10Lon, f10Lat) || + !grid->valueAt(indx.lam, indx.phi + 1, f01Lon, f01Lat) || + !grid->valueAt(indx.lam + 1, indx.phi + 1, f11Lon, f11Lat) ) + { + return val; } + + double m10 = frct.lam; + double m11 = m10; + double m01 = 1. - frct.lam; + double m00 = m01; + m11 *= frct.phi; + m01 *= frct.phi; + frct.phi = 1. - frct.phi; + m00 *= frct.phi; + m10 *= frct.phi; + val.lam = m00 * f00Lon + m10 * f10Lon + + m01 * f01Lon + m11 * f11Lon; + val.phi = m00 * f00Lat + m10 * f10Lat + + m01 * f01Lat + m11 * f11Lat; + return val; +} + + +#define MAX_ITERATIONS 10 +#define TOL 1e-12 + +static +PJ_LP nad_cvt(PJ_LP in, int inverse, const HorizontalShiftGrid* grid) { + PJ_LP t, tb,del, dif; + int i = MAX_ITERATIONS; + const double toltol = TOL*TOL; + + if (in.lam == HUGE_VAL) + return in; + + /* normalize input to ll origin */ + tb = in; + const auto& extent = grid->extentAndRes(); + tb.lam -= extent.westLon; + tb.phi -= extent.southLat; + + tb.lam = adjlon (tb.lam - M_PI) + M_PI; + + t = nad_intr (tb, grid); + if (t.lam == HUGE_VAL) + return t; + + if (!inverse) { + in.lam -= t.lam; + in.phi += t.phi; + return in; } - if (P->gridlist_count == 0) { - proj_errno_set(P, PJD_ERR_FAILED_TO_LOAD_GRID); + t.lam = tb.lam + t.lam; + t.phi = tb.phi - t.phi; + + do { + del = nad_intr(t, grid); + + /* This case used to return failure, but I have + changed it to return the first order approximation + of the inverse shift. This avoids cases where the + grid shift *into* this grid came from another grid. + While we aren't returning optimally correct results + I feel a close result in this case is better than + no result. NFW + To demonstrate use -112.5839956 49.4914451 against + the NTv2 grid shift file from Canada. */ + if (del.lam == HUGE_VAL) + break; + + dif.lam = t.lam - del.lam - tb.lam; + dif.phi = t.phi + del.phi - tb.phi; + t.lam -= dif.lam; + t.phi -= dif.phi; + + } while (--i && (dif.lam*dif.lam + dif.phi*dif.phi > toltol)); /* prob. slightly faster than hypot() */ + + if (i==0) { + /* If we had access to a context, this should go through pj_log, and we should set ctx->errno */ + if (getenv ("PROJ_DEBUG")) + fprintf( stderr, "Inverse grid shift iterator failed to converge.\n" ); + t.lam = t.phi = HUGE_VAL; + return t; } - pj_dealloc(sgrids); - return P->gridlist_count; + /* and again: pj_log and ctx->errno */ + if (del.lam==HUGE_VAL && getenv ("PROJ_DEBUG")) + fprintf (stderr, "Inverse grid shift iteration failed, presumably at grid edge.\nUsing first approximation.\n"); + + in.lam = adjlon (t.lam + extent.westLon); + in.phi = t.phi + extent.southLat; + return in; } /********************************************/ @@ -306,22 +268,22 @@ int proj_hgrid_init(PJ* P, const char *grids) { /* 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); + const auto grid = findGrid(P->hgrids, lp); + if( !grid ) { + 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; + const auto& extent = grid->extentAndRes(); + lp.lam -= extent.westLon; + lp.phi -= extent.southLat; lp.lam = adjlon(lp.lam - M_PI) + M_PI; - out = nad_intr(lp, ct); + out = nad_intr(lp, grid); if (out.lam == HUGE_VAL || out.phi == HUGE_VAL) { pj_ctx_set_errno(P->ctx, PJD_ERR_GRID_AREA); @@ -330,33 +292,89 @@ PJ_LP proj_hgrid_value(PJ *P, PJ_LP lp) { return out; } -PJ_LP proj_hgrid_apply(PJ *P, PJ_LP lp, PJ_DIRECTION direction) { - struct CTABLE *ct; - int inverse; +static +PJ_LP proj_hgrid_apply_internal(PJ_CONTEXT* ctx, + PJ_LP lp, + PJ_DIRECTION direction, + const std::vector<std::unique_ptr<HorizontalShiftGridSet>>& grids) +{ PJ_LP out; - out.lam = HUGE_VAL; out.phi = HUGE_VAL; + 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 ); - } + const auto grid = findGrid(grids, lp); + if( !grid ) { + pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); return out; } + if( grid->isNullGrid() ) + { + return lp; + } - inverse = direction == PJ_FWD ? 0 : 1; - out = nad_cvt(lp, inverse, ct); + int inverse = direction == PJ_FWD ? 0 : 1; + out = nad_cvt(lp, inverse, grid); if (out.lam == HUGE_VAL || out.phi == HUGE_VAL) - pj_ctx_set_errno(P->ctx, PJD_ERR_GRID_AREA); + pj_ctx_set_errno(ctx, PJD_ERR_GRID_AREA); return out; +} +PJ_LP proj_hgrid_apply(PJ *P, PJ_LP lp, PJ_DIRECTION direction) { + return proj_hgrid_apply_internal(P->ctx, lp, direction, P->hgrids); +} + +/************************************************************************/ +/* 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 */ ) + +{ + auto hgrids = getListOfGridSets(ctx, nadgrids); + if( hgrids.empty() ) + { + pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); + return 1; + } + + + for( long i = 0; i < point_count; i++ ) + { + PJ_LP input; + + long io = i * point_offset; + input.phi = y[io]; + input.lam = x[io]; + + auto output = proj_hgrid_apply_internal(ctx, input, inverse ? PJ_INV : PJ_FWD, hgrids); + + if ( output.lam != HUGE_VAL ) + { + y[io] = output.phi; + x[io] = output.lam; + } + else + { + 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 ); + } + } + } + + return 0; } diff --git a/src/gridinfo.cpp b/src/gridinfo.cpp deleted file mode 100644 index 01d392b3..00000000 --- a/src/gridinfo.cpp +++ /dev/null @@ -1,834 +0,0 @@ -/****************************************************************************** - * Project: PROJ.4 - * Purpose: Functions for handling individual PJ_GRIDINFO's. Includes - * loaders for all formats but CTABLE (in nad_init.c). - * 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 <errno.h> -#include <math.h> -#include <stddef.h> -#include <stdio.h> -#include <stdlib.h> -#include <string.h> - -#include "proj_internal.h" - -/************************************************************************/ -/* swap_words() */ -/* */ -/* Convert the byte order of the given word(s) in place. */ -/************************************************************************/ - -static const int byte_order_test = 1; -#define IS_LSB (1 == ((const unsigned char *) (&byte_order_test))[0]) - -static void swap_words( unsigned char *data, int word_size, int word_count ) - -{ - int word; - - for( word = 0; word < word_count; word++ ) - { - int i; - - for( i = 0; i < word_size/2; i++ ) - { - unsigned char t; - - t = data[i]; - data[i] = data[word_size-i-1]; - data[word_size-i-1] = t; - } - - data += word_size; - } -} - -/************************************************************************/ -/* to_double() */ -/* */ -/* Returns a double from the pointed data. */ -/************************************************************************/ - -static double to_double( unsigned char* data ) -{ - double d; - memcpy(&d, data, sizeof(d)); - return d; -} - -/************************************************************************/ -/* pj_gridinfo_free() */ -/************************************************************************/ - -void pj_gridinfo_free( projCtx ctx, PJ_GRIDINFO *gi ) - -{ - if( gi == nullptr ) - return; - - if( gi->child != nullptr ) - { - PJ_GRIDINFO *child, *next; - - for( child = gi->child; child != nullptr; child=next) - { - next=child->next; - pj_gridinfo_free( ctx, child ); - } - } - - if( gi->ct != nullptr ) - nad_free( gi->ct ); - - free( gi->gridname ); - if( gi->filename != nullptr ) - free( gi->filename ); - - pj_dalloc( gi ); -} - -/************************************************************************/ -/* pj_gridinfo_load() */ -/* */ -/* This function is intended to implement delayed loading of */ -/* the data contents of a grid file. The header and related */ -/* stuff are loaded by pj_gridinfo_init(). */ -/************************************************************************/ - -int pj_gridinfo_load( projCtx_t* ctx, PJ_GRIDINFO *gi ) - -{ - struct CTABLE ct_tmp; - - if( gi == nullptr || gi->ct == nullptr ) - return 0; - - pj_acquire_lock(); - if( gi->ct->cvs != nullptr ) - { - pj_release_lock(); - return 1; - } - - memcpy(&ct_tmp, gi->ct, sizeof(struct CTABLE)); - -/* -------------------------------------------------------------------- */ -/* Original platform specific CTable format. */ -/* -------------------------------------------------------------------- */ - if( strcmp(gi->format,"ctable") == 0 ) - { - PAFile fid; - int result; - - fid = pj_open_lib( ctx, gi->filename, "rb" ); - - if( fid == nullptr ) - { - pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); - pj_release_lock(); - return 0; - } - - result = nad_ctable_load( ctx, &ct_tmp, (struct projFileAPI_t*)fid ); - - pj_ctx_fclose( ctx, fid ); - - gi->ct->cvs = ct_tmp.cvs; - pj_release_lock(); - - return result; - } - -/* -------------------------------------------------------------------- */ -/* CTable2 format. */ -/* -------------------------------------------------------------------- */ - else if( strcmp(gi->format,"ctable2") == 0 ) - { - PAFile fid; - int result; - - fid = pj_open_lib( ctx, gi->filename, "rb" ); - - if( fid == nullptr ) - { - pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); - pj_release_lock(); - return 0; - } - - result = nad_ctable2_load( ctx, &ct_tmp, (struct projFileAPI_t*)fid ); - - pj_ctx_fclose( ctx, fid ); - - gi->ct->cvs = ct_tmp.cvs; - - pj_release_lock(); - return result; - } - -/* -------------------------------------------------------------------- */ -/* NTv1 format. */ -/* We process one line at a time. Note that the array storage */ -/* direction (e-w) is different in the NTv1 file and what */ -/* the CTABLE is supposed to have. The phi/lam are also */ -/* reversed, and we have to be aware of byte swapping. */ -/* -------------------------------------------------------------------- */ - else if( strcmp(gi->format,"ntv1") == 0 ) - { - double *row_buf; - int row; - PAFile fid; - - fid = pj_open_lib( ctx, gi->filename, "rb" ); - - if( fid == nullptr ) - { - pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); - pj_release_lock(); - return 0; - } - - pj_ctx_fseek( ctx, fid, gi->grid_offset, SEEK_SET ); - - row_buf = (double *) pj_malloc(gi->ct->lim.lam * sizeof(double) * 2); - ct_tmp.cvs = (FLP *) pj_malloc(gi->ct->lim.lam*gi->ct->lim.phi*sizeof(FLP)); - if( row_buf == nullptr || ct_tmp.cvs == nullptr ) - { - pj_dalloc( row_buf ); - pj_dalloc( ct_tmp.cvs ); - pj_ctx_set_errno( ctx, ENOMEM ); - pj_release_lock(); - return 0; - } - - for( row = 0; row < gi->ct->lim.phi; row++ ) - { - int i; - FLP *cvs; - double *diff_seconds; - - if( pj_ctx_fread( ctx, row_buf, - sizeof(double), gi->ct->lim.lam * 2, fid ) - != (size_t)( 2 * gi->ct->lim.lam ) ) - { - pj_dalloc( row_buf ); - pj_dalloc( ct_tmp.cvs ); - pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); - pj_release_lock(); - return 0; - } - - if( IS_LSB ) - swap_words( (unsigned char *) row_buf, 8, gi->ct->lim.lam*2 ); - - /* convert seconds to radians */ - diff_seconds = row_buf; - - for( i = 0; i < gi->ct->lim.lam; i++ ) - { - cvs = ct_tmp.cvs + (row) * gi->ct->lim.lam - + (gi->ct->lim.lam - i - 1); - - cvs->phi = (float) (*(diff_seconds++) * ((M_PI/180.0) / 3600.0)); - cvs->lam = (float) (*(diff_seconds++) * ((M_PI/180.0) / 3600.0)); - } - } - - pj_dalloc( row_buf ); - - pj_ctx_fclose( ctx, fid ); - - gi->ct->cvs = ct_tmp.cvs; - pj_release_lock(); - - return 1; - } - -/* -------------------------------------------------------------------- */ -/* NTv2 format. */ -/* We process one line at a time. Note that the array storage */ -/* direction (e-w) is different in the NTv2 file and what */ -/* the CTABLE is supposed to have. The phi/lam are also */ -/* reversed, and we have to be aware of byte swapping. */ -/* -------------------------------------------------------------------- */ - else if( strcmp(gi->format,"ntv2") == 0 ) - { - float *row_buf; - int row; - PAFile fid; - - pj_log( ctx, PJ_LOG_DEBUG_MINOR, - "NTv2 - loading grid %s", gi->ct->id ); - - fid = pj_open_lib( ctx, gi->filename, "rb" ); - - if( fid == nullptr ) - { - pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); - pj_release_lock(); - return 0; - } - - pj_ctx_fseek( ctx, fid, gi->grid_offset, SEEK_SET ); - - row_buf = (float *) pj_malloc(gi->ct->lim.lam * sizeof(float) * 4); - ct_tmp.cvs = (FLP *) pj_malloc(gi->ct->lim.lam*gi->ct->lim.phi*sizeof(FLP)); - if( row_buf == nullptr || ct_tmp.cvs == nullptr ) - { - pj_dalloc( row_buf ); - pj_dalloc( ct_tmp.cvs ); - pj_ctx_set_errno( ctx, ENOMEM ); - pj_release_lock(); - return 0; - } - - for( row = 0; row < gi->ct->lim.phi; row++ ) - { - int i; - FLP *cvs; - float *diff_seconds; - - if( pj_ctx_fread( ctx, row_buf, sizeof(float), - gi->ct->lim.lam*4, fid ) - != (size_t)( 4 * gi->ct->lim.lam ) ) - { - pj_dalloc( row_buf ); - pj_dalloc( ct_tmp.cvs ); - pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); - pj_release_lock(); - return 0; - } - - if( gi->must_swap ) - swap_words( (unsigned char *) row_buf, 4, - gi->ct->lim.lam*4 ); - - /* convert seconds to radians */ - diff_seconds = row_buf; - - for( i = 0; i < gi->ct->lim.lam; i++ ) - { - cvs = ct_tmp.cvs + (row) * gi->ct->lim.lam - + (gi->ct->lim.lam - i - 1); - - cvs->phi = (float) (*(diff_seconds++) * ((M_PI/180.0) / 3600.0)); - cvs->lam = (float) (*(diff_seconds++) * ((M_PI/180.0) / 3600.0)); - diff_seconds += 2; /* skip accuracy values */ - } - } - - pj_dalloc( row_buf ); - - pj_ctx_fclose( ctx, fid ); - - gi->ct->cvs = ct_tmp.cvs; - - pj_release_lock(); - return 1; - } - - else - { - pj_release_lock(); - return 0; - } -} - -/************************************************************************/ -/* gridinfo_parent() */ -/* */ -/* Seek a parent grid file by name from a grid list */ -/************************************************************************/ - -static PJ_GRIDINFO* gridinfo_parent( PJ_GRIDINFO *gilist, - const char *name, int length ) -{ - while( gilist ) - { - if( strncmp(gilist->ct->id,name,length) == 0 ) return gilist; - if( gilist->child ) - { - PJ_GRIDINFO *parent=gridinfo_parent( gilist->child, name, length ); - if( parent ) return parent; - } - gilist=gilist->next; - } - return gilist; -} - -/************************************************************************/ -/* pj_gridinfo_init_ntv2() */ -/* */ -/* Load a ntv2 (.gsb) file. */ -/************************************************************************/ - -static int pj_gridinfo_init_ntv2( projCtx ctx, PAFile fid, PJ_GRIDINFO *gilist ) -{ - unsigned char header[11*16]; - int num_subfiles, subfile; - int must_swap; - - /* cppcheck-suppress sizeofCalculation */ - STATIC_ASSERT( sizeof(pj_int32) == 4 ); - /* cppcheck-suppress sizeofCalculation */ - STATIC_ASSERT( sizeof(double) == 8 ); - -/* -------------------------------------------------------------------- */ -/* Read the overview header. */ -/* -------------------------------------------------------------------- */ - if( pj_ctx_fread( ctx, header, sizeof(header), 1, fid ) != 1 ) - { - pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); - return 0; - } - - if( header[8] == 11 ) - must_swap = !IS_LSB; - else - must_swap = IS_LSB; - -/* -------------------------------------------------------------------- */ -/* Byte swap interesting fields if needed. */ -/* -------------------------------------------------------------------- */ - if( must_swap ) - { - swap_words( header+8, 4, 1 ); - swap_words( header+8+16, 4, 1 ); - swap_words( header+8+32, 4, 1 ); - swap_words( header+8+7*16, 8, 1 ); - swap_words( header+8+8*16, 8, 1 ); - swap_words( header+8+9*16, 8, 1 ); - swap_words( header+8+10*16, 8, 1 ); - } - -/* -------------------------------------------------------------------- */ -/* Get the subfile count out ... all we really use for now. */ -/* -------------------------------------------------------------------- */ - memcpy( &num_subfiles, header+8+32, 4 ); - -/* ==================================================================== */ -/* Step through the subfiles, creating a PJ_GRIDINFO for each. */ -/* ==================================================================== */ - for( subfile = 0; subfile < num_subfiles; subfile++ ) - { - struct CTABLE *ct; - PJ_LP ur; - int gs_count; - PJ_GRIDINFO *gi; - -/* -------------------------------------------------------------------- */ -/* Read header. */ -/* -------------------------------------------------------------------- */ - if( pj_ctx_fread( ctx, header, sizeof(header), 1, fid ) != 1 ) - { - pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); - return 0; - } - - if( strncmp((const char *) header,"SUB_NAME",8) != 0 ) - { - pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); - return 0; - } - -/* -------------------------------------------------------------------- */ -/* Byte swap interesting fields if needed. */ -/* -------------------------------------------------------------------- */ - if( must_swap ) - { - swap_words( header+8+16*4, 8, 1 ); - swap_words( header+8+16*5, 8, 1 ); - swap_words( header+8+16*6, 8, 1 ); - swap_words( header+8+16*7, 8, 1 ); - swap_words( header+8+16*8, 8, 1 ); - swap_words( header+8+16*9, 8, 1 ); - swap_words( header+8+16*10, 4, 1 ); - } - -/* -------------------------------------------------------------------- */ -/* Initialize a corresponding "ct" structure. */ -/* -------------------------------------------------------------------- */ - ct = (struct CTABLE *) pj_malloc(sizeof(struct CTABLE)); - if (!ct) { - pj_ctx_set_errno(ctx, ENOMEM); - return 0; - } - strncpy( ct->id, (const char *) header + 8, 8 ); - ct->id[8] = '\0'; - - ct->ll.lam = - to_double(header+7*16+8); /* W_LONG */ - ct->ll.phi = to_double(header+4*16+8); /* S_LAT */ - - ur.lam = - to_double(header+6*16+8); /* E_LONG */ - ur.phi = to_double(header+5*16+8); /* N_LAT */ - - ct->del.lam = to_double(header+9*16+8); - ct->del.phi = to_double(header+8*16+8); - - ct->lim.lam = (pj_int32) (fabs(ur.lam-ct->ll.lam)/ct->del.lam + 0.5) + 1; - ct->lim.phi = (pj_int32) (fabs(ur.phi-ct->ll.phi)/ct->del.phi + 0.5) + 1; - - pj_log( ctx, PJ_LOG_DEBUG_MINOR, - "NTv2 %s %dx%d: LL=(%.9g,%.9g) UR=(%.9g,%.9g)", - ct->id, - ct->lim.lam, ct->lim.phi, - ct->ll.lam/3600.0, ct->ll.phi/3600.0, - ur.lam/3600.0, ur.phi/3600.0 ); - - ct->ll.lam *= DEG_TO_RAD/3600.0; - ct->ll.phi *= DEG_TO_RAD/3600.0; - ct->del.lam *= DEG_TO_RAD/3600.0; - ct->del.phi *= DEG_TO_RAD/3600.0; - - memcpy( &gs_count, header + 8 + 16*10, 4 ); - if( gs_count != ct->lim.lam * ct->lim.phi ) - { - pj_log( ctx, PJ_LOG_ERROR, - "GS_COUNT(%d) does not match expected cells (%dx%d=%d)", - gs_count, ct->lim.lam, ct->lim.phi, - ct->lim.lam * ct->lim.phi ); - pj_dalloc(ct); - pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); - return 0; - } - - ct->cvs = nullptr; - -/* -------------------------------------------------------------------- */ -/* Create a new gridinfo for this if we aren't processing the */ -/* 1st subfile, and initialize our grid info. */ -/* -------------------------------------------------------------------- */ - if( subfile == 0 ) - gi = gilist; - else - { - gi = (PJ_GRIDINFO *) pj_calloc(1, sizeof(PJ_GRIDINFO)); - if (!gi) { - pj_dalloc(ct); - pj_gridinfo_free(ctx, gilist); - pj_ctx_set_errno(ctx, ENOMEM); - return 0; - } - - gi->gridname = pj_strdup( gilist->gridname ); - gi->filename = pj_strdup( gilist->filename ); - if (!gi->gridname || !gi->filename) { - pj_gridinfo_free(ctx, gi); - pj_dalloc(ct); - pj_gridinfo_free(ctx, gilist); - pj_ctx_set_errno(ctx, ENOMEM); - return 0; - } - gi->next = nullptr; - } - - gi->must_swap = must_swap; - gi->ct = ct; - gi->format = "ntv2"; - gi->grid_offset = pj_ctx_ftell( ctx, fid ); - -/* -------------------------------------------------------------------- */ -/* 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 != nullptr; lnk = lnk->next ) {} - lnk->next = gi; - } - } - - else - { - PJ_GRIDINFO *lnk; - PJ_GRIDINFO *gp = gridinfo_parent(gilist, - (const char*)header+24,8); - - if( gp == nullptr ) - { - pj_log( ctx, PJ_LOG_ERROR, - "pj_gridinfo_init_ntv2(): " - "failed to find parent %8.8s for %s.", - (const char *) header+24, gi->ct->id ); - - for( lnk = gilist; lnk->next != nullptr; lnk = lnk->next ) {} - lnk->next = gi; - } - else - { - if( gp->child == nullptr ) - { - gp->child = gi; - } - else - { - for( lnk = gp->child; lnk->next != nullptr; lnk = lnk->next ) {} - lnk->next = gi; - } - } - } - -/* -------------------------------------------------------------------- */ -/* Seek past the data. */ -/* -------------------------------------------------------------------- */ - pj_ctx_fseek( ctx, fid, gs_count * 16, SEEK_CUR ); - } - - return 1; -} - -/************************************************************************/ -/* pj_gridinfo_init_ntv1() */ -/* */ -/* Load an NTv1 style Canadian grid shift file. */ -/************************************************************************/ - -static int pj_gridinfo_init_ntv1( projCtx ctx, PAFile fid, PJ_GRIDINFO *gi ) - -{ - unsigned char header[192]; /* 12 records of 16 bytes */ - struct CTABLE *ct; - PJ_LP ur; - - /* cppcheck-suppress sizeofCalculation */ - STATIC_ASSERT( sizeof(pj_int32) == 4 ); - /* cppcheck-suppress sizeofCalculation */ - STATIC_ASSERT( sizeof(double) == 8 ); - -/* -------------------------------------------------------------------- */ -/* Read the header. */ -/* -------------------------------------------------------------------- */ - if( pj_ctx_fread( ctx, header, sizeof(header), 1, fid ) != 1 ) - { - pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); - return 0; - } - -/* -------------------------------------------------------------------- */ -/* Regularize fields of interest. */ -/* -------------------------------------------------------------------- */ - if( IS_LSB ) - { - swap_words( header+8, 4, 1 ); - swap_words( header+24, 8, 1 ); - swap_words( header+40, 8, 1 ); - swap_words( header+56, 8, 1 ); - swap_words( header+72, 8, 1 ); - swap_words( header+88, 8, 1 ); - swap_words( header+104, 8, 1 ); - } - - if( *((int *) (header+8)) != 12 ) - { - pj_log( ctx, PJ_LOG_ERROR, - "NTv1 grid shift file has wrong record count, corrupt?" ); - pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); - return 0; - } - -/* -------------------------------------------------------------------- */ -/* Fill in CTABLE structure. */ -/* -------------------------------------------------------------------- */ - ct = (struct CTABLE *) pj_malloc(sizeof(struct CTABLE)); - if (!ct) { - pj_ctx_set_errno(ctx, ENOMEM); - return 0; - } - strcpy( ct->id, "NTv1 Grid Shift File" ); - - ct->ll.lam = - to_double(header+72); - ct->ll.phi = to_double(header+24); - ur.lam = - to_double(header+56); - ur.phi = to_double(header+40); - ct->del.lam = to_double(header+104); - ct->del.phi = to_double(header+88); - ct->lim.lam = (pj_int32) (fabs(ur.lam-ct->ll.lam)/ct->del.lam + 0.5) + 1; - ct->lim.phi = (pj_int32) (fabs(ur.phi-ct->ll.phi)/ct->del.phi + 0.5) + 1; - - pj_log( ctx, PJ_LOG_DEBUG_MINOR, - "NTv1 %dx%d: LL=(%.9g,%.9g) UR=(%.9g,%.9g)", - ct->lim.lam, ct->lim.phi, - ct->ll.lam, ct->ll.phi, ur.lam, ur.phi ); - - ct->ll.lam *= DEG_TO_RAD; - ct->ll.phi *= DEG_TO_RAD; - ct->del.lam *= DEG_TO_RAD; - ct->del.phi *= DEG_TO_RAD; - ct->cvs = nullptr; - - gi->ct = ct; - gi->grid_offset = (long) sizeof(header); - gi->format = "ntv1"; - - return 1; -} - -/************************************************************************/ -/* pj_gridinfo_init() */ -/* */ -/* Open and parse header details from a datum gridshift file */ -/* returning a list of PJ_GRIDINFOs for the grids in that */ -/* file. This superceeds use of nad_init() for modern */ -/* applications. */ -/************************************************************************/ - -PJ_GRIDINFO *pj_gridinfo_init( projCtx ctx, const char *gridname ) - -{ - PJ_GRIDINFO *gilist; - PAFile fp; - char header[160]; - size_t header_size = 0; - - errno = pj_errno = 0; - ctx->last_errno = 0; - -/* -------------------------------------------------------------------- */ -/* Initialize a GRIDINFO with stub info we would use if it */ -/* cannot be loaded. */ -/* -------------------------------------------------------------------- */ - gilist = (PJ_GRIDINFO *) pj_calloc(1, sizeof(PJ_GRIDINFO)); - if (!gilist) { - pj_ctx_set_errno(ctx, ENOMEM); - return nullptr; - } - - gilist->gridname = pj_strdup( gridname ); - if (!gilist->gridname) { - pj_dalloc(gilist); - pj_ctx_set_errno(ctx, ENOMEM); - return nullptr; - } - gilist->filename = nullptr; - gilist->format = "missing"; - gilist->grid_offset = 0; - gilist->ct = nullptr; - gilist->next = nullptr; - -/* -------------------------------------------------------------------- */ -/* Open the file using the usual search rules. */ -/* -------------------------------------------------------------------- */ - if (!(fp = pj_open_lib(ctx, gridname, "rb"))) { - ctx->last_errno = 0; /* don't treat as a persistent error */ - return gilist; - } - - gilist->filename = pj_strdup(gridname); - if (!gilist->filename) { - pj_dalloc(gilist->gridname); - pj_dalloc(gilist); - pj_ctx_set_errno(ctx, ENOMEM); - return nullptr; - } - -/* -------------------------------------------------------------------- */ -/* Load a header, to determine the file type. */ -/* -------------------------------------------------------------------- */ - if( (header_size = pj_ctx_fread( ctx, header, 1, - sizeof(header), fp ) ) != sizeof(header) ) - { - /* some files may be smaller that sizeof(header), eg 160, so */ - ctx->last_errno = 0; /* don't treat as a persistent error */ - pj_log( ctx, PJ_LOG_DEBUG_MAJOR, - "pj_gridinfo_init: short header read of %d bytes", - (int)header_size ); - } - - pj_ctx_fseek( ctx, fp, SEEK_SET, 0 ); - -/* -------------------------------------------------------------------- */ -/* Determine file type. */ -/* -------------------------------------------------------------------- */ - if( header_size >= 144 + 16 - && strncmp(header + 0, "HEADER", 6) == 0 - && strncmp(header + 96, "W GRID", 6) == 0 - && strncmp(header + 144, "TO NAD83 ", 16) == 0 ) - { - pj_gridinfo_init_ntv1( ctx, fp, gilist ); - } - - else if( header_size >= 48 + 7 - && strncmp(header + 0, "NUM_OREC", 8) == 0 - && strncmp(header + 48, "GS_TYPE", 7) == 0 ) - { - pj_gridinfo_init_ntv2( ctx, fp, gilist ); - } - - else if( header_size >= 9 && strncmp(header + 0,"CTABLE V2",9) == 0 ) - { - struct CTABLE *ct = nad_ctable2_init( ctx, (struct projFileAPI_t*)fp ); - - gilist->format = "ctable2"; - gilist->ct = ct; - - if (ct == nullptr) - { - pj_log( ctx, PJ_LOG_DEBUG_MAJOR, - "CTABLE V2 ct is NULL."); - } - else - { - pj_log( ctx, PJ_LOG_DEBUG_MAJOR, - "Ctable2 %s %dx%d: LL=(%.9g,%.9g) UR=(%.9g,%.9g)", - ct->id, - ct->lim.lam, ct->lim.phi, - ct->ll.lam * RAD_TO_DEG, ct->ll.phi * RAD_TO_DEG, - (ct->ll.lam + (ct->lim.lam-1)*ct->del.lam) * RAD_TO_DEG, - (ct->ll.phi + (ct->lim.phi-1)*ct->del.phi) * RAD_TO_DEG ); - } - } - - else - { - struct CTABLE *ct = nad_ctable_init( ctx, (struct projFileAPI_t*)fp ); - if (ct == nullptr) - { - pj_log( ctx, PJ_LOG_DEBUG_MAJOR, - "CTABLE ct is NULL."); - } else - { - gilist->format = "ctable"; - gilist->ct = ct; - - pj_log( ctx, PJ_LOG_DEBUG_MAJOR, - "Ctable %s %dx%d: LL=(%.9g,%.9g) UR=(%.9g,%.9g)", - ct->id, - ct->lim.lam, ct->lim.phi, - ct->ll.lam * RAD_TO_DEG, ct->ll.phi * RAD_TO_DEG, - (ct->ll.lam + (ct->lim.lam-1)*ct->del.lam) * RAD_TO_DEG, - (ct->ll.phi + (ct->lim.phi-1)*ct->del.phi) * RAD_TO_DEG ); - } - } - - pj_ctx_fclose(ctx, fp); - - return gilist; -} diff --git a/src/gridlist.cpp b/src/gridlist.cpp deleted file mode 100644 index c540b134..00000000 --- a/src/gridlist.cpp +++ /dev/null @@ -1,220 +0,0 @@ -/****************************************************************************** - * Project: PROJ.4 - * Purpose: Code to manage the list of currently loaded (cached) PJ_GRIDINFOs - * See pj_gridinfo.c for details of loading individual grids. - * 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 <errno.h> -#include <stddef.h> -#include <string.h> - -#include "proj.h" -#include "proj_internal.h" - -static PJ_GRIDINFO *grid_list = nullptr; -#define PJ_MAX_PATH_LENGTH 1024 - -/************************************************************************/ -/* pj_deallocate_grids() */ -/* */ -/* Deallocate all loaded grids. */ -/************************************************************************/ - -void pj_deallocate_grids() - -{ - while( grid_list != nullptr ) - { - PJ_GRIDINFO *item = grid_list; - grid_list = grid_list->next; - item->next = nullptr; - - pj_gridinfo_free( pj_get_default_ctx(), item ); - } -} - -/************************************************************************/ -/* pj_gridlist_merge_grid() */ -/* */ -/* Find/load the named gridfile and merge it into the */ -/* last_nadgrids_list. */ -/************************************************************************/ - -static int pj_gridlist_merge_gridfile( projCtx ctx, - const char *gridname, - PJ_GRIDINFO ***p_gridlist, - int *p_gridcount, - int *p_gridmax ) - -{ - int got_match=0; - PJ_GRIDINFO *this_grid, *tail = nullptr; - -/* -------------------------------------------------------------------- */ -/* Try to find in the existing list of loaded grids. Add all */ -/* matching grids as with NTv2 we can get many grids from one */ -/* file (one shared gridname). */ -/* -------------------------------------------------------------------- */ - for( this_grid = grid_list; this_grid != nullptr; this_grid = this_grid->next) - { - if( strcmp(this_grid->gridname,gridname) == 0 ) - { - got_match = 1; - - /* don't add to the list if it is invalid. */ - if( this_grid->ct == nullptr ) - return 0; - - /* do we need to grow the list? */ - if( *p_gridcount >= *p_gridmax - 2 ) - { - PJ_GRIDINFO **new_list; - int new_max = *p_gridmax + 20; - - new_list = (PJ_GRIDINFO **) pj_calloc(new_max, sizeof(void *)); - if (!new_list) { - pj_ctx_set_errno( ctx, ENOMEM ); - return 0; - } - if( *p_gridlist != nullptr ) - { - memcpy( new_list, *p_gridlist, - sizeof(void *) * (*p_gridmax) ); - pj_dalloc( *p_gridlist ); - } - - *p_gridlist = new_list; - *p_gridmax = new_max; - } - - /* add to the list */ - (*p_gridlist)[(*p_gridcount)++] = this_grid; - (*p_gridlist)[*p_gridcount] = nullptr; - } - - tail = this_grid; - } - - if( got_match ) - return 1; - -/* -------------------------------------------------------------------- */ -/* Try to load the named grid. */ -/* -------------------------------------------------------------------- */ - this_grid = pj_gridinfo_init( ctx, gridname ); - - if( this_grid == nullptr ) - { - return 0; - } - - if( tail != nullptr ) - tail->next = this_grid; - else - grid_list = this_grid; - -/* -------------------------------------------------------------------- */ -/* Recurse to add the grid now that it is loaded. */ -/* -------------------------------------------------------------------- */ - return pj_gridlist_merge_gridfile( ctx, gridname, p_gridlist, - p_gridcount, p_gridmax ); -} - -/************************************************************************/ -/* pj_gridlist_from_nadgrids() */ -/* */ -/* This functions loads the list of grids corresponding to a */ -/* particular nadgrids string into a list, and returns it. The */ -/* list is kept around till a request is made with a different */ -/* string in order to cut down on the string parsing cost, and */ -/* the cost of building the list of tables each time. */ -/************************************************************************/ - -PJ_GRIDINFO **pj_gridlist_from_nadgrids( projCtx ctx, const char *nadgrids, - int *grid_count) - -{ - const char *s; - PJ_GRIDINFO **gridlist = nullptr; - int grid_max = 0; - - pj_errno = 0; - *grid_count = 0; - - pj_acquire_lock(); - -/* -------------------------------------------------------------------- */ -/* Loop processing names out of nadgrids one at a time. */ -/* -------------------------------------------------------------------- */ - for( s = nadgrids; *s != '\0'; ) - { - size_t end_char; - int required = 1; - char name[PJ_MAX_PATH_LENGTH]; - - if( *s == '@' ) - { - required = 0; - s++; - } - - for( end_char = 0; - s[end_char] != '\0' && s[end_char] != ','; - end_char++ ) {} - - if( end_char >= sizeof(name) ) - { - pj_dalloc( gridlist ); - pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); - pj_release_lock(); - return nullptr; - } - - strncpy( name, s, end_char ); - name[end_char] = '\0'; - - s += end_char; - if( *s == ',' ) - s++; - - if( !pj_gridlist_merge_gridfile( ctx, name, &gridlist, grid_count, - &grid_max) - && required ) - { - pj_dalloc( gridlist ); - pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); - pj_release_lock(); - return nullptr; - } - else - pj_errno = 0; - } - - pj_release_lock(); - - return gridlist; -} diff --git a/src/grids.cpp b/src/grids.cpp index 1a18f975..54dae555 100644 --- a/src/grids.cpp +++ b/src/grids.cpp @@ -302,4 +302,582 @@ const VerticalShiftGrid *VerticalShiftGridSet::gridAt(double lon, return nullptr; } +// --------------------------------------------------------------------------- + +HorizontalShiftGrid::HorizontalShiftGrid(int widthIn, int heightIn, + const ExtentAndRes &extentIn) + : m_width(widthIn), m_height(heightIn), m_extent(extentIn) {} + +// --------------------------------------------------------------------------- + +HorizontalShiftGrid::~HorizontalShiftGrid() = default; + +// --------------------------------------------------------------------------- + +HorizontalShiftGridSet::HorizontalShiftGridSet() = default; + +// --------------------------------------------------------------------------- + +HorizontalShiftGridSet::~HorizontalShiftGridSet() = default; + +// --------------------------------------------------------------------------- + +class NullHorizontalShiftGrid : public HorizontalShiftGrid { + + public: + NullHorizontalShiftGrid() : HorizontalShiftGrid(3, 3, globalExtent()) {} + + bool isNullGrid() const override { return true; } + + bool valueAt(int, int, float &lonShift, float &latShift) const override; +}; + +// --------------------------------------------------------------------------- + +bool NullHorizontalShiftGrid::valueAt(int, int, float &lonShift, + float &latShift) const { + lonShift = 0.0f; + latShift = 0.0f; + return true; +} + +// --------------------------------------------------------------------------- + +static double to_double(const void *data) { + double d; + memcpy(&d, data, sizeof(d)); + return d; +} + +// --------------------------------------------------------------------------- + +class NTv1Grid : public HorizontalShiftGrid { + PJ_CONTEXT *m_ctx; + PAFile m_fp; + + NTv1Grid(const NTv1Grid &) = delete; + NTv1Grid &operator=(const NTv1Grid &) = delete; + + public: + NTv1Grid(PJ_CONTEXT *ctx, PAFile fp, int widthIn, int heightIn, + const ExtentAndRes &extentIn) + : HorizontalShiftGrid(widthIn, heightIn, extentIn), m_ctx(ctx), + m_fp(fp) {} + + ~NTv1Grid() override; + + bool valueAt(int, int, float &lonShift, float &latShift) const override; + + static NTv1Grid *open(PJ_CONTEXT *ctx, PAFile fp, + const std::string &filename); +}; + +// --------------------------------------------------------------------------- + +NTv1Grid::~NTv1Grid() { pj_ctx_fclose(m_ctx, m_fp); } + +// --------------------------------------------------------------------------- + +NTv1Grid *NTv1Grid::open(PJ_CONTEXT *ctx, PAFile fp, + const std::string &filename) { + unsigned char header[192]; + + /* -------------------------------------------------------------------- */ + /* Read the header. */ + /* -------------------------------------------------------------------- */ + if (pj_ctx_fread(ctx, header, sizeof(header), 1, fp) != 1) { + pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID); + return nullptr; + } + + /* -------------------------------------------------------------------- */ + /* Regularize fields of interest. */ + /* -------------------------------------------------------------------- */ + if (IS_LSB) { + swap_words(header + 8, sizeof(int), 1); + swap_words(header + 24, sizeof(double), 1); + swap_words(header + 40, sizeof(double), 1); + swap_words(header + 56, sizeof(double), 1); + swap_words(header + 72, sizeof(double), 1); + swap_words(header + 88, sizeof(double), 1); + swap_words(header + 104, sizeof(double), 1); + } + + if (*((int *)(header + 8)) != 12) { + pj_log(ctx, PJ_LOG_ERROR, + "NTv1 grid shift file has wrong record count, corrupt?"); + pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID); + return nullptr; + } + + ExtentAndRes extent; + extent.westLon = -to_double(header + 72) * DEG_TO_RAD; + extent.southLat = to_double(header + 24) * DEG_TO_RAD; + extent.eastLon = -to_double(header + 56) * DEG_TO_RAD; + extent.northLat = to_double(header + 40) * DEG_TO_RAD; + extent.resLon = to_double(header + 104) * DEG_TO_RAD; + extent.resLat = to_double(header + 88) * DEG_TO_RAD; + if (!(fabs(extent.westLon) <= 4 * M_PI && + fabs(extent.eastLon) <= 4 * M_PI && + fabs(extent.northLat) <= M_PI + 1e-5 && + fabs(extent.southLat) <= M_PI + 1e-5 && + extent.westLon < extent.eastLon && + extent.southLat < extent.northLat && extent.resLon > 1e-10 && + extent.resLat > 1e-10)) { + pj_log(ctx, PJ_LOG_ERROR, "Inconsistent georeferencing for %s", + filename.c_str()); + pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID); + return nullptr; + } + const int columns = static_cast<int>( + fabs((extent.eastLon - extent.westLon) / extent.resLon + 0.5) + 1); + const int rows = static_cast<int>( + fabs((extent.northLat - extent.southLat) / extent.resLat + 0.5) + 1); + + return new NTv1Grid(ctx, fp, columns, rows, extent); +} + +// --------------------------------------------------------------------------- + +bool NTv1Grid::valueAt(int x, int y, float &lonShift, float &latShift) const { + assert(x >= 0 && y >= 0 && x < m_width && y < m_height); + + double two_doubles[2]; + // NTv1 is organized from east to west ! + pj_ctx_fseek(m_ctx, m_fp, + 192 + 2 * sizeof(double) * (y * m_width + m_width - 1 - x), + SEEK_SET); + if (pj_ctx_fread(m_ctx, &two_doubles[0], sizeof(two_doubles), 1, m_fp) != + 1) { + pj_ctx_set_errno(m_ctx, PJD_ERR_FAILED_TO_LOAD_GRID); + return false; + } + if (IS_LSB) { + swap_words(&two_doubles[0], sizeof(double), 2); + } + /* convert seconds to radians */ + latShift = static_cast<float>(two_doubles[0] * ((M_PI / 180.0) / 3600.0)); + lonShift = static_cast<float>(two_doubles[1] * ((M_PI / 180.0) / 3600.0)); + + return true; +} + +// --------------------------------------------------------------------------- + +class CTable2Grid : public HorizontalShiftGrid { + PJ_CONTEXT *m_ctx; + PAFile m_fp; + + CTable2Grid(const CTable2Grid &) = delete; + CTable2Grid &operator=(const CTable2Grid &) = delete; + + public: + CTable2Grid(PJ_CONTEXT *ctx, PAFile fp, int widthIn, int heightIn, + const ExtentAndRes &extentIn) + : HorizontalShiftGrid(widthIn, heightIn, extentIn), m_ctx(ctx), + m_fp(fp) {} + + ~CTable2Grid() override; + + bool valueAt(int, int, float &lonShift, float &latShift) const override; + + static CTable2Grid *open(PJ_CONTEXT *ctx, PAFile fp, + const std::string &filename); +}; + +// --------------------------------------------------------------------------- + +CTable2Grid::~CTable2Grid() { pj_ctx_fclose(m_ctx, m_fp); } + +// --------------------------------------------------------------------------- + +CTable2Grid *CTable2Grid::open(PJ_CONTEXT *ctx, PAFile fp, + const std::string &filename) { + unsigned char header[160]; + + /* -------------------------------------------------------------------- */ + /* Read the header. */ + /* -------------------------------------------------------------------- */ + if (pj_ctx_fread(ctx, header, sizeof(header), 1, fp) != 1) { + pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID); + return nullptr; + } + + /* -------------------------------------------------------------------- */ + /* Regularize fields of interest. */ + /* -------------------------------------------------------------------- */ + if (!IS_LSB) { + swap_words(header + 96, sizeof(double), 4); + swap_words(header + 128, sizeof(int), 2); + } + + ExtentAndRes extent; + static_assert(sizeof(extent.westLon) == 8, "wrong sizeof"); + static_assert(sizeof(extent.southLat) == 8, "wrong sizeof"); + static_assert(sizeof(extent.resLon) == 8, "wrong sizeof"); + static_assert(sizeof(extent.resLat) == 8, "wrong sizeof"); + memcpy(&extent.westLon, header + 96, 8); + memcpy(&extent.southLat, header + 104, 8); + memcpy(&extent.resLon, header + 112, 8); + memcpy(&extent.resLat, header + 120, 8); + if (!(fabs(extent.westLon) <= 4 * M_PI && + fabs(extent.southLat) <= M_PI + 1e-5 && extent.resLon > 1e-10 && + extent.resLat > 1e-10)) { + pj_log(ctx, PJ_LOG_ERROR, "Inconsistent georeferencing for %s", + filename.c_str()); + pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID); + return nullptr; + } + int width; + int height; + memcpy(&width, header + 128, 4); + memcpy(&height, header + 132, 4); + if (width <= 0 || height <= 0) { + pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID); + return nullptr; + } + extent.eastLon = extent.westLon + (width - 1) * extent.resLon; + extent.northLat = extent.southLat + (height - 1) * extent.resLon; + + return new CTable2Grid(ctx, fp, width, height, extent); +} + +// --------------------------------------------------------------------------- + +bool CTable2Grid::valueAt(int x, int y, float &lonShift, + float &latShift) const { + assert(x >= 0 && y >= 0 && x < m_width && y < m_height); + + float two_floats[2]; + pj_ctx_fseek(m_ctx, m_fp, 160 + 2 * sizeof(float) * (y * m_width + x), + SEEK_SET); + if (pj_ctx_fread(m_ctx, &two_floats[0], sizeof(two_floats), 1, m_fp) != 1) { + pj_ctx_set_errno(m_ctx, PJD_ERR_FAILED_TO_LOAD_GRID); + return false; + } + if (!IS_LSB) { + swap_words(&two_floats[0], sizeof(float), 2); + } + + latShift = two_floats[1]; + lonShift = two_floats[0]; + + return true; +} + +// --------------------------------------------------------------------------- + +class NTv2GridSet : public HorizontalShiftGridSet { + PJ_CONTEXT *m_ctx; + PAFile m_fp; + + NTv2GridSet(const NTv2GridSet &) = delete; + NTv2GridSet &operator=(const NTv2GridSet &) = delete; + + NTv2GridSet(PJ_CONTEXT *ctx, PAFile fp) : m_ctx(ctx), m_fp(fp) {} + + public: + ~NTv2GridSet() override; + + static std::unique_ptr<NTv2GridSet> open(PJ_CONTEXT *ctx, PAFile fp, + const std::string &filename); +}; + +// --------------------------------------------------------------------------- + +class NTv2Grid : public HorizontalShiftGrid { + friend class NTv2GridSet; + + std::string m_name; + PJ_CONTEXT *m_ctx; // owned by the parent NTv2GridSet + PAFile m_fp; // owned by the parent NTv2GridSet + unsigned long long m_offset; + bool m_mustSwap; + + NTv2Grid(const NTv2Grid &) = delete; + NTv2Grid &operator=(const NTv2Grid &) = delete; + + public: + NTv2Grid(const std::string &nameIn, PJ_CONTEXT *ctx, PAFile fp, + unsigned long long offsetIn, bool mustSwapIn, int widthIn, + int heightIn, const ExtentAndRes &extentIn) + : HorizontalShiftGrid(widthIn, heightIn, extentIn), m_name(nameIn), + m_ctx(ctx), m_fp(fp), m_offset(offsetIn), m_mustSwap(mustSwapIn) {} + + bool valueAt(int, int, float &lonShift, float &latShift) const override; +}; + +// --------------------------------------------------------------------------- + +bool NTv2Grid::valueAt(int x, int y, float &lonShift, float &latShift) const { + assert(x >= 0 && y >= 0 && x < m_width && y < m_height); + + float two_float[2]; + // NTv2 is organized from east to west ! + // there are 4 components: lat shift, lon shift, lat error, lon error + pj_ctx_fseek( + m_ctx, m_fp, + // FIXME when fseek support unsigned long long + static_cast<long>(m_offset + + 4 * sizeof(float) * (y * m_width + m_width - 1 - x)), + SEEK_SET); + if (pj_ctx_fread(m_ctx, &two_float[0], sizeof(two_float), 1, m_fp) != 1) { + pj_ctx_set_errno(m_ctx, PJD_ERR_FAILED_TO_LOAD_GRID); + return false; + } + if (m_mustSwap) { + swap_words(&two_float[0], sizeof(float), 2); + } + /* convert seconds to radians */ + latShift = static_cast<float>(two_float[0] * ((M_PI / 180.0) / 3600.0)); + lonShift = static_cast<float>(two_float[1] * ((M_PI / 180.0) / 3600.0)); + return true; +} + +// --------------------------------------------------------------------------- + +NTv2GridSet::~NTv2GridSet() { pj_ctx_fclose(m_ctx, m_fp); } + +// --------------------------------------------------------------------------- + +std::unique_ptr<NTv2GridSet> NTv2GridSet::open(PJ_CONTEXT *ctx, PAFile fp, + const std::string &filename) { + auto set = std::unique_ptr<NTv2GridSet>(new NTv2GridSet(ctx, fp)); + set->m_name = filename; + set->m_format = "ntv2"; + + char header[11 * 16]; + + /* -------------------------------------------------------------------- */ + /* Read the header. */ + /* -------------------------------------------------------------------- */ + if (pj_ctx_fread(ctx, header, sizeof(header), 1, fp) != 1) { + pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID); + return nullptr; + } + + const bool must_swap = (header[8] == 11) ? !IS_LSB : IS_LSB; + if (must_swap) { + // swap_words( header+8, 4, 1 ); + // swap_words( header+8+16, 4, 1 ); + swap_words(header + 8 + 32, 4, 1); + // swap_words( header+8+7*16, 8, 1 ); + // swap_words( header+8+8*16, 8, 1 ); + // swap_words( header+8+9*16, 8, 1 ); + // swap_words( header+8+10*16, 8, 1 ); + } + + /* -------------------------------------------------------------------- */ + /* Get the subfile count out ... all we really use for now. */ + /* -------------------------------------------------------------------- */ + unsigned int num_subfiles; + memcpy(&num_subfiles, header + 8 + 32, 4); + + std::map<std::string, NTv2Grid *> mapGrids; + + /* ==================================================================== */ + /* Step through the subfiles, creating a grid for each. */ + /* ==================================================================== */ + for (unsigned subfile = 0; subfile < num_subfiles; subfile++) { + // Read header + if (pj_ctx_fread(ctx, header, sizeof(header), 1, fp) != 1) { + pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID); + return nullptr; + } + + if (strncmp(header, "SUB_NAME", 8) != 0) { + pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID); + return nullptr; + } + + // Byte swap interesting fields if needed. + if (must_swap) { + swap_words(header + 8 + 16 * 4, sizeof(double), 6); + swap_words(header + 8 + 16 * 10, sizeof(int), 1); + } + + std::string gridName; + gridName.append(header + 8, 8); + + ExtentAndRes extent; + extent.westLon = + -to_double(header + 7 * 16 + 8) * DEG_TO_RAD / 3600.0; /* W_LONG */ + extent.southLat = + to_double(header + 4 * 16 + 8) * DEG_TO_RAD / 3600.0; /* S_LAT */ + extent.eastLon = + -to_double(header + 6 * 16 + 8) * DEG_TO_RAD / 3600.0; /* E_LONG */ + extent.northLat = + to_double(header + 5 * 16 + 8) * DEG_TO_RAD / 3600.0; /* N_LAT */ + extent.resLon = to_double(header + 9 * 16 + 8) * DEG_TO_RAD / 3600.0; + extent.resLat = to_double(header + 8 * 16 + 8) * DEG_TO_RAD / 3600.0; + + if (!(fabs(extent.westLon) <= 4 * M_PI && + fabs(extent.eastLon) <= 4 * M_PI && + fabs(extent.northLat) <= M_PI + 1e-5 && + fabs(extent.southLat) <= M_PI + 1e-5 && + extent.westLon < extent.eastLon && + extent.southLat < extent.northLat && extent.resLon > 1e-10 && + extent.resLat > 1e-10)) { + pj_log(ctx, PJ_LOG_ERROR, "Inconsistent georeferencing for %s", + filename.c_str()); + pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID); + return nullptr; + } + const int columns = static_cast<int>( + fabs((extent.eastLon - extent.westLon) / extent.resLon + 0.5) + 1); + const int rows = static_cast<int>( + fabs((extent.northLat - extent.southLat) / extent.resLat + 0.5) + + 1); + + pj_log(ctx, PJ_LOG_DEBUG_MINOR, + "NTv2 %s %dx%d: LL=(%.9g,%.9g) UR=(%.9g,%.9g)", gridName.c_str(), + columns, rows, extent.westLon * RAD_TO_DEG, + extent.southLat * RAD_TO_DEG, extent.eastLon * RAD_TO_DEG, + extent.northLat * RAD_TO_DEG); + + unsigned int gs_count; + memcpy(&gs_count, header + 8 + 16 * 10, 4); + if (gs_count / columns != static_cast<unsigned>(rows)) { + pj_log(ctx, PJ_LOG_ERROR, + "GS_COUNT(%u) does not match expected cells (%dx%d)", + gs_count, columns, rows); + pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID); + return nullptr; + } + + auto offset = pj_ctx_ftell(ctx, fp); + auto grid = std::unique_ptr<NTv2Grid>(new NTv2Grid( + gridName, ctx, fp, offset, must_swap, columns, rows, extent)); + std::string parentName; + parentName.assign(header + 24, 8); + auto iter = mapGrids.find(parentName); + auto gridPtr = grid.get(); + if (iter == mapGrids.end()) { + set->m_grids.emplace_back(std::move(grid)); + } else { + iter->second->m_children.emplace_back(std::move(grid)); + } + mapGrids[gridName] = gridPtr; + + // Skip grid data. 4 components of size float + pj_ctx_fseek(ctx, fp, gs_count * 4 * 4, SEEK_CUR); + } + return set; +} + +// --------------------------------------------------------------------------- + +std::unique_ptr<HorizontalShiftGridSet> +HorizontalShiftGridSet::open(PJ_CONTEXT *ctx, const std::string &filename) { + if (filename == "null") { + auto set = std::unique_ptr<HorizontalShiftGridSet>( + new HorizontalShiftGridSet()); + set->m_name = filename; + set->m_format = "null"; + set->m_grids.push_back(std::unique_ptr<NullHorizontalShiftGrid>( + new NullHorizontalShiftGrid())); + return set; + } + + PAFile fp; + if (!(fp = pj_open_lib(ctx, filename.c_str(), "rb"))) { + ctx->last_errno = 0; /* don't treat as a persistent error */ + return nullptr; + } + + char header[160]; + /* -------------------------------------------------------------------- */ + /* Load a header, to determine the file type. */ + /* -------------------------------------------------------------------- */ + size_t header_size; + if ((header_size = pj_ctx_fread(ctx, header, 1, sizeof(header), fp)) != + sizeof(header)) { + /* some files may be smaller that sizeof(header), eg 160, so */ + ctx->last_errno = 0; /* don't treat as a persistent error */ + pj_log(ctx, PJ_LOG_DEBUG_MAJOR, + "pj_gridinfo_init: short header read of %d bytes", + (int)header_size); + } + + pj_ctx_fseek(ctx, fp, SEEK_SET, 0); + + /* -------------------------------------------------------------------- */ + /* Determine file type. */ + /* -------------------------------------------------------------------- */ + if (header_size >= 144 + 16 && strncmp(header + 0, "HEADER", 6) == 0 && + strncmp(header + 96, "W GRID", 6) == 0 && + strncmp(header + 144, "TO NAD83 ", 16) == 0) { + auto grid = NTv1Grid::open(ctx, fp, filename); + if (!grid) { + pj_ctx_fclose(ctx, fp); + return nullptr; + } + auto set = std::unique_ptr<HorizontalShiftGridSet>( + new HorizontalShiftGridSet()); + set->m_name = filename; + set->m_format = "ntv1"; + set->m_grids.push_back(std::unique_ptr<HorizontalShiftGrid>(grid)); + return set; + } else if (header_size >= 9 && strncmp(header + 0, "CTABLE V2", 9) == 0) { + auto grid = CTable2Grid::open(ctx, fp, filename); + if (!grid) { + pj_ctx_fclose(ctx, fp); + return nullptr; + } + auto set = std::unique_ptr<HorizontalShiftGridSet>( + new HorizontalShiftGridSet()); + set->m_name = filename; + set->m_format = "ctable2"; + set->m_grids.push_back(std::unique_ptr<HorizontalShiftGrid>(grid)); + return set; + } else if (header_size >= 48 + 7 && + strncmp(header + 0, "NUM_OREC", 8) == 0 && + strncmp(header + 48, "GS_TYPE", 7) == 0) { + return NTv2GridSet::open(ctx, fp, filename); + } + + pj_log(ctx, PJ_LOG_DEBUG_MAJOR, "Unrecognized horizontal grid format"); + pj_ctx_fclose(ctx, fp); + return nullptr; +} + +// --------------------------------------------------------------------------- + +const HorizontalShiftGrid *HorizontalShiftGrid::gridAt(double lon, + double lat) const { + for (const auto &child : m_children) { + const auto &extentChild = child->extentAndRes(); + const double epsilon = + (extentChild.resLon + extentChild.resLat) / 10000.0; + if ((extentChild.fullWorldLongitude() || + (lon + epsilon >= extentChild.westLon && + lon - epsilon <= extentChild.eastLon)) && + lat + epsilon >= extentChild.southLat && + lat - epsilon <= extentChild.northLat) { + return child->gridAt(lon, lat); + } + } + return this; +} +// --------------------------------------------------------------------------- + +const HorizontalShiftGrid *HorizontalShiftGridSet::gridAt(double lon, + double lat) const { + for (const auto &grid : m_grids) { + if (dynamic_cast<NullHorizontalShiftGrid *>(grid.get())) { + return grid.get(); + } + const auto &extent = grid->extentAndRes(); + const double epsilon = (extent.resLon + extent.resLat) / 10000.0; + if ((extent.fullWorldLongitude() || + (lon + epsilon >= extent.westLon && + lon - epsilon <= extent.eastLon)) && + lat + epsilon >= extent.southLat && + lat - epsilon <= extent.northLat) { + return grid->gridAt(lon, lat); + } + } + return nullptr; +} + NS_PROJ_END diff --git a/src/grids.hpp b/src/grids.hpp index 0e3ccffb..0c4af09f 100644 --- a/src/grids.hpp +++ b/src/grids.hpp @@ -98,6 +98,57 @@ class VerticalShiftGridSet { const VerticalShiftGrid *gridAt(double lon, double lat) const; }; +// --------------------------------------------------------------------------- + +class HorizontalShiftGrid { + protected: + int m_width; + int m_height; + ExtentAndRes m_extent; + std::vector<std::unique_ptr<HorizontalShiftGrid>> m_children{}; + + public: + HorizontalShiftGrid(int widthIn, int heightIn, + const ExtentAndRes &extentIn); + virtual ~HorizontalShiftGrid(); + + int width() const { return m_width; } + int height() const { return m_height; } + const ExtentAndRes &extentAndRes() const { return m_extent; } + + const HorizontalShiftGrid *gridAt(double lon, double lat) const; + + virtual bool isNullGrid() const { return false; } + + // x = 0 is western-most column, y = 0 is southern-most line + virtual bool valueAt(int x, int y, float &lonShift, + float &latShift) const = 0; +}; + +// --------------------------------------------------------------------------- + +class HorizontalShiftGridSet { + protected: + std::string m_name{}; + std::string m_format{}; + std::vector<std::unique_ptr<HorizontalShiftGrid>> m_grids{}; + + HorizontalShiftGridSet(); + + public: + virtual ~HorizontalShiftGridSet(); + + static std::unique_ptr<HorizontalShiftGridSet> + open(PJ_CONTEXT *ctx, const std::string &filename); + + const std::string &name() const { return m_name; } + const std::string &format() const { return m_format; } + const std::vector<std::unique_ptr<HorizontalShiftGrid>> &grids() const { + return m_grids; + } + const HorizontalShiftGrid *gridAt(double lon, double lat) const; +}; + NS_PROJ_END #endif // GRIDS_HPP_INCLUDED diff --git a/src/init.cpp b/src/init.cpp index 69b4ec5f..19fcf47b 100644 --- a/src/init.cpp +++ b/src/init.cpp @@ -646,9 +646,6 @@ pj_init_ctx_with_allow_init_epsg(projCtx ctx, int argc, char **argv, int allow_i PIN->long_wrap_center = 0.0; strcpy( PIN->axis, "enu" ); - PIN->gridlist = nullptr; - PIN->gridlist_count = 0; - /* Set datum parameters. Similarly to +init parameters we want to expand */ /* +datum parameters as late as possible when dealing with pipelines. */ /* otherwise only the first occurrence of +datum will be expanded and that */ diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake index 1cc3b90c..f6112aef 100644 --- a/src/lib_proj.cmake +++ b/src/lib_proj.cmake @@ -237,8 +237,6 @@ set(SRC_LIBPROJ_CORE geocent.cpp geocent.h geodesic.c - gridinfo.cpp - gridlist.cpp init.cpp initcache.cpp internal.cpp @@ -250,9 +248,6 @@ set(SRC_LIBPROJ_CORE mlfn.cpp msfn.cpp mutex.cpp - nad_cvt.cpp - nad_init.cpp - nad_intr.cpp open_lib.cpp param.cpp phi2.cpp diff --git a/src/malloc.cpp b/src/malloc.cpp index 050721e6..b5fe76fc 100644 --- a/src/malloc.cpp +++ b/src/malloc.cpp @@ -225,9 +225,6 @@ PJ *pj_default_destructor (PJ *P, int errlev) { /* Destructor */ pj_dealloc(P->def_spherification); pj_dealloc(P->def_ellps); - /* free grid lists */ - pj_dealloc( P->gridlist ); - /* We used to call pj_dalloc( P->catalog ), but this will leak */ /* memory. The safe way to clear catalog and grid is to call */ /* pj_gc_unloadall(pj_get_default_ctx()); and pj_deallocate_grids(); */ diff --git a/src/nad_cvt.cpp b/src/nad_cvt.cpp deleted file mode 100644 index 79441d0a..00000000 --- a/src/nad_cvt.cpp +++ /dev/null @@ -1,77 +0,0 @@ -#define PJ_LIB__ - -#include <stdio.h> -#include <stdlib.h> - -#include "proj.h" -#include "proj_internal.h" -#include <math.h> - -#define MAX_ITERATIONS 10 -#define TOL 1e-12 - -PJ_LP nad_cvt(PJ_LP in, int inverse, struct CTABLE *ct) { - PJ_LP t, tb,del, dif; - int i = MAX_ITERATIONS; - const double toltol = TOL*TOL; - - if (in.lam == HUGE_VAL) - return in; - - /* normalize input to ll origin */ - tb = in; - tb.lam -= ct->ll.lam; - tb.phi -= ct->ll.phi; - tb.lam = adjlon (tb.lam - M_PI) + M_PI; - - t = nad_intr (tb, ct); - if (t.lam == HUGE_VAL) - return t; - - if (!inverse) { - in.lam -= t.lam; - in.phi += t.phi; - return in; - } - - t.lam = tb.lam + t.lam; - t.phi = tb.phi - t.phi; - - do { - del = nad_intr(t, ct); - - /* This case used to return failure, but I have - changed it to return the first order approximation - of the inverse shift. This avoids cases where the - grid shift *into* this grid came from another grid. - While we aren't returning optimally correct results - I feel a close result in this case is better than - no result. NFW - To demonstrate use -112.5839956 49.4914451 against - the NTv2 grid shift file from Canada. */ - if (del.lam == HUGE_VAL) - break; - - dif.lam = t.lam - del.lam - tb.lam; - dif.phi = t.phi + del.phi - tb.phi; - t.lam -= dif.lam; - t.phi -= dif.phi; - - } while (--i && (dif.lam*dif.lam + dif.phi*dif.phi > toltol)); /* prob. slightly faster than hypot() */ - - if (i==0) { - /* If we had access to a context, this should go through pj_log, and we should set ctx->errno */ - if (getenv ("PROJ_DEBUG")) - fprintf( stderr, "Inverse grid shift iterator failed to converge.\n" ); - t.lam = t.phi = HUGE_VAL; - return t; - } - - /* and again: pj_log and ctx->errno */ - if (del.lam==HUGE_VAL && getenv ("PROJ_DEBUG")) - fprintf (stderr, "Inverse grid shift iteration failed, presumably at grid edge.\nUsing first approximation.\n"); - - in.lam = adjlon (t.lam + ct->ll.lam); - in.phi = t.phi + ct->ll.phi; - return in; -} diff --git a/src/nad_init.cpp b/src/nad_init.cpp deleted file mode 100644 index 315318be..00000000 --- a/src/nad_init.cpp +++ /dev/null @@ -1,308 +0,0 @@ -/****************************************************************************** - * Project: PROJ.4 - * Purpose: Load datum shift files into memory. - * Author: Frank Warmerdam, warmerdam@pobox.com - * - ****************************************************************************** - * Copyright (c) 2000, Frank Warmerdam - * - * 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 <errno.h> -#include <stddef.h> -#include <stdio.h> -#include <stdlib.h> -#include <string.h> - -#include "proj.h" -#include "proj_internal.h" - -/************************************************************************/ -/* swap_words() */ -/* */ -/* Convert the byte order of the given word(s) in place. */ -/************************************************************************/ - -static const int byte_order_test = 1; -#define IS_LSB (((const unsigned char *) (&byte_order_test))[0] == 1) - -static void swap_words( void *data_in, int word_size, int word_count ) - -{ - int word; - unsigned char *data = (unsigned char *) data_in; - - for( word = 0; word < word_count; word++ ) - { - int i; - - for( i = 0; i < word_size/2; i++ ) - { - unsigned char t; - - t = data[i]; - data[i] = data[word_size-i-1]; - data[word_size-i-1] = t; - } - - data += word_size; - } -} - -/************************************************************************/ -/* nad_ctable_load() */ -/* */ -/* Load the data portion of a ctable formatted grid. */ -/************************************************************************/ - -int nad_ctable_load( projCtx ctx, struct CTABLE *ct, struct projFileAPI_t* fileapi ) - -{ - PAFile fid = (PAFile)fileapi; - size_t a_size; - - pj_ctx_fseek( ctx, fid, sizeof(struct CTABLE), SEEK_SET ); - - /* read all the actual shift values */ - a_size = ct->lim.lam * ct->lim.phi; - ct->cvs = (FLP *) pj_malloc(sizeof(FLP) * a_size); - if( ct->cvs == nullptr - || pj_ctx_fread(ctx, ct->cvs, sizeof(FLP), a_size, fid) != a_size ) - { - pj_dalloc( ct->cvs ); - ct->cvs = nullptr; - - pj_log( ctx, PJ_LOG_ERROR, - "ctable loading failed on fread() - binary incompatible?" ); - pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); - return 0; - } - - return 1; -} - -/************************************************************************/ -/* nad_ctable_init() */ -/* */ -/* Read the header portion of a "ctable" format grid. */ -/************************************************************************/ - -struct CTABLE *nad_ctable_init( projCtx ctx, struct projFileAPI_t* fileapi ) -{ - PAFile fid = (PAFile)fileapi; - struct CTABLE *ct; - int id_end; - - /* read the table header */ - ct = (struct CTABLE *) pj_malloc(sizeof(struct CTABLE)); - if( ct == nullptr - || pj_ctx_fread( ctx, ct, sizeof(struct CTABLE), 1, fid ) != 1 ) - { - pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); - pj_dalloc( ct ); - return nullptr; - } - - /* do some minimal validation to ensure the structure isn't corrupt */ - if( ct->lim.lam < 1 || ct->lim.lam > 100000 - || ct->lim.phi < 1 || ct->lim.phi > 100000 ) - { - pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); - pj_dalloc( ct ); - return nullptr; - } - - /* trim white space and newlines off id */ - for( id_end = (int)strlen(ct->id)-1; id_end > 0; id_end-- ) - { - if( ct->id[id_end] == '\n' || ct->id[id_end] == ' ' ) - ct->id[id_end] = '\0'; - else - break; - } - - ct->cvs = nullptr; - - return ct; -} - -/************************************************************************/ -/* nad_ctable2_load() */ -/* */ -/* Load the data portion of a ctable2 formatted grid. */ -/************************************************************************/ - -int nad_ctable2_load( projCtx ctx, struct CTABLE *ct, struct projFileAPI_t* fileapi ) - -{ - PAFile fid = (PAFile)fileapi; - size_t a_size; - - pj_ctx_fseek( ctx, fid, 160, SEEK_SET ); - - /* read all the actual shift values */ - a_size = ct->lim.lam * ct->lim.phi; - ct->cvs = (FLP *) pj_malloc(sizeof(FLP) * a_size); - if( ct->cvs == nullptr - || pj_ctx_fread(ctx, ct->cvs, sizeof(FLP), a_size, fid) != a_size ) - { - pj_dalloc( ct->cvs ); - ct->cvs = nullptr; - - if( getenv("PROJ_DEBUG") != nullptr ) - { - fprintf( stderr, - "ctable2 loading failed on fread() - binary incompatible?\n" ); - } - - pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); - return 0; - } - - if( !IS_LSB ) - { - swap_words( ct->cvs, 4, (int)a_size * 2 ); - } - - return 1; -} - -/************************************************************************/ -/* nad_ctable2_init() */ -/* */ -/* Read the header portion of a "ctable2" format grid. */ -/************************************************************************/ - -struct CTABLE *nad_ctable2_init( projCtx ctx, struct projFileAPI_t* fileapi ) -{ - PAFile fid = (PAFile)fileapi; - struct CTABLE *ct; - int id_end; - char header[160]; - - if( pj_ctx_fread( ctx, header, sizeof(header), 1, fid ) != 1 ) - { - pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); - return nullptr; - } - - if( !IS_LSB ) - { - swap_words( header + 96, 8, 4 ); - swap_words( header + 128, 4, 2 ); - } - - if( strncmp(header,"CTABLE V2",9) != 0 ) - { - pj_log( ctx, PJ_LOG_ERROR, "ctable2 - wrong header!" ); - pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); - return nullptr; - } - - /* read the table header */ - ct = (struct CTABLE *) pj_malloc(sizeof(struct CTABLE)); - if( ct == nullptr ) - { - pj_ctx_set_errno( ctx, ENOMEM ); - return nullptr; - } - - memcpy( ct->id, header + 16, 80 ); - memcpy( &ct->ll.lam, header + 96, 8 ); - memcpy( &ct->ll.phi, header + 104, 8 ); - memcpy( &ct->del.lam, header + 112, 8 ); - memcpy( &ct->del.phi, header + 120, 8 ); - memcpy( &ct->lim.lam, header + 128, 4 ); - memcpy( &ct->lim.phi, header + 132, 4 ); - - /* do some minimal validation to ensure the structure isn't corrupt */ - if( ct->lim.lam < 1 || ct->lim.lam > 100000 - || ct->lim.phi < 1 || ct->lim.phi > 100000 ) - { - pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); - pj_dalloc( ct ); - return nullptr; - } - - /* trim white space and newlines off id */ - for( id_end = (int)strlen(ct->id)-1; id_end > 0; id_end-- ) - { - if( ct->id[id_end] == '\n' || ct->id[id_end] == ' ' ) - ct->id[id_end] = '\0'; - else - break; - } - - ct->cvs = nullptr; - - return ct; -} - -/************************************************************************/ -/* nad_init() */ -/* */ -/* Read a datum shift file in any of the supported binary formats. */ -/************************************************************************/ - -struct CTABLE *nad_init(projCtx ctx, char *name) -{ - struct CTABLE *ct; - PAFile fid; - - ctx->last_errno = 0; - -/* -------------------------------------------------------------------- */ -/* Open the file using the usual search rules. */ -/* -------------------------------------------------------------------- */ - if (!(fid = pj_open_lib(ctx, name, "rb"))) { - return nullptr; - } - - ct = nad_ctable_init( ctx, (struct projFileAPI_t*)fid ); - if( ct != nullptr ) - { - if( !nad_ctable_load( ctx, ct, (struct projFileAPI_t*)fid ) ) - { - nad_free( ct ); - ct = nullptr; - } - } - - pj_ctx_fclose(ctx, fid); - return ct; -} - -/************************************************************************/ -/* nad_free() */ -/* */ -/* Free a CTABLE grid shift structure produced by nad_init(). */ -/************************************************************************/ - -void nad_free(struct CTABLE *ct) -{ - if (ct) { - if( ct->cvs != nullptr ) - pj_dalloc(ct->cvs); - - pj_dalloc(ct); - } -} diff --git a/src/nad_intr.cpp b/src/nad_intr.cpp deleted file mode 100644 index 36ab0a99..00000000 --- a/src/nad_intr.cpp +++ /dev/null @@ -1,67 +0,0 @@ -/* Determine nad table correction value */ -#define PJ_LIB__ -#include "proj.h" -#include "proj_internal.h" -#include <math.h> - -PJ_LP nad_intr(PJ_LP t, struct CTABLE *ct) { - PJ_LP val, frct; - ILP indx; - double m00, m10, m01, m11; - FLP *f00, *f10, *f01, *f11; - long index; - int in; - - t.lam /= ct->del.lam; - indx.lam = isnan(t.lam) ? 0 : (pj_int32)lround(floor(t.lam)); - t.phi /= ct->del.phi; - indx.phi = isnan(t.phi) ? 0 : (pj_int32)lround(floor(t.phi)); - - frct.lam = t.lam - indx.lam; - frct.phi = t.phi - indx.phi; - val.lam = val.phi = HUGE_VAL; - if (indx.lam < 0) { - if (indx.lam == -1 && frct.lam > 0.99999999999) { - ++indx.lam; - frct.lam = 0.; - } else - return val; - } else if ((in = indx.lam + 1) >= ct->lim.lam) { - if (in == ct->lim.lam && frct.lam < 1e-11) { - --indx.lam; - frct.lam = 1.; - } else - return val; - } - if (indx.phi < 0) { - if (indx.phi == -1 && frct.phi > 0.99999999999) { - ++indx.phi; - frct.phi = 0.; - } else - return val; - } else if ((in = indx.phi + 1) >= ct->lim.phi) { - if (in == ct->lim.phi && frct.phi < 1e-11) { - --indx.phi; - frct.phi = 1.; - } else - return val; - } - index = indx.phi * ct->lim.lam + indx.lam; - f00 = ct->cvs + index++; - f10 = ct->cvs + index; - index += ct->lim.lam; - f11 = ct->cvs + index--; - f01 = ct->cvs + index; - m11 = m10 = frct.lam; - m00 = m01 = 1. - frct.lam; - m11 *= frct.phi; - m01 *= frct.phi; - frct.phi = 1. - frct.phi; - m00 *= frct.phi; - m10 *= frct.phi; - val.lam = m00 * f00->lam + m10 * f10->lam + - m01 * f01->lam + m11 * f11->lam; - val.phi = m00 * f00->phi + m10 * f10->phi + - m01 * f01->phi + m11 * f11->phi; - return val; -} diff --git a/src/proj_internal.h b/src/proj_internal.h index 1e43e066..b7a84f87 100644 --- a/src/proj_internal.h +++ b/src/proj_internal.h @@ -479,8 +479,7 @@ struct PJconsts { int datum_type = PJD_UNKNOWN; /* PJD_UNKNOWN/3PARAM/7PARAM/GRIDSHIFT/WGS84 */ double datum_params[7] = {0,0,0,0,0,0,0}; /* Parameters for 3PARAM and 7PARAM */ - struct _pj_gi **gridlist = nullptr; /* TODO: Description needed */ - int gridlist_count = 0; + std::vector<std::unique_ptr<NS_PROJ::HorizontalShiftGridSet>> hgrids{}; int has_geoid_vgrids = 0; /* TODO: Description needed */ std::vector<std::unique_ptr<NS_PROJ::VerticalShiftGridSet>> vgrids{}; @@ -752,35 +751,6 @@ PJ *pj_projection_specific_setup_##name (PJ *P) #endif /* def PJ_LIB__ */ - -#define MAX_TAB_ID 80 -typedef struct { float lam, phi; } FLP; -typedef struct { pj_int32 lam, phi; } ILP; - -struct CTABLE { - char id[MAX_TAB_ID]; /* ascii info */ - PJ_LP ll; /* lower left corner coordinates */ - PJ_LP del; /* size of cells */ - ILP lim; /* limits of conversion matrix */ - FLP *cvs; /* conversion matrix */ -}; - -typedef struct _pj_gi { - char *gridname; /* identifying name of grid, eg "conus" or ntv2_0.gsb */ - char *filename; /* full path to filename */ - - const char *format; /* format of this grid, ie "ctable", "ntv1", - "ntv2" or "missing". */ - - long grid_offset; /* offset in file, for delayed loading */ - int must_swap; /* only for NTv2 */ - - struct CTABLE *ct; - - struct _pj_gi *next; - struct _pj_gi *child; -} PJ_GRIDINFO; - /* procedure prototypes */ double PROJ_DLL dmstor(const char *, char **); double dmstor_ctx(projCtx_t *ctx, const char *, char **); @@ -827,32 +797,6 @@ COMPLEX pj_zpolyd1(COMPLEX, const COMPLEX *, int, COMPLEX *); int pj_deriv(PJ_LP, double, const PJ *, struct DERIVS *); int pj_factors(PJ_LP, const PJ *, double, struct FACTORS *); -/* nadcon related protos */ -PJ_LP nad_intr(PJ_LP, struct CTABLE *); -PJ_LP nad_cvt(PJ_LP, int, struct CTABLE *); -struct CTABLE *nad_init(projCtx_t *ctx, char *); -struct CTABLE *nad_ctable_init( projCtx_t *ctx, struct projFileAPI_t* fid ); -int nad_ctable_load( projCtx_t *ctx, struct CTABLE *, struct projFileAPI_t* fid ); -struct CTABLE *nad_ctable2_init( projCtx_t *ctx, struct projFileAPI_t* fid ); -int nad_ctable2_load( projCtx_t *ctx, struct CTABLE *, struct projFileAPI_t* fid ); -void nad_free(struct CTABLE *); - -/* higher level handling of datum grid shift files */ - -int pj_apply_gridshift_2( PJ *defn, int inverse, - long point_count, int point_offset, - double *x, double *y, double *z ); -int pj_apply_gridshift_3( projCtx_t *ctx, - PJ_GRIDINFO **gridlist, int gridlist_count, - int inverse, long point_count, int point_offset, - double *x, double *y, double *z ); - -PJ_GRIDINFO **pj_gridlist_from_nadgrids( projCtx_t *, const char *, int * ); - -PJ_GRIDINFO *pj_gridinfo_init( projCtx_t *, const char * ); -int pj_gridinfo_load( projCtx_t *, PJ_GRIDINFO * ); -void pj_gridinfo_free( projCtx_t *, PJ_GRIDINFO * ); - void *proj_mdist_ini(double); double proj_mdist(double, double, double, const void *); double proj_inv_mdist(projCtx_t *ctx, double, const void *); diff --git a/src/transform.cpp b/src/transform.cpp index 863da605..f6002b70 100644 --- a/src/transform.cpp +++ b/src/transform.cpp @@ -867,6 +867,58 @@ int pj_geocentric_from_wgs84( PJ *defn, return 0; } + +/************************************************************************/ +/* 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. */ +/************************************************************************/ +static +int pj_apply_gridshift_2( PJ *defn, int inverse, + long point_count, int point_offset, + double *x, double *y, double * /*z*/ ) + +{ + if( defn->hgrids.empty() ) + { + proj_hgrid_init(defn, "nadgrids"); + } + + for( long i = 0; i < point_count; i++ ) + { + PJ_LP input; + + long io = i * point_offset; + input.phi = y[io]; + input.lam = x[io]; + + auto output = proj_hgrid_apply(defn, input, inverse ? PJ_INV : PJ_FWD); + + if ( output.lam != HUGE_VAL ) + { + y[io] = output.phi; + x[io] = output.lam; + } + else + { + if( defn->ctx->debug_level >= PJ_LOG_DEBUG_MAJOR ) + { + pj_log( defn->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 ); + } + } + } + + return 0; +} + + /************************************************************************/ /* pj_datum_transform() */ /* */ @@ -1107,3 +1159,8 @@ static int adjust_axis( projCtx ctx, return 0; } +// --------------------------------------------------------------------------- + +void pj_deallocate_grids() +{ +} diff --git a/src/transformations/hgridshift.cpp b/src/transformations/hgridshift.cpp index 90633939..465bfe98 100644 --- a/src/transformations/hgridshift.cpp +++ b/src/transformations/hgridshift.cpp @@ -20,7 +20,7 @@ static PJ_XYZ forward_3d(PJ_LPZ lpz, PJ *P) { PJ_COORD point = {{0,0,0,0}}; point.lpz = lpz; - if (P->gridlist != nullptr) { + if (!P->hgrids.empty()) { /* Only try the gridshift if at least one grid is loaded, * otherwise just pass the coordinate through unchanged. */ point.lp = proj_hgrid_apply(P, point.lp, PJ_FWD); @@ -34,7 +34,7 @@ static PJ_LPZ reverse_3d(PJ_XYZ xyz, PJ *P) { PJ_COORD point = {{0,0,0,0}}; point.xyz = xyz; - if (P->gridlist != nullptr) { + if (!P->hgrids.empty()) { /* Only try the gridshift if at least one grid is loaded, * otherwise just pass the coordinate through unchanged. */ point.lp = proj_hgrid_apply(P, point.lp, PJ_INV); |