diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/4D_api.cpp | 11 | ||||
| -rw-r--r-- | src/Makefile.am | 4 | ||||
| -rw-r--r-- | src/apply_gridshift.cpp | 397 | ||||
| -rw-r--r-- | src/apply_vgridshift.cpp | 217 | ||||
| -rw-r--r-- | src/filemanager.cpp | 595 | ||||
| -rw-r--r-- | src/filemanager.hpp | 9 | ||||
| -rw-r--r-- | src/grids.cpp | 787 | ||||
| -rw-r--r-- | src/grids.hpp | 31 | ||||
| -rw-r--r-- | src/lib_proj.cmake | 2 | ||||
| -rw-r--r-- | src/open_lib.cpp | 38 | ||||
| -rw-r--r-- | src/proj.h | 12 | ||||
| -rw-r--r-- | src/proj_internal.h | 4 | ||||
| -rw-r--r-- | src/sqlite3.cpp | 11 | ||||
| -rw-r--r-- | src/transform.cpp | 8 | ||||
| -rw-r--r-- | src/transformations/deformation.cpp | 86 | ||||
| -rw-r--r-- | src/transformations/hgridshift.cpp | 10 | ||||
| -rw-r--r-- | src/transformations/vgridshift.cpp | 10 | ||||
| -rw-r--r-- | src/transformations/xyzgridshift.cpp | 80 |
18 files changed, 1444 insertions, 868 deletions
diff --git a/src/4D_api.cpp b/src/4D_api.cpp index cee8262e..9107723d 100644 --- a/src/4D_api.cpp +++ b/src/4D_api.cpp @@ -48,6 +48,7 @@ #include <math.h> #include "geodesic.h" #include "grids.hpp" +#include "filemanager.hpp" #include "proj/common.hpp" #include "proj/coordinateoperation.hpp" @@ -1450,6 +1451,16 @@ PJ_INFO proj_info (void) { /* build search path string */ auto ctx = pj_get_default_ctx(); if (!ctx || ctx->search_paths.empty()) { + // Env var mostly for testing purposes and being independent from + // an existing installation + const char* ignoreUserWritableDirectory = + getenv("PROJ_IGNORE_USER_WRITABLE_DIRECTORY"); + if( ignoreUserWritableDirectory == nullptr || + ignoreUserWritableDirectory[0] == '\0' ) { + buf = path_append(buf, + pj_context_get_user_writable_directory(ctx, false).c_str(), + &buf_size); + } const char *envPROJ_LIB = getenv("PROJ_LIB"); buf = path_append(buf, envPROJ_LIB, &buf_size); #ifdef PROJ_LIB diff --git a/src/Makefile.am b/src/Makefile.am index d29eb976..39667509 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -196,9 +196,9 @@ libproj_la_SOURCES = \ release.cpp gauss.cpp \ fileapi.cpp \ \ - apply_gridshift.cpp datums.cpp datum_set.cpp transform.cpp \ + datums.cpp datum_set.cpp transform.cpp \ geocent.cpp geocent.h utils.cpp \ - jniproj.cpp mutex.cpp initcache.cpp apply_vgridshift.cpp geodesic.c \ + jniproj.cpp mutex.cpp initcache.cpp geodesic.c \ strtod.cpp \ \ 4D_api.cpp pipeline.cpp \ diff --git a/src/apply_gridshift.cpp b/src/apply_gridshift.cpp deleted file mode 100644 index 4ef86fc0..00000000 --- a/src/apply_gridshift.cpp +++ /dev/null @@ -1,397 +0,0 @@ -/****************************************************************************** - * Project: PROJ.4 - * Purpose: Apply datum shifts based on grid shift files (normally NAD27 to - * NAD83 or the reverse). This module is responsible for keeping - * a list of loaded grids, and calling with each one that is - * allowed for a given datum (expressed as the nadgrids= parameter). - * Author: Frank Warmerdam, warmerdam@pobox.com - * - ****************************************************************************** - * Copyright (c) 2000, Frank Warmerdam <warmerdam@pobox.com> - * - * Permission is hereby granted, free of charge, to any person obtaining a - * copy of this software and associated documentation files (the "Software"), - * to deal in the Software without restriction, including without limitation - * the rights to use, copy, modify, merge, publish, distribute, sublicense, - * and/or sell copies of the Software, and to permit persons to whom the - * Software is furnished to do so, subject to the following conditions: - * - * The above copyright notice and this permission notice shall be included - * in all copies or substantial portions of the Software. - * - * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS - * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, - * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL - * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER - * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING - * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER - * DEALINGS IN THE SOFTWARE. - *****************************************************************************/ - -#ifndef FROM_PROJ_CPP -#define FROM_PROJ_CPP -#endif - -#include <math.h> -#include <stdio.h> -#include <string.h> - -#include <limits> - -#include "proj.h" -#include "proj_internal.h" -#include "proj/internal/internal.hpp" -#include "grids.hpp" - -NS_PROJ_START - -// --------------------------------------------------------------------------- - -static const HorizontalShiftGrid* findGrid(const ListOfHGrids& grids, - PJ_LP input) -{ - for( const auto& gridset: grids ) - { - auto grid = gridset->gridAt(input.lam, input.phi); - if( grid ) - return grid; - } - return nullptr; -} - -// --------------------------------------------------------------------------- - -static ListOfHGrids getListOfGridSets(PJ_CONTEXT* ctx, const char* grids) -{ - ListOfHGrids list; - auto listOfGrids = internal::split(std::string(grids), ','); - for( const auto& grid: listOfGrids ) - { - const char* gridname = grid.c_str(); - bool canFail = false; - if( gridname[0] == '@' ) - { - canFail = true; - gridname ++; - } - auto gridSet = HorizontalShiftGridSet::open(ctx, gridname); - if( !gridSet ) - { - if( !canFail ) - { - pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); - return {}; - } - } - else - { - list.emplace_back(std::move(gridSet)); - } - } - return list; -} - - -/**********************************************/ -ListOfHGrids proj_hgrid_init(PJ* P, const char *gridkey) { -/********************************************** - - Initizalize and populate list of horizontal - grids. - - Takes a PJ-object and the plus-parameter - name that is used in the proj-string to - specify the grids to load, e.g. "+grids". - The + should be left out here. - - Returns the number of loaded grids. - -***********************************************/ - - std::string key("s"); - key += gridkey; - const char* grids = pj_param(P->ctx, P->params, key.c_str()).s; - if( grids == nullptr ) - return {}; - - return getListOfGridSets(P->ctx, grids); -} - -typedef struct { pj_int32 lam, phi; } ILP; - -static PJ_LP nad_intr(PJ_LP t, const HorizontalShiftGrid* grid, bool compensateNTConvention) { - 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, compensateNTConvention, f00Lon, f00Lat) || - !grid->valueAt(indx.lam + 1, indx.phi, compensateNTConvention, f10Lon, f10Lat) || - !grid->valueAt(indx.lam, indx.phi + 1, compensateNTConvention, f01Lon, f01Lat) || - !grid->valueAt(indx.lam + 1, indx.phi + 1, compensateNTConvention, 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(projCtx ctx, PJ_LP in, int inverse, const HorizontalShiftGrid* grid, - const ListOfHGrids& grids) { - 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, true); - 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, grid, true); - - /* We can possibly go outside of the initial guessed grid, so try */ - /* to fetch a new grid into which iterate... */ - if (del.lam == HUGE_VAL) - { - PJ_LP lp; - lp.lam = t.lam + extent->westLon; - lp.phi = t.phi + extent->southLat; - auto newGrid = findGrid(grids, lp); - if( newGrid == nullptr || newGrid == grid || newGrid->isNullGrid() ) - break; - pj_log(ctx, PJ_LOG_DEBUG_MINOR, "Switching from grid %s to grid %s", - grid->name().c_str(), - newGrid->name().c_str()); - grid = newGrid; - extent = &(grid->extentAndRes()); - t.lam = lp.lam - extent->westLon; - t.phi = lp.phi - extent->southLat; - tb = in; - tb.lam -= extent->westLon; - tb.phi -= extent->southLat; - tb.lam = adjlon (tb.lam - M_PI) + M_PI; - dif.lam = std::numeric_limits<double>::max(); - dif.phi = std::numeric_limits<double>::max(); - continue; - } - - 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 + extent->westLon); - in.phi = t.phi + extent->southLat; - return in; -} - -/********************************************/ -/* proj_hgrid_value() */ -/* */ -/* Return coordinate offset in grid */ -/********************************************/ -PJ_LP proj_hgrid_value(PJ *P, const ListOfHGrids& grids, PJ_LP lp) { - PJ_LP out = proj_coord_error().lp; - - const auto grid = findGrid(grids, lp); - if( !grid ) { - pj_ctx_set_errno( P->ctx, PJD_ERR_GRID_AREA ); - return out; - } - - /* normalize input to ll origin */ - 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, grid, false); - - if (out.lam == HUGE_VAL || out.phi == HUGE_VAL) { - pj_ctx_set_errno(P->ctx, PJD_ERR_GRID_AREA); - } - - return out; -} - -static -PJ_LP proj_hgrid_apply_internal(PJ_CONTEXT* ctx, - PJ_LP lp, - PJ_DIRECTION direction, - const ListOfHGrids& grids) -{ - PJ_LP out; - - out.lam = HUGE_VAL; - out.phi = HUGE_VAL; - - 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; - } - - int inverse = direction == PJ_FWD ? 0 : 1; - out = nad_cvt(ctx, lp, inverse, grid, grids); - - if (out.lam == HUGE_VAL || out.phi == HUGE_VAL) - pj_ctx_set_errno(ctx, PJD_ERR_GRID_AREA); - - return out; -} - -PJ_LP proj_hgrid_apply(PJ *P, const ListOfHGrids& grids, PJ_LP lp, PJ_DIRECTION direction) { - return proj_hgrid_apply_internal(P->ctx, lp, direction, grids); -} - - -NS_PROJ_END - -/************************************************************************/ -/* 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 = NS_PROJ::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/apply_vgridshift.cpp b/src/apply_vgridshift.cpp deleted file mode 100644 index b0136e5c..00000000 --- a/src/apply_vgridshift.cpp +++ /dev/null @@ -1,217 +0,0 @@ -/****************************************************************************** - * Project: PROJ.4 - * Purpose: Apply vertical datum shifts based on grid shift files, normally - * geoid grids mapping WGS84 to NAVD88 or something similar. - * Author: Frank Warmerdam, warmerdam@pobox.com - * - ****************************************************************************** - * Copyright (c) 2010, 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. - *****************************************************************************/ - -#ifndef FROM_PROJ_CPP -#define FROM_PROJ_CPP -#endif - -#include <assert.h> -#include <stdio.h> -#include <string.h> - -#include <math.h> -#include "proj_internal.h" -#include "proj/internal/internal.hpp" -#include "grids.hpp" - -NS_PROJ_START - -static double read_vgrid_value(const ListOfVGrids& grids, PJ_LP input, double vmultiplier) { - - /* do not deal with NaN coordinates */ - /* cppcheck-suppress duplicateExpression */ - if( isnan(input.phi) || isnan(input.lam) ) - { - return HUGE_VAL; - } - - const VerticalShiftGrid* grid = nullptr; - for( const auto& gridset: grids ) - { - grid = gridset->gridAt(input.lam, input.phi); - if( grid ) - break; - } - if( !grid ) - { - return HUGE_VAL; - } - - const auto& extent = grid->extentAndRes(); - - /* Interpolation a location within the grid */ - double grid_x = (input.lam - extent.westLon) / extent.resLon; - if( extent.fullWorldLongitude() ) { - // The first fmod goes to ]-lim, lim[ range - // So we add lim again to be in ]0, 2*lim[ and fmod again - grid_x = fmod( - fmod(grid_x + grid->width(), grid->width()) + grid->width(), - grid->width()); - } - double grid_y = (input.phi - extent.southLat) / extent.resLat; - int grid_ix = static_cast<int>(lround(floor(grid_x))); - assert(grid_ix >= 0 && grid_ix < grid->width()); - int grid_iy = static_cast<int>(lround(floor(grid_y))); - assert(grid_iy >= 0 && grid_iy < grid->height()); - grid_x -= grid_ix; - grid_y -= grid_iy; - - int grid_ix2 = grid_ix + 1; - if( grid_ix2 >= grid->width() ) { - if( extent.fullWorldLongitude() ) { - grid_ix2 = 0; - } else { - grid_ix2 = grid->width() - 1; - } - } - int grid_iy2 = grid_iy + 1; - if( grid_iy2 >= grid->height() ) - grid_iy2 = grid->height() - 1; - - float value_a = 0; - float value_b = 0; - float value_c = 0; - float value_d = 0; - if( !grid->valueAt(grid_ix, grid_iy, value_a) || - !grid->valueAt(grid_ix2, grid_iy, value_b) || - !grid->valueAt(grid_ix, grid_iy2, value_c) || - !grid->valueAt(grid_ix2, grid_iy2, value_d) ) - { - return HUGE_VAL; - } - - double total_weight = 0.0; - int n_weights = 0; - double value = 0.0f; - - if( !grid->isNodata(value_a, vmultiplier) ) - { - double weight = (1.0-grid_x) * (1.0-grid_y); - value += value_a * weight; - total_weight += weight; - n_weights ++; - } - if( !grid->isNodata(value_b, vmultiplier) ) - { - double weight = (grid_x) * (1.0-grid_y); - value += value_b * weight; - total_weight += weight; - n_weights ++; - } - if( !grid->isNodata(value_c, vmultiplier) ) - { - double weight = (1.0-grid_x) * (grid_y); - value += value_c * weight; - total_weight += weight; - n_weights ++; - } - if( !grid->isNodata(value_d, vmultiplier) ) - { - double weight = (grid_x) * (grid_y); - value += value_d * weight; - total_weight += weight; - n_weights ++; - } - if( n_weights == 0 ) - value = HUGE_VAL; - else if( n_weights != 4 ) - value /= total_weight; - - return value * vmultiplier; -} - -/**********************************************/ -ListOfVGrids proj_vgrid_init(PJ* P, const char *gridkey) { -/********************************************** - - Initizalize and populate gridlist. - - Takes a PJ-object and the plus-parameter - name that is used in the proj-string to - specify the grids to load, e.g. "+grids". - The + should be left out here. - - Returns the number of loaded grids. - -***********************************************/ - - std::string key("s"); - key += gridkey; - const char* gridnames = pj_param(P->ctx, P->params, key.c_str()).s; - if( gridnames == nullptr ) - return {}; - - auto listOfGridNames = internal::split(std::string(gridnames), ','); - ListOfVGrids grids; - for( const auto& gridnameStr: listOfGridNames ) - { - const char* gridname = gridnameStr.c_str(); - bool canFail = false; - if( gridname[0] == '@' ) - { - canFail = true; - gridname ++; - } - auto gridSet = VerticalShiftGridSet::open(P->ctx, gridname); - if( !gridSet ) - { - if( !canFail ) - { - pj_ctx_set_errno( P->ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); - return {}; - } - } - else - { - grids.emplace_back(std::move(gridSet)); - } - } - - return grids; -} - -/***********************************************/ -double proj_vgrid_value(PJ *P, const ListOfVGrids& grids, PJ_LP lp, double vmultiplier){ -/*********************************************** - - Read grid value at position lp in grids loaded - with proj_grid_init. - - Returns the grid value of the given coordinate. - -************************************************/ - - double value; - - value = read_vgrid_value(grids, lp, vmultiplier); - proj_log_trace(P, "proj_vgrid_value: (%f, %f) = %f", lp.lam*RAD_TO_DEG, lp.phi*RAD_TO_DEG, value); - - return value; -} - -NS_PROJ_END diff --git a/src/filemanager.cpp b/src/filemanager.cpp index 68910a94..9acea83e 100644 --- a/src/filemanager.cpp +++ b/src/filemanager.cpp @@ -135,6 +135,9 @@ class FileStdio : public File { unsigned long long tell() override; void reassign_context(PJ_CONTEXT *ctx) override { m_ctx = ctx; } + // We may lie, but the real use case is only for network files + bool hasChanged() const override { return false; } + static std::unique_ptr<File> open(PJ_CONTEXT *ctx, const char *filename); }; @@ -198,6 +201,9 @@ class FileLegacyAdapter : public File { unsigned long long tell() override; void reassign_context(PJ_CONTEXT *ctx) override { m_ctx = ctx; } + // We may lie, but the real use case is only for network files + bool hasChanged() const override { return false; } + static std::unique_ptr<File> open(PJ_CONTEXT *ctx, const char *filename); }; @@ -462,6 +468,13 @@ static const char *cache_db_structure_sql = " lastModified TEXT," " etag TEXT" ");" + "CREATE TABLE downloaded_file_properties(" + " url TEXT PRIMARY KEY NOT NULL," + " lastChecked TIMESTAMP NOT NULL," + " fileSize INTEGER NOT NULL," + " lastModified TEXT," + " etag TEXT" + ");" "CREATE TABLE chunk_data(" " id INTEGER PRIMARY KEY AUTOINCREMENT CHECK (id > 0)," " data BLOB NOT NULL" @@ -672,6 +685,7 @@ void DiskChunkCache::closeAndUnlink() { if (hDB_) { sqlite3_exec(hDB_, "COMMIT", nullptr, nullptr, nullptr); sqlite3_close(hDB_); + hDB_ = nullptr; } if (vfs_) { vfs_->raw()->xDelete(vfs_->raw(), path_.c_str(), 0); @@ -1336,8 +1350,9 @@ class NetworkFile : public File { unsigned long long m_pos = 0; size_t m_nBlocksToDownload = 1; unsigned long long m_lastDownloadedOffset; - unsigned long long m_filesize; + FileProperties m_props; proj_network_close_cbk_type m_closeCbk; + bool m_hasChanged = false; NetworkFile(const NetworkFile &) = delete; NetworkFile &operator=(const NetworkFile &) = delete; @@ -1346,9 +1361,9 @@ class NetworkFile : public File { NetworkFile(PJ_CONTEXT *ctx, const std::string &url, PROJ_NETWORK_HANDLE *handle, unsigned long long lastDownloadOffset, - unsigned long long filesize) + const FileProperties &props) : File(url), m_ctx(ctx), m_url(url), m_handle(handle), - m_lastDownloadedOffset(lastDownloadOffset), m_filesize(filesize), + m_lastDownloadedOffset(lastDownloadOffset), m_props(props), m_closeCbk(ctx->networking.close) {} public: @@ -1358,18 +1373,51 @@ class NetworkFile : public File { bool seek(unsigned long long offset, int whence) override; unsigned long long tell() override; void reassign_context(PJ_CONTEXT *ctx) override; + bool hasChanged() const override { return m_hasChanged; } static std::unique_ptr<File> open(PJ_CONTEXT *ctx, const char *filename); + + static bool get_props_from_headers(PJ_CONTEXT *ctx, + PROJ_NETWORK_HANDLE *handle, + FileProperties &props); }; // --------------------------------------------------------------------------- +bool NetworkFile::get_props_from_headers(PJ_CONTEXT *ctx, + PROJ_NETWORK_HANDLE *handle, + FileProperties &props) { + const char *contentRange = ctx->networking.get_header_value( + ctx, handle, "Content-Range", ctx->networking.user_data); + if (contentRange) { + const char *slash = strchr(contentRange, '/'); + if (slash) { + props.size = std::stoull(slash + 1); + + const char *lastModified = ctx->networking.get_header_value( + ctx, handle, "Last-Modified", ctx->networking.user_data); + if (lastModified) + props.lastModified = lastModified; + + const char *etag = ctx->networking.get_header_value( + ctx, handle, "ETag", ctx->networking.user_data); + if (etag) + props.etag = etag; + + return true; + } + } + return false; +} + +// --------------------------------------------------------------------------- + std::unique_ptr<File> NetworkFile::open(PJ_CONTEXT *ctx, const char *filename) { FileProperties props; if (gNetworkChunkCache.get(ctx, filename, 0, props)) { return std::unique_ptr<File>(new NetworkFile( ctx, filename, nullptr, - std::numeric_limits<unsigned long long>::max(), props.size)); + std::numeric_limits<unsigned long long>::max(), props)); } else { std::vector<unsigned char> buffer(DOWNLOAD_CHUNK_SIZE); size_t size_read = 0; @@ -1386,40 +1434,18 @@ std::unique_ptr<File> NetworkFile::open(PJ_CONTEXT *ctx, const char *filename) { errorBuffer.c_str()); } - unsigned long long filesize = 0; + bool ok = false; if (handle) { - const char *contentRange = ctx->networking.get_header_value( - ctx, handle, "Content-Range", ctx->networking.user_data); - if (contentRange) { - const char *slash = strchr(contentRange, '/'); - if (slash) { - filesize = std::stoull(slash + 1); - - props.size = filesize; - - const char *lastModified = ctx->networking.get_header_value( - ctx, handle, "Last-Modified", - ctx->networking.user_data); - if (lastModified) - props.lastModified = lastModified; - - const char *etag = ctx->networking.get_header_value( - ctx, handle, "ETag", ctx->networking.user_data); - if (etag) - props.etag = etag; - - gNetworkFileProperties.insert(ctx, filename, props); - } - } - if (filesize != 0) { + if (get_props_from_headers(ctx, handle, props)) { + ok = true; + gNetworkFileProperties.insert(ctx, filename, props); gNetworkChunkCache.insert(ctx, filename, 0, std::move(buffer)); } } return std::unique_ptr<File>( - handle != nullptr && filesize != 0 - ? new NetworkFile(ctx, filename, handle, size_read, filesize) - : nullptr); + ok ? new NetworkFile(ctx, filename, handle, size_read, props) + : nullptr); } } @@ -1504,6 +1530,20 @@ size_t NetworkFile::read(void *buffer, size_t sizeBytes) { } return 0; } + + if (!m_hasChanged) { + FileProperties props; + if (get_props_from_headers(m_ctx, m_handle, props)) { + if (props.size != m_props.size || + props.lastModified != m_props.lastModified || + props.etag != m_props.etag) { + gNetworkFileProperties.insert(m_ctx, m_url, props); + gNetworkChunkCache.clearMemoryCache(); + m_hasChanged = true; + } + } + } + region.resize(nRead); m_lastDownloadedOffset = offsetToDownload + nRead; @@ -1546,7 +1586,7 @@ bool NetworkFile::seek(unsigned long long offset, int whence) { } else { if (offset != 0) return false; - m_pos = m_filesize; + m_pos = m_props.size; } return true; } @@ -1870,6 +1910,7 @@ static size_t pj_curl_read_range(PJ_CONTEXT *ctx, auto hCurlHandle = handle->m_handle; double oldDelay = MIN_RETRY_DELAY_MS; + std::string headers; std::string body; char szBuffer[128]; @@ -1879,6 +1920,12 @@ static size_t pj_curl_read_range(PJ_CONTEXT *ctx, while (true) { curl_easy_setopt(hCurlHandle, CURLOPT_RANGE, szBuffer); + headers.clear(); + headers.reserve(16 * 1024); + curl_easy_setopt(hCurlHandle, CURLOPT_HEADERDATA, &headers); + curl_easy_setopt(hCurlHandle, CURLOPT_HEADERFUNCTION, + pj_curl_write_func); + body.clear(); body.reserve(size_to_read); curl_easy_setopt(hCurlHandle, CURLOPT_WRITEDATA, &body); @@ -1929,6 +1976,8 @@ static size_t pj_curl_read_range(PJ_CONTEXT *ctx, if (!body.empty()) { memcpy(buffer, body.data(), std::min(size_to_read, body.size())); } + handle->m_headers = std::move(headers); + return std::min(size_to_read, body.size()); } @@ -2252,50 +2301,468 @@ static void CreateDirectory(const std::string &path) { // --------------------------------------------------------------------------- -std::string pj_context_get_grid_cache_filename(PJ_CONTEXT *ctx) { - pj_load_ini(ctx); - if (!ctx->gridChunkCache.filename.empty()) { - return ctx->gridChunkCache.filename; +std::string pj_context_get_user_writable_directory(PJ_CONTEXT *ctx, + bool create) { + if (ctx->user_writable_directory.empty()) { + // For testing purposes only + const char *env_var_PROJ_USER_WRITABLE_DIRECTORY = + getenv("PROJ_USER_WRITABLE_DIRECTORY"); + if (env_var_PROJ_USER_WRITABLE_DIRECTORY && + env_var_PROJ_USER_WRITABLE_DIRECTORY[0] != '\0') { + ctx->user_writable_directory = env_var_PROJ_USER_WRITABLE_DIRECTORY; + } } - std::string path; + if (ctx->user_writable_directory.empty()) { + std::string path; #ifdef _WIN32 - std::wstring wPath; - wPath.resize(MAX_PATH); - if (SHGetFolderPathW(nullptr, CSIDL_LOCAL_APPDATA, nullptr, 0, &wPath[0]) == - S_OK) { - wPath.resize(wcslen(wPath.data())); - path = WStringToUTF8(wPath); - } else { - const char *local_app_data = getenv("LOCALAPPDATA"); - if (!local_app_data) { - local_app_data = getenv("TEMP"); + std::wstring wPath; + wPath.resize(MAX_PATH); + if (SHGetFolderPathW(nullptr, CSIDL_LOCAL_APPDATA, nullptr, 0, + &wPath[0]) == S_OK) { + wPath.resize(wcslen(wPath.data())); + path = WStringToUTF8(wPath); + } else { + const char *local_app_data = getenv("LOCALAPPDATA"); if (!local_app_data) { - local_app_data = "c:/users"; + local_app_data = getenv("TEMP"); + if (!local_app_data) { + local_app_data = "c:/users"; + } } + path = local_app_data; } - path = local_app_data; - } #else - const char *xdg_data_home = getenv("XDG_DATA_HOME"); - if (xdg_data_home != nullptr) { - path = xdg_data_home; - } else { - const char *home = getenv("HOME"); - if (home) { + const char *xdg_data_home = getenv("XDG_DATA_HOME"); + if (xdg_data_home != nullptr) { + path = xdg_data_home; + } else { + const char *home = getenv("HOME"); + if (home) { #if defined(__MACH__) && defined(__APPLE__) - path = std::string(home) + "/Library/Logs"; + path = std::string(home) + "/Library/Logs"; #else - path = std::string(home) + "/.local/share"; + path = std::string(home) + "/.local/share"; #endif - } else { - path = "/tmp"; + } else { + path = "/tmp"; + } } - } #endif - path += "/proj"; - CreateDirectory(path); + path += "/proj"; + ctx->user_writable_directory = path; + } + if (create) { + CreateDirectory(ctx->user_writable_directory); + } + return ctx->user_writable_directory; +} + +// --------------------------------------------------------------------------- + +std::string pj_context_get_grid_cache_filename(PJ_CONTEXT *ctx) { + pj_load_ini(ctx); + if (!ctx->gridChunkCache.filename.empty()) { + return ctx->gridChunkCache.filename; + } + const std::string path(pj_context_get_user_writable_directory(ctx, true)); ctx->gridChunkCache.filename = path + "/cache.db"; return ctx->gridChunkCache.filename; } +// --------------------------------------------------------------------------- + +#ifdef WIN32 +static const char dir_chars[] = "/\\"; +#else +static const char dir_chars[] = "/"; +#endif + +static bool is_tilde_slash(const char *name) { + return *name == '~' && strchr(dir_chars, name[1]); +} + +static bool is_rel_or_absolute_filename(const char *name) { + return strchr(dir_chars, *name) || + (*name == '.' && strchr(dir_chars, name[1])) || + (!strncmp(name, "..", 2) && strchr(dir_chars, name[2])) || + (name[0] != '\0' && name[1] == ':' && strchr(dir_chars, name[2])); +} + +static std::string build_url(PJ_CONTEXT *ctx, const char *name) { + if (!is_tilde_slash(name) && !is_rel_or_absolute_filename(name) && + !starts_with(name, "http://") && !starts_with(name, "https://")) { + std::string remote_file(pj_context_get_url_endpoint(ctx)); + if (!remote_file.empty()) { + if (remote_file.back() != '/') { + remote_file += '/'; + } + remote_file += name; + auto pos = remote_file.rfind('.'); + if (pos + 4 == remote_file.size()) { + remote_file = remote_file.substr(0, pos) + ".tif"; + } else { + // For example for resource files like 'alaska' + remote_file += ".tif"; + } + } + return remote_file; + } + return name; +} + //! @endcond + +// --------------------------------------------------------------------------- + +/** Return if a file must be downloaded or is already available in the + * PROJ user-writable directory. + * + * The file will be determinted to have to be downloaded if it does not exist + * yet in the user-writable directory, or if it is determined that a more recent + * version exists. To determine if a more recent version exists, PROJ will + * use the "downloaded_file_properties" table of its grid cache database. + * Consequently files manually placed in the user-writable + * directory without using this function would be considered as + * non-existing/obsolete and would be unconditionnaly downloaded again. + * + * This function can only be used if networking is enabled, and either + * the default curl network API or a custom one have been installed. + * + * @param ctx PROJ context, or NULL + * @param url_or_filename URL or filename (without directory component) + * @param ignore_ttl_setting If set to FALSE, PROJ will only check the + * recentness of an already downloaded file, if + * the delay between the last time it has been + * verified and the current time exceeds the TTL + * setting. This can save network accesses. + * If set to TRUE, PROJ will unconditionnally + * check from the server the recentness of the file. + * @return TRUE if the file must be downloaded with proj_download_file() + * @since 7.0 + */ + +int proj_is_download_needed(PJ_CONTEXT *ctx, const char *url_or_filename, + int ignore_ttl_setting) { + if (ctx == nullptr) { + ctx = pj_get_default_ctx(); + } + if (!pj_context_is_network_enabled(ctx)) { + pj_log(ctx, PJ_LOG_ERROR, "Networking capabilities are not enabled"); + return false; + } + + const auto url(build_url(ctx, url_or_filename)); + const char *filename = strrchr(url.c_str(), '/'); + if (filename == nullptr) + return false; + const auto localFilename( + pj_context_get_user_writable_directory(ctx, false) + filename); + + auto f = NS_PROJ::FileManager::open(ctx, localFilename.c_str()); + if (!f) { + return true; + } + f.reset(); + + auto diskCache = NS_PROJ::DiskChunkCache::open(ctx); + if (!diskCache) + return false; + auto stmt = + diskCache->prepare("SELECT lastChecked, fileSize, lastModified, etag " + "FROM downloaded_file_properties WHERE url = ?"); + if (!stmt) + return true; + stmt->bindText(url.c_str()); + if (stmt->execute() != SQLITE_ROW) { + return true; + } + + NS_PROJ::FileProperties cachedProps; + cachedProps.lastChecked = stmt->getInt64(); + cachedProps.size = stmt->getInt64(); + const char *lastModified = stmt->getText(); + cachedProps.lastModified = lastModified ? lastModified : std::string(); + const char *etag = stmt->getText(); + cachedProps.etag = etag ? etag : std::string(); + + if (!ignore_ttl_setting) { + const auto ttl = NS_PROJ::pj_context_get_grid_cache_ttl(ctx); + if (ttl > 0) { + time_t curTime; + time(&curTime); + if (curTime > cachedProps.lastChecked + ttl) { + + unsigned char dummy; + size_t size_read = 0; + std::string errorBuffer; + errorBuffer.resize(1024); + auto handle = ctx->networking.open( + ctx, url.c_str(), 0, 1, &dummy, &size_read, + errorBuffer.size(), &errorBuffer[0], + ctx->networking.user_data); + if (!handle) { + errorBuffer.resize(strlen(errorBuffer.data())); + pj_log(ctx, PJ_LOG_ERROR, "Cannot open %s: %s", url.c_str(), + errorBuffer.c_str()); + return false; + } + NS_PROJ::FileProperties props; + if (!NS_PROJ::NetworkFile::get_props_from_headers(ctx, handle, + props)) { + ctx->networking.close(ctx, handle, + ctx->networking.user_data); + return false; + } + ctx->networking.close(ctx, handle, ctx->networking.user_data); + + if (props.size != cachedProps.size || + props.lastModified != cachedProps.lastModified || + props.etag != cachedProps.etag) { + return true; + } + + stmt = diskCache->prepare( + "UPDATE downloaded_file_properties SET lastChecked = ? " + "WHERE url = ?"); + if (!stmt) + return false; + stmt->bindInt64(curTime); + stmt->bindText(url.c_str()); + if (stmt->execute() != SQLITE_DONE) { + auto hDB = diskCache->handle(); + pj_log(ctx, PJ_LOG_ERROR, "%s", sqlite3_errmsg(hDB)); + return false; + } + } + } + } + + return false; +} + +// --------------------------------------------------------------------------- + +/** Download a file in the PROJ user-writable directory. + * + * The file will only be downloaded if it does not exist yet in the + * user-writable directory, or if it is determined that a more recent + * version exists. To determine if a more recent version exists, PROJ will + * use the "downloaded_file_properties" table of its grid cache database. + * Consequently files manually placed in the user-writable + * directory without using this function would be considered as + * non-existing/obsolete and would be unconditionnaly downloaded again. + * + * This function can only be used if networking is enabled, and either + * the default curl network API or a custom one have been installed. + * + * @param ctx PROJ context, or NULL + * @param url_or_filename URL or filename (without directory component) + * @param ignore_ttl_setting If set to FALSE, PROJ will only check the + * recentness of an already downloaded file, if + * the delay between the last time it has been + * verified and the current time exceeds the TTL + * setting. This can save network accesses. + * If set to TRUE, PROJ will unconditionnally + * check from the server the recentness of the file. + * @param progress_cbk Progress callback, or NULL. + * The passed percentage is in the [0, 1] range. + * The progress callback must return TRUE + * if download must be continued. + * @param user_data User data to provide to the progress callback, or NULL + * @return TRUE if the download was successful (or not needed) + * @since 7.0 + */ + +int proj_download_file(PJ_CONTEXT *ctx, const char *url_or_filename, + int ignore_ttl_setting, + int (*progress_cbk)(PJ_CONTEXT *, double pct, + void *user_data), + void *user_data) { + if (ctx == nullptr) { + ctx = pj_get_default_ctx(); + } + if (!pj_context_is_network_enabled(ctx)) { + pj_log(ctx, PJ_LOG_ERROR, "Networking capabilities are not enabled"); + return false; + } + if (!proj_is_download_needed(ctx, url_or_filename, ignore_ttl_setting)) { + return true; + } + + const auto url(build_url(ctx, url_or_filename)); + const char *filename = strrchr(url.c_str(), '/'); + if (filename == nullptr) + return false; + const auto localFilename(pj_context_get_user_writable_directory(ctx, true) + + filename); + +#ifdef _WIN32 + const int nPID = GetCurrentProcessId(); +#else + const int nPID = getpid(); +#endif + char szUniqueSuffix[128]; + snprintf(szUniqueSuffix, sizeof(szUniqueSuffix), "%d_%p", nPID, &url); + const auto localFilenameTmp(localFilename + szUniqueSuffix); + FILE *f = fopen(localFilenameTmp.c_str(), "wb"); + if (!f) { + pj_log(ctx, PJ_LOG_ERROR, "Cannot create %s", localFilenameTmp.c_str()); + return false; + } + + constexpr size_t FULL_FILE_CHUNK_SIZE = 1024 * 1024; + std::vector<unsigned char> buffer(FULL_FILE_CHUNK_SIZE); + // For testing purposes only + const char *env_var_PROJ_FULL_FILE_CHUNK_SIZE = + getenv("PROJ_FULL_FILE_CHUNK_SIZE"); + if (env_var_PROJ_FULL_FILE_CHUNK_SIZE && + env_var_PROJ_FULL_FILE_CHUNK_SIZE[0] != '\0') { + buffer.resize(atoi(env_var_PROJ_FULL_FILE_CHUNK_SIZE)); + } + size_t size_read = 0; + std::string errorBuffer; + errorBuffer.resize(1024); + auto handle = ctx->networking.open( + ctx, url.c_str(), 0, buffer.size(), &buffer[0], &size_read, + errorBuffer.size(), &errorBuffer[0], ctx->networking.user_data); + if (!handle) { + errorBuffer.resize(strlen(errorBuffer.data())); + pj_log(ctx, PJ_LOG_ERROR, "Cannot open %s: %s", url.c_str(), + errorBuffer.c_str()); + fclose(f); + unlink(localFilenameTmp.c_str()); + return false; + } + + time_t curTime; + time(&curTime); + NS_PROJ::FileProperties props; + if (!NS_PROJ::NetworkFile::get_props_from_headers(ctx, handle, props)) { + ctx->networking.close(ctx, handle, ctx->networking.user_data); + fclose(f); + unlink(localFilenameTmp.c_str()); + return false; + } + + if (size_read < + std::min(static_cast<unsigned long long>(buffer.size()), props.size)) { + pj_log(ctx, PJ_LOG_ERROR, "Did not get as many bytes as expected"); + ctx->networking.close(ctx, handle, ctx->networking.user_data); + fclose(f); + unlink(localFilenameTmp.c_str()); + return false; + } + if (fwrite(buffer.data(), size_read, 1, f) != 1) { + pj_log(ctx, PJ_LOG_ERROR, "Write error"); + ctx->networking.close(ctx, handle, ctx->networking.user_data); + fclose(f); + unlink(localFilenameTmp.c_str()); + return false; + } + + unsigned long long totalDownloaded = size_read; + while (totalDownloaded < props.size) { + if (totalDownloaded + buffer.size() > props.size) { + buffer.resize(static_cast<size_t>(props.size - totalDownloaded)); + } + errorBuffer.resize(1024); + size_read = ctx->networking.read_range( + ctx, handle, totalDownloaded, buffer.size(), &buffer[0], + errorBuffer.size(), &errorBuffer[0], ctx->networking.user_data); + + if (size_read < buffer.size()) { + pj_log(ctx, PJ_LOG_ERROR, "Did not get as many bytes as expected"); + ctx->networking.close(ctx, handle, ctx->networking.user_data); + fclose(f); + unlink(localFilenameTmp.c_str()); + return false; + } + if (fwrite(buffer.data(), size_read, 1, f) != 1) { + pj_log(ctx, PJ_LOG_ERROR, "Write error"); + ctx->networking.close(ctx, handle, ctx->networking.user_data); + fclose(f); + unlink(localFilenameTmp.c_str()); + return false; + } + + totalDownloaded += size_read; + if (progress_cbk && + !progress_cbk(ctx, double(totalDownloaded) / props.size, + user_data)) { + ctx->networking.close(ctx, handle, ctx->networking.user_data); + fclose(f); + unlink(localFilenameTmp.c_str()); + return false; + } + } + + ctx->networking.close(ctx, handle, ctx->networking.user_data); + fclose(f); + + unlink(localFilename.c_str()); + if (rename(localFilenameTmp.c_str(), localFilename.c_str()) != 0) { + pj_log(ctx, PJ_LOG_ERROR, "Cannot rename %s to %s", + localFilenameTmp.c_str(), localFilename.c_str()); + return false; + } + + auto diskCache = NS_PROJ::DiskChunkCache::open(ctx); + if (!diskCache) + return false; + auto stmt = + diskCache->prepare("SELECT lastChecked, fileSize, lastModified, etag " + "FROM downloaded_file_properties WHERE url = ?"); + if (!stmt) + return false; + stmt->bindText(url.c_str()); + + props.lastChecked = curTime; + auto hDB = diskCache->handle(); + + if (stmt->execute() == SQLITE_ROW) { + stmt = diskCache->prepare( + "UPDATE downloaded_file_properties SET lastChecked = ?, " + "fileSize = ?, lastModified = ?, etag = ? " + "WHERE url = ?"); + if (!stmt) + return false; + stmt->bindInt64(props.lastChecked); + stmt->bindInt64(props.size); + if (props.lastModified.empty()) + stmt->bindNull(); + else + stmt->bindText(props.lastModified.c_str()); + if (props.etag.empty()) + stmt->bindNull(); + else + stmt->bindText(props.etag.c_str()); + stmt->bindText(url.c_str()); + if (stmt->execute() != SQLITE_DONE) { + pj_log(ctx, PJ_LOG_ERROR, "%s", sqlite3_errmsg(hDB)); + return false; + } + } else { + stmt = diskCache->prepare( + "INSERT INTO downloaded_file_properties (url, lastChecked, " + "fileSize, lastModified, etag) VALUES " + "(?,?,?,?,?)"); + if (!stmt) + return false; + stmt->bindText(url.c_str()); + stmt->bindInt64(props.lastChecked); + stmt->bindInt64(props.size); + if (props.lastModified.empty()) + stmt->bindNull(); + else + stmt->bindText(props.lastModified.c_str()); + if (props.etag.empty()) + stmt->bindNull(); + else + stmt->bindText(props.etag.c_str()); + if (stmt->execute() != SQLITE_DONE) { + pj_log(ctx, PJ_LOG_ERROR, "%s", sqlite3_errmsg(hDB)); + return false; + } + } + return true; +} diff --git a/src/filemanager.hpp b/src/filemanager.hpp index 993048a7..9793267c 100644 --- a/src/filemanager.hpp +++ b/src/filemanager.hpp @@ -33,10 +33,10 @@ #include "proj.h" #include "proj/util.hpp" -NS_PROJ_START - //! @cond Doxygen_Suppress +NS_PROJ_START + class File; class FileManager { @@ -69,12 +69,13 @@ class File { virtual bool seek(unsigned long long offset, int whence = SEEK_SET) = 0; virtual unsigned long long tell() = 0; virtual void reassign_context(PJ_CONTEXT *ctx) = 0; + virtual bool hasChanged() const = 0; const std::string &name() const { return name_; } }; -//! @endcond Doxygen_Suppress - NS_PROJ_END +//! @endcond Doxygen_Suppress + #endif // FILEMANAGER_HPP_INCLUDED
\ No newline at end of file diff --git a/src/grids.cpp b/src/grids.cpp index 5a99106b..3007fedc 100644 --- a/src/grids.cpp +++ b/src/grids.cpp @@ -140,6 +140,7 @@ class NullVerticalShiftGrid : public VerticalShiftGrid { bool valueAt(int, int, float &out) const override; bool isNodata(float, double) const override { return false; } void reassign_context(PJ_CONTEXT *) override {} + bool hasChanged() const override { return false; } }; // --------------------------------------------------------------------------- @@ -177,6 +178,8 @@ class GTXVerticalShiftGrid : public VerticalShiftGrid { m_ctx = ctx; m_fp->reassign_context(ctx); } + + bool hasChanged() const override { return m_fp->hasChanged(); } }; // --------------------------------------------------------------------------- @@ -373,6 +376,7 @@ class GTiffGrid : public Grid { PJ_CONTEXT *m_ctx; // owned by the belonging GTiffDataset TIFF *m_hTIFF; // owned by the belonging GTiffDataset BlockCache &m_cache; // owned by the belonging GTiffDataset + File *m_fp; // owned by the belonging GTiffDataset uint32 m_ifdIdx; TIFFDataType m_dt; uint16 m_samplesPerPixel; @@ -402,9 +406,9 @@ class GTiffGrid : public Grid { uint32 offsetInBlock, uint16 sample) const; public: - GTiffGrid(PJ_CONTEXT *ctx, TIFF *hTIFF, BlockCache &cache, uint32 ifdIdx, - const std::string &nameIn, int widthIn, int heightIn, - const ExtentAndRes &extentIn, TIFFDataType dtIn, + GTiffGrid(PJ_CONTEXT *ctx, TIFF *hTIFF, BlockCache &cache, File *fp, + uint32 ifdIdx, const std::string &nameIn, int widthIn, + int heightIn, const ExtentAndRes &extentIn, TIFFDataType dtIn, uint16 samplesPerPixelIn, uint16 planarConfig, bool bottomUpIn); ~GTiffGrid() override; @@ -420,17 +424,19 @@ class GTiffGrid : public Grid { uint32 subfileType() const { return m_subfileType; } void reassign_context(PJ_CONTEXT *ctx) { m_ctx = ctx; } + + bool hasChanged() const override { return m_fp->hasChanged(); } }; // --------------------------------------------------------------------------- -GTiffGrid::GTiffGrid(PJ_CONTEXT *ctx, TIFF *hTIFF, BlockCache &cache, +GTiffGrid::GTiffGrid(PJ_CONTEXT *ctx, TIFF *hTIFF, BlockCache &cache, File *fp, uint32 ifdIdx, const std::string &nameIn, int widthIn, int heightIn, const ExtentAndRes &extentIn, TIFFDataType dtIn, uint16 samplesPerPixelIn, uint16 planarConfig, bool bottomUpIn) : Grid(nameIn, widthIn, heightIn, extentIn), m_ctx(ctx), m_hTIFF(hTIFF), - m_cache(cache), m_ifdIdx(ifdIdx), m_dt(dtIn), + m_cache(cache), m_fp(fp), m_ifdIdx(ifdIdx), m_dt(dtIn), m_samplesPerPixel(samplesPerPixelIn), m_planarConfig(planarConfig), m_bottomUp(bottomUpIn), m_dirOffset(TIFFCurrentDirOffset(hTIFF)), m_tiled(TIFFIsTiled(hTIFF) != 0) { @@ -1048,8 +1054,8 @@ std::unique_ptr<GTiffGrid> GTiffDataset::nextGrid() { } auto ret = std::unique_ptr<GTiffGrid>(new GTiffGrid( - m_ctx, m_hTIFF, m_cache, m_ifdIdx, m_filename, width, height, extent, - dt, samplesPerPixel, planarConfig, vRes < 0)); + m_ctx, m_hTIFF, m_cache, m_fp.get(), m_ifdIdx, m_filename, width, + height, extent, dt, samplesPerPixel, planarConfig, vRes < 0)); m_ifdIdx++; m_hasNextGrid = TIFFReadDirectory(m_hTIFF) != 0; m_nextDirOffset = TIFFCurrentDirOffset(m_hTIFF); @@ -1058,10 +1064,12 @@ std::unique_ptr<GTiffGrid> GTiffDataset::nextGrid() { // --------------------------------------------------------------------------- -class GTiffVGridShiftSet : public VerticalShiftGridSet, public GTiffDataset { +class GTiffVGridShiftSet : public VerticalShiftGridSet { + + std::unique_ptr<GTiffDataset> m_GTiffDataset; GTiffVGridShiftSet(PJ_CONTEXT *ctx, std::unique_ptr<File> &&fp) - : GTiffDataset(ctx, std::move(fp)) {} + : m_GTiffDataset(new GTiffDataset(ctx, std::move(fp))) {} public: ~GTiffVGridShiftSet() override; @@ -1072,7 +1080,26 @@ class GTiffVGridShiftSet : public VerticalShiftGridSet, public GTiffDataset { void reassign_context(PJ_CONTEXT *ctx) override { VerticalShiftGridSet::reassign_context(ctx); - GTiffDataset::reassign_context(ctx); + if (m_GTiffDataset) { + m_GTiffDataset->reassign_context(ctx); + } + } + + bool reopen(PJ_CONTEXT *ctx) override { + pj_log(ctx, PJ_LOG_DEBUG_MAJOR, "Grid %s has changed. Re-loading it", + m_name.c_str()); + m_grids.clear(); + m_GTiffDataset.reset(); + auto fp = FileManager::open_resource_file(ctx, m_name.c_str()); + if (!fp) { + return false; + } + auto newGS = open(ctx, std::move(fp), m_name); + if (newGS) { + m_grids = std::move(newGS->m_grids); + m_GTiffDataset = std::move(newGS->m_GTiffDataset); + } + return !m_grids.empty(); } }; @@ -1172,6 +1199,8 @@ class GTiffVGrid : public VerticalShiftGrid { void reassign_context(PJ_CONTEXT *ctx) override { m_grid->reassign_context(ctx); } + + bool hasChanged() const override { return m_grid->hasChanged(); } }; // --------------------------------------------------------------------------- @@ -1221,14 +1250,14 @@ GTiffVGridShiftSet::open(PJ_CONTEXT *ctx, std::unique_ptr<File> fp, new GTiffVGridShiftSet(ctx, std::move(fp))); set->m_name = filename; set->m_format = "gtiff"; - if (!set->openTIFF(filename)) { + if (!set->m_GTiffDataset->openTIFF(filename)) { return nullptr; } uint16 idxSample = 0; std::map<std::string, GTiffVGrid *> mapGrids; for (int ifd = 0;; ++ifd) { - auto grid = set->nextGrid(); + auto grid = set->m_GTiffDataset->nextGrid(); if (!grid) { if (ifd == 0) { return nullptr; @@ -1364,6 +1393,19 @@ VerticalShiftGridSet::open(PJ_CONTEXT *ctx, const std::string &filename) { // --------------------------------------------------------------------------- +bool VerticalShiftGridSet::reopen(PJ_CONTEXT *ctx) { + pj_log(ctx, PJ_LOG_DEBUG_MAJOR, "Grid %s has changed. Re-loading it", + m_name.c_str()); + auto newGS = open(ctx, m_name); + m_grids.clear(); + if (newGS) { + m_grids = std::move(newGS->m_grids); + } + return !m_grids.empty(); +} + +// --------------------------------------------------------------------------- + const VerticalShiftGrid *VerticalShiftGrid::gridAt(double lon, double lat) const { for (const auto &child : m_children) { @@ -1435,6 +1477,8 @@ class NullHorizontalShiftGrid : public HorizontalShiftGrid { float &latShift) const override; void reassign_context(PJ_CONTEXT *) override {} + + bool hasChanged() const override { return false; } }; // --------------------------------------------------------------------------- @@ -1482,6 +1526,8 @@ class NTv1Grid : public HorizontalShiftGrid { m_ctx = ctx; m_fp->reassign_context(ctx); } + + bool hasChanged() const override { return m_fp->hasChanged(); } }; // --------------------------------------------------------------------------- @@ -1603,6 +1649,8 @@ class CTable2Grid : public HorizontalShiftGrid { m_ctx = ctx; m_fp->reassign_context(ctx); } + + bool hasChanged() const override { return m_fp->hasChanged(); } }; // --------------------------------------------------------------------------- @@ -1737,6 +1785,8 @@ class NTv2Grid : public HorizontalShiftGrid { m_ctx = ctx; m_fp->reassign_context(ctx); } + + bool hasChanged() const override { return m_fp->hasChanged(); } }; // --------------------------------------------------------------------------- @@ -1790,6 +1840,11 @@ std::unique_ptr<NTv2GridSet> NTv2GridSet::open(PJ_CONTEXT *ctx, pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID); return nullptr; } + if (memcmp(header + 56, "SECONDS", 7) != 0) { + pj_log(ctx, PJ_LOG_ERROR, "Only GS_TYPE=SECONDS is supported"); + 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) { @@ -1906,10 +1961,12 @@ std::unique_ptr<NTv2GridSet> NTv2GridSet::open(PJ_CONTEXT *ctx, // --------------------------------------------------------------------------- -class GTiffHGridShiftSet : public HorizontalShiftGridSet, public GTiffDataset { +class GTiffHGridShiftSet : public HorizontalShiftGridSet { + + std::unique_ptr<GTiffDataset> m_GTiffDataset; GTiffHGridShiftSet(PJ_CONTEXT *ctx, std::unique_ptr<File> &&fp) - : GTiffDataset(ctx, std::move(fp)) {} + : m_GTiffDataset(new GTiffDataset(ctx, std::move(fp))) {} public: ~GTiffHGridShiftSet() override; @@ -1920,7 +1977,26 @@ class GTiffHGridShiftSet : public HorizontalShiftGridSet, public GTiffDataset { void reassign_context(PJ_CONTEXT *ctx) override { HorizontalShiftGridSet::reassign_context(ctx); - GTiffDataset::reassign_context(ctx); + if (m_GTiffDataset) { + m_GTiffDataset->reassign_context(ctx); + } + } + + bool reopen(PJ_CONTEXT *ctx) override { + pj_log(ctx, PJ_LOG_DEBUG_MAJOR, "Grid %s has changed. Re-loading it", + m_name.c_str()); + m_grids.clear(); + m_GTiffDataset.reset(); + auto fp = FileManager::open_resource_file(ctx, m_name.c_str()); + if (!fp) { + return false; + } + auto newGS = open(ctx, std::move(fp), m_name); + if (newGS) { + m_grids = std::move(newGS->m_grids); + m_GTiffDataset = std::move(newGS->m_GTiffDataset); + } + return !m_grids.empty(); } }; @@ -1954,6 +2030,8 @@ class GTiffHGrid : public HorizontalShiftGrid { void reassign_context(PJ_CONTEXT *ctx) override { m_grid->reassign_context(ctx); } + + bool hasChanged() const override { return m_grid->hasChanged(); } }; // --------------------------------------------------------------------------- @@ -2024,7 +2102,7 @@ GTiffHGridShiftSet::open(PJ_CONTEXT *ctx, std::unique_ptr<File> fp, new GTiffHGridShiftSet(ctx, std::move(fp))); set->m_name = filename; set->m_format = "gtiff"; - if (!set->openTIFF(filename)) { + if (!set->m_GTiffDataset->openTIFF(filename)) { return nullptr; } @@ -2037,7 +2115,7 @@ GTiffHGridShiftSet::open(PJ_CONTEXT *ctx, std::unique_ptr<File> fp, std::map<std::string, GTiffHGrid *> mapGrids; for (int ifd = 0;; ++ifd) { - auto grid = set->nextGrid(); + auto grid = set->m_GTiffDataset->nextGrid(); if (!grid) { if (ifd == 0) { return nullptr; @@ -2269,6 +2347,19 @@ HorizontalShiftGridSet::open(PJ_CONTEXT *ctx, const std::string &filename) { // --------------------------------------------------------------------------- +bool HorizontalShiftGridSet::reopen(PJ_CONTEXT *ctx) { + pj_log(ctx, PJ_LOG_DEBUG_MAJOR, "Grid %s has changed. Re-loading it", + m_name.c_str()); + auto newGS = open(ctx, m_name); + m_grids.clear(); + if (newGS) { + m_grids = std::move(newGS->m_grids); + } + return !m_grids.empty(); +} + +// --------------------------------------------------------------------------- + const HorizontalShiftGrid *HorizontalShiftGrid::gridAt(double lon, double lat) const { for (const auto &child : m_children) { @@ -2317,11 +2408,12 @@ void HorizontalShiftGridSet::reassign_context(PJ_CONTEXT *ctx) { #ifdef TIFF_ENABLED // --------------------------------------------------------------------------- -class GTiffGenericGridShiftSet : public GenericShiftGridSet, - public GTiffDataset { +class GTiffGenericGridShiftSet : public GenericShiftGridSet { + + std::unique_ptr<GTiffDataset> m_GTiffDataset; GTiffGenericGridShiftSet(PJ_CONTEXT *ctx, std::unique_ptr<File> &&fp) - : GTiffDataset(ctx, std::move(fp)) {} + : m_GTiffDataset(new GTiffDataset(ctx, std::move(fp))) {} public: ~GTiffGenericGridShiftSet() override; @@ -2332,7 +2424,26 @@ class GTiffGenericGridShiftSet : public GenericShiftGridSet, void reassign_context(PJ_CONTEXT *ctx) override { GenericShiftGridSet::reassign_context(ctx); - GTiffDataset::reassign_context(ctx); + if (m_GTiffDataset) { + m_GTiffDataset->reassign_context(ctx); + } + } + + bool reopen(PJ_CONTEXT *ctx) override { + pj_log(ctx, PJ_LOG_DEBUG_MAJOR, "Grid %s has changed. Re-loading it", + m_name.c_str()); + m_grids.clear(); + m_GTiffDataset.reset(); + auto fp = FileManager::open_resource_file(ctx, m_name.c_str()); + if (!fp) { + return false; + } + auto newGS = open(ctx, std::move(fp), m_name); + if (newGS) { + m_grids = std::move(newGS->m_grids); + m_GTiffDataset = std::move(newGS->m_GTiffDataset); + } + return !m_grids.empty(); } }; @@ -2375,6 +2486,8 @@ class GTiffGenericGrid : public GenericShiftGrid { void reassign_context(PJ_CONTEXT *ctx) override { m_grid->reassign_context(ctx); } + + bool hasChanged() const override { return m_grid->hasChanged(); } }; // --------------------------------------------------------------------------- @@ -2446,6 +2559,8 @@ class NullGenericShiftGrid : public GenericShiftGrid { } void reassign_context(PJ_CONTEXT *) override {} + + bool hasChanged() const override { return false; } }; // --------------------------------------------------------------------------- @@ -2466,13 +2581,13 @@ GTiffGenericGridShiftSet::open(PJ_CONTEXT *ctx, std::unique_ptr<File> fp, new GTiffGenericGridShiftSet(ctx, std::move(fp))); set->m_name = filename; set->m_format = "gtiff"; - if (!set->openTIFF(filename)) { + if (!set->m_GTiffDataset->openTIFF(filename)) { return nullptr; } std::map<std::string, GTiffGenericGrid *> mapGrids; for (int ifd = 0;; ++ifd) { - auto grid = set->nextGrid(); + auto grid = set->m_GTiffDataset->nextGrid(); if (!grid) { if (ifd == 0) { return nullptr; @@ -2574,6 +2689,19 @@ GenericShiftGridSet::open(PJ_CONTEXT *ctx, const std::string &filename) { // --------------------------------------------------------------------------- +bool GenericShiftGridSet::reopen(PJ_CONTEXT *ctx) { + pj_log(ctx, PJ_LOG_DEBUG_MAJOR, "Grid %s has changed. Re-loading it", + m_name.c_str()); + auto newGS = open(ctx, m_name); + m_grids.clear(); + if (newGS) { + m_grids = std::move(newGS->m_grids); + } + return !m_grids.empty(); +} + +// --------------------------------------------------------------------------- + const GenericShiftGrid *GenericShiftGrid::gridAt(double lon, double lat) const { for (const auto &child : m_children) { const auto &extentChild = child->extentAndRes(); @@ -2614,7 +2742,7 @@ void GenericShiftGridSet::reassign_context(PJ_CONTEXT *ctx) { // --------------------------------------------------------------------------- -ListOfGenericGrids proj_generic_grid_init(PJ *P, const char *gridkey) { +ListOfGenericGrids pj_generic_grid_init(PJ *P, const char *gridkey) { std::string key("s"); key += gridkey; const char *gridnames = pj_param(P->ctx, P->params, key.c_str()).s; @@ -2644,4 +2772,615 @@ ListOfGenericGrids proj_generic_grid_init(PJ *P, const char *gridkey) { return grids; } +// --------------------------------------------------------------------------- + +static const HorizontalShiftGrid * +findGrid(const ListOfHGrids &grids, const PJ_LP &input, + HorizontalShiftGridSet *&gridSetOut) { + for (const auto &gridset : grids) { + auto grid = gridset->gridAt(input.lam, input.phi); + if (grid) { + gridSetOut = gridset.get(); + return grid; + } + } + return nullptr; +} + +// --------------------------------------------------------------------------- + +static ListOfHGrids getListOfGridSets(PJ_CONTEXT *ctx, const char *grids) { + ListOfHGrids list; + auto listOfGrids = internal::split(std::string(grids), ','); + for (const auto &grid : listOfGrids) { + const char *gridname = grid.c_str(); + bool canFail = false; + if (gridname[0] == '@') { + canFail = true; + gridname++; + } + auto gridSet = HorizontalShiftGridSet::open(ctx, gridname); + if (!gridSet) { + if (!canFail) { + pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID); + return {}; + } + } else { + list.emplace_back(std::move(gridSet)); + } + } + return list; +} + +/**********************************************/ +ListOfHGrids pj_hgrid_init(PJ *P, const char *gridkey) { + /********************************************** + + Initizalize and populate list of horizontal + grids. + + Takes a PJ-object and the plus-parameter + name that is used in the proj-string to + specify the grids to load, e.g. "+grids". + The + should be left out here. + + Returns the number of loaded grids. + + ***********************************************/ + + std::string key("s"); + key += gridkey; + const char *grids = pj_param(P->ctx, P->params, key.c_str()).s; + if (grids == nullptr) + return {}; + + return getListOfGridSets(P->ctx, grids); +} + +// --------------------------------------------------------------------------- + +typedef struct { pj_int32 lam, phi; } ILP; + +static PJ_LP pj_hgrid_interpolate(PJ_LP t, const HorizontalShiftGrid *grid, + bool compensateNTConvention) { + PJ_LP val, frct; + ILP indx; + int in; + + const auto &extent = grid->extentAndRes(); + t.lam /= extent.resLon; + indx.lam = std::isnan(t.lam) ? 0 : (pj_int32)lround(floor(t.lam)); + t.phi /= extent.resLat; + indx.phi = std::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, compensateNTConvention, f00Lon, + f00Lat) || + !grid->valueAt(indx.lam + 1, indx.phi, compensateNTConvention, f10Lon, + f10Lat) || + !grid->valueAt(indx.lam, indx.phi + 1, compensateNTConvention, f01Lon, + f01Lat) || + !grid->valueAt(indx.lam + 1, indx.phi + 1, compensateNTConvention, + 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 pj_hgrid_apply_internal(projCtx ctx, PJ_LP in, + PJ_DIRECTION direction, + const HorizontalShiftGrid *grid, + HorizontalShiftGridSet *gridset, + const ListOfHGrids &grids, + bool &shouldRetry) { + PJ_LP t, tb, del, dif; + int i = MAX_ITERATIONS; + const double toltol = TOL * TOL; + + shouldRetry = false; + 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 = pj_hgrid_interpolate(tb, grid, true); + if (grid->hasChanged()) { + shouldRetry = gridset->reopen(ctx); + return t; + } + if (t.lam == HUGE_VAL) + return t; + + if (direction == PJ_FWD) { + in.lam += t.lam; + in.phi += t.phi; + return in; + } + + t.lam = tb.lam - t.lam; + t.phi = tb.phi - t.phi; + + do { + del = pj_hgrid_interpolate(t, grid, true); + if (grid->hasChanged()) { + shouldRetry = gridset->reopen(ctx); + return t; + } + + /* We can possibly go outside of the initial guessed grid, so try */ + /* to fetch a new grid into which iterate... */ + if (del.lam == HUGE_VAL) { + PJ_LP lp; + lp.lam = t.lam + extent->westLon; + lp.phi = t.phi + extent->southLat; + auto newGrid = findGrid(grids, lp, gridset); + if (newGrid == nullptr || newGrid == grid || newGrid->isNullGrid()) + break; + pj_log(ctx, PJ_LOG_DEBUG_MINOR, "Switching from grid %s to grid %s", + grid->name().c_str(), newGrid->name().c_str()); + grid = newGrid; + extent = &(grid->extentAndRes()); + t.lam = lp.lam - extent->westLon; + t.phi = lp.phi - extent->southLat; + tb = in; + tb.lam -= extent->westLon; + tb.phi -= extent->southLat; + tb.lam = adjlon(tb.lam - M_PI) + M_PI; + dif.lam = std::numeric_limits<double>::max(); + dif.phi = std::numeric_limits<double>::max(); + continue; + } + + 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 + extent->westLon); + in.phi = t.phi + extent->southLat; + return in; +} + +// --------------------------------------------------------------------------- + +PJ_LP pj_hgrid_apply(PJ_CONTEXT *ctx, const ListOfHGrids &grids, PJ_LP lp, + PJ_DIRECTION direction) { + PJ_LP out; + + out.lam = HUGE_VAL; + out.phi = HUGE_VAL; + + while (true) { + HorizontalShiftGridSet *gridset = nullptr; + const auto grid = findGrid(grids, lp, gridset); + if (!grid) { + pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID); + return out; + } + if (grid->isNullGrid()) { + return lp; + } + + bool shouldRetry = false; + out = pj_hgrid_apply_internal(ctx, lp, direction, grid, gridset, grids, + shouldRetry); + if (!shouldRetry) { + break; + } + } + + if (out.lam == HUGE_VAL || out.phi == HUGE_VAL) + pj_ctx_set_errno(ctx, PJD_ERR_GRID_AREA); + + return out; +} + +/********************************************/ +/* proj_hgrid_value() */ +/* */ +/* Return coordinate offset in grid */ +/********************************************/ +PJ_LP pj_hgrid_value(PJ *P, const ListOfHGrids &grids, PJ_LP lp) { + PJ_LP out = proj_coord_error().lp; + + HorizontalShiftGridSet *gridset = nullptr; + const auto grid = findGrid(grids, lp, gridset); + if (!grid) { + pj_ctx_set_errno(P->ctx, PJD_ERR_GRID_AREA); + return out; + } + + /* normalize input to ll origin */ + const auto &extent = grid->extentAndRes(); + lp.lam -= extent.westLon; + lp.phi -= extent.southLat; + + lp.lam = adjlon(lp.lam - M_PI) + M_PI; + + out = pj_hgrid_interpolate(lp, grid, false); + if (grid->hasChanged()) { + if (gridset->reopen(P->ctx)) { + return pj_hgrid_value(P, grids, lp); + } + out.lam = HUGE_VAL; + out.phi = HUGE_VAL; + } + + if (out.lam == HUGE_VAL || out.phi == HUGE_VAL) { + pj_ctx_set_errno(P->ctx, PJD_ERR_GRID_AREA); + } + + return out; +} + +// --------------------------------------------------------------------------- + +static double read_vgrid_value(PJ_CONTEXT *ctx, const ListOfVGrids &grids, + const PJ_LP &input, const double vmultiplier) { + + /* do not deal with NaN coordinates */ + /* cppcheck-suppress duplicateExpression */ + if (std::isnan(input.phi) || std::isnan(input.lam)) { + return HUGE_VAL; + } + + VerticalShiftGridSet *curGridset = nullptr; + const VerticalShiftGrid *grid = nullptr; + for (const auto &gridset : grids) { + grid = gridset->gridAt(input.lam, input.phi); + if (grid) { + curGridset = gridset.get(); + break; + } + } + if (!grid) { + return HUGE_VAL; + } + + const auto &extent = grid->extentAndRes(); + + /* Interpolation a location within the grid */ + double grid_x = (input.lam - extent.westLon) / extent.resLon; + if (extent.fullWorldLongitude()) { + // The first fmod goes to ]-lim, lim[ range + // So we add lim again to be in ]0, 2*lim[ and fmod again + grid_x = + fmod(fmod(grid_x + grid->width(), grid->width()) + grid->width(), + grid->width()); + } + double grid_y = (input.phi - extent.southLat) / extent.resLat; + int grid_ix = static_cast<int>(lround(floor(grid_x))); + assert(grid_ix >= 0 && grid_ix < grid->width()); + int grid_iy = static_cast<int>(lround(floor(grid_y))); + assert(grid_iy >= 0 && grid_iy < grid->height()); + grid_x -= grid_ix; + grid_y -= grid_iy; + + int grid_ix2 = grid_ix + 1; + if (grid_ix2 >= grid->width()) { + if (extent.fullWorldLongitude()) { + grid_ix2 = 0; + } else { + grid_ix2 = grid->width() - 1; + } + } + int grid_iy2 = grid_iy + 1; + if (grid_iy2 >= grid->height()) + grid_iy2 = grid->height() - 1; + + float value_a = 0; + float value_b = 0; + float value_c = 0; + float value_d = 0; + bool error = (!grid->valueAt(grid_ix, grid_iy, value_a) || + !grid->valueAt(grid_ix2, grid_iy, value_b) || + !grid->valueAt(grid_ix, grid_iy2, value_c) || + !grid->valueAt(grid_ix2, grid_iy2, value_d)); + if (grid->hasChanged()) { + if (curGridset->reopen(ctx)) { + return read_vgrid_value(ctx, grids, input, vmultiplier); + } + error = true; + } + + if (error) { + return HUGE_VAL; + } + + double total_weight = 0.0; + int n_weights = 0; + double value = 0.0f; + + if (!grid->isNodata(value_a, vmultiplier)) { + double weight = (1.0 - grid_x) * (1.0 - grid_y); + value += value_a * weight; + total_weight += weight; + n_weights++; + } + if (!grid->isNodata(value_b, vmultiplier)) { + double weight = (grid_x) * (1.0 - grid_y); + value += value_b * weight; + total_weight += weight; + n_weights++; + } + if (!grid->isNodata(value_c, vmultiplier)) { + double weight = (1.0 - grid_x) * (grid_y); + value += value_c * weight; + total_weight += weight; + n_weights++; + } + if (!grid->isNodata(value_d, vmultiplier)) { + double weight = (grid_x) * (grid_y); + value += value_d * weight; + total_weight += weight; + n_weights++; + } + if (n_weights == 0) + value = HUGE_VAL; + else if (n_weights != 4) + value /= total_weight; + + return value * vmultiplier; +} + +/**********************************************/ +ListOfVGrids pj_vgrid_init(PJ *P, const char *gridkey) { + /********************************************** + + Initizalize and populate gridlist. + + Takes a PJ-object and the plus-parameter + name that is used in the proj-string to + specify the grids to load, e.g. "+grids". + The + should be left out here. + + Returns the number of loaded grids. + + ***********************************************/ + + std::string key("s"); + key += gridkey; + const char *gridnames = pj_param(P->ctx, P->params, key.c_str()).s; + if (gridnames == nullptr) + return {}; + + auto listOfGridNames = internal::split(std::string(gridnames), ','); + ListOfVGrids grids; + for (const auto &gridnameStr : listOfGridNames) { + const char *gridname = gridnameStr.c_str(); + bool canFail = false; + if (gridname[0] == '@') { + canFail = true; + gridname++; + } + auto gridSet = VerticalShiftGridSet::open(P->ctx, gridname); + if (!gridSet) { + if (!canFail) { + pj_ctx_set_errno(P->ctx, PJD_ERR_FAILED_TO_LOAD_GRID); + return {}; + } + } else { + grids.emplace_back(std::move(gridSet)); + } + } + + return grids; +} + +/***********************************************/ +double pj_vgrid_value(PJ *P, const ListOfVGrids &grids, PJ_LP lp, + double vmultiplier) { + /*********************************************** + + Read grid value at position lp in grids loaded + with proj_grid_init. + + Returns the grid value of the given coordinate. + + ************************************************/ + + double value; + + value = read_vgrid_value(P->ctx, grids, lp, vmultiplier); + proj_log_trace(P, "proj_vgrid_value: (%f, %f) = %f", lp.lam * RAD_TO_DEG, + lp.phi * RAD_TO_DEG, value); + + return value; +} + +// --------------------------------------------------------------------------- + +const GenericShiftGrid *pj_find_generic_grid(const ListOfGenericGrids &grids, + const PJ_LP &input, + GenericShiftGridSet *&gridSetOut) { + for (const auto &gridset : grids) { + auto grid = gridset->gridAt(input.lam, input.phi); + if (grid) { + gridSetOut = gridset.get(); + return grid; + } + } + return nullptr; +} + +// --------------------------------------------------------------------------- + +bool pj_bilinear_interpolation_three_samples(const GenericShiftGrid *grid, + const PJ_LP &lp, int idx1, + int idx2, int idx3, double &v1, + double &v2, double &v3, + bool &must_retry) { + must_retry = false; + + const auto &extent = grid->extentAndRes(); + double grid_x = (lp.lam - extent.westLon) / extent.resLon; + double grid_y = (lp.phi - extent.southLat) / extent.resLat; + int ix = static_cast<int>(grid_x); + int iy = static_cast<int>(grid_y); + int ix2 = std::min(ix + 1, grid->width() - 1); + int iy2 = std::min(iy + 1, grid->height() - 1); + + float dx1 = 0.0f, dy1 = 0.0f, dz1 = 0.0f; + float dx2 = 0.0f, dy2 = 0.0f, dz2 = 0.0f; + float dx3 = 0.0f, dy3 = 0.0f, dz3 = 0.0f; + float dx4 = 0.0f, dy4 = 0.0f, dz4 = 0.0f; + bool error = (!grid->valueAt(ix, iy, idx1, dx1) || + !grid->valueAt(ix, iy, idx2, dy1) || + !grid->valueAt(ix, iy, idx3, dz1) || + !grid->valueAt(ix2, iy, idx1, dx2) || + !grid->valueAt(ix2, iy, idx2, dy2) || + !grid->valueAt(ix2, iy, idx3, dz2) || + !grid->valueAt(ix, iy2, idx1, dx3) || + !grid->valueAt(ix, iy2, idx2, dy3) || + !grid->valueAt(ix, iy2, idx3, dz3) || + !grid->valueAt(ix2, iy2, idx1, dx4) || + !grid->valueAt(ix2, iy2, idx2, dy4) || + !grid->valueAt(ix2, iy2, idx3, dz4)); + if (grid->hasChanged()) { + must_retry = true; + return false; + } + if (error) { + return false; + } + + double frct_lam = grid_x - ix; + double frct_phi = grid_y - iy; + 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; + + v1 = m00 * dx1 + m10 * dx2 + m01 * dx3 + m11 * dx4; + v2 = m00 * dy1 + m10 * dy2 + m01 * dy3 + m11 * dy4; + v3 = m00 * dz1 + m10 * dz2 + m01 * dz3 + m11 * dz4; + return true; +} + NS_PROJ_END + +/************************************************************************/ +/* 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 = NS_PROJ::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 = + pj_hgrid_apply(ctx, hgrids, input, inverse ? PJ_INV : PJ_FWD); + + 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/grids.hpp b/src/grids.hpp index aa852ef6..6b6ee0d2 100644 --- a/src/grids.hpp +++ b/src/grids.hpp @@ -70,6 +70,7 @@ class Grid { const std::string &name() const { return m_name; } virtual bool isNullGrid() const { return false; } + virtual bool hasChanged() const = 0; }; // --------------------------------------------------------------------------- @@ -116,6 +117,7 @@ class VerticalShiftGridSet { const VerticalShiftGrid *gridAt(double lon, double lat) const; virtual void reassign_context(PJ_CONTEXT *ctx); + virtual bool reopen(PJ_CONTEXT *ctx); }; // --------------------------------------------------------------------------- @@ -162,6 +164,7 @@ class HorizontalShiftGridSet { const HorizontalShiftGrid *gridAt(double lon, double lat) const; virtual void reassign_context(PJ_CONTEXT *ctx); + virtual bool reopen(PJ_CONTEXT *ctx); }; // --------------------------------------------------------------------------- @@ -217,6 +220,7 @@ class GenericShiftGridSet { const GenericShiftGrid *gridAt(double lon, double lat) const; virtual void reassign_context(PJ_CONTEXT *ctx); + virtual bool reopen(PJ_CONTEXT *ctx); }; // --------------------------------------------------------------------------- @@ -225,15 +229,24 @@ typedef std::vector<std::unique_ptr<HorizontalShiftGridSet>> ListOfHGrids; typedef std::vector<std::unique_ptr<VerticalShiftGridSet>> ListOfVGrids; typedef std::vector<std::unique_ptr<GenericShiftGridSet>> ListOfGenericGrids; -ListOfVGrids proj_vgrid_init(PJ *P, const char *grids); -ListOfHGrids proj_hgrid_init(PJ *P, const char *grids); -ListOfGenericGrids proj_generic_grid_init(PJ *P, const char *grids); - -double proj_vgrid_value(PJ *P, const ListOfVGrids &, PJ_LP lp, - double vmultiplier); -PJ_LP proj_hgrid_value(PJ *P, const ListOfHGrids &, PJ_LP lp); -PJ_LP proj_hgrid_apply(PJ *P, const ListOfHGrids &, PJ_LP lp, - PJ_DIRECTION direction); +ListOfVGrids pj_vgrid_init(PJ *P, const char *grids); +ListOfHGrids pj_hgrid_init(PJ *P, const char *grids); +ListOfGenericGrids pj_generic_grid_init(PJ *P, const char *grids); + +PJ_LP pj_hgrid_value(PJ *P, const ListOfHGrids &grids, PJ_LP lp); +double pj_vgrid_value(PJ *P, const ListOfVGrids &, PJ_LP lp, + double vmultiplier); +PJ_LP pj_hgrid_apply(PJ_CONTEXT *ctx, const ListOfHGrids &grids, PJ_LP lp, + PJ_DIRECTION direction); + +const GenericShiftGrid *pj_find_generic_grid(const ListOfGenericGrids &grids, + const PJ_LP &input, + GenericShiftGridSet *&gridSetOut); +bool pj_bilinear_interpolation_three_samples(const GenericShiftGrid *grid, + const PJ_LP &lp, int idx1, + int idx2, int idx3, double &v1, + double &v2, double &v3, + bool &must_retry); NS_PROJ_END diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake index 18c91021..fdb59434 100644 --- a/src/lib_proj.cmake +++ b/src/lib_proj.cmake @@ -220,8 +220,6 @@ set(SRC_LIBPROJ_CORE 4D_api.cpp aasincos.cpp adjlon.cpp - apply_gridshift.cpp - apply_vgridshift.cpp auth.cpp ctx.cpp datum_set.cpp diff --git a/src/open_lib.cpp b/src/open_lib.cpp index 23ee79a0..ae387281 100644 --- a/src/open_lib.cpp +++ b/src/open_lib.cpp @@ -224,6 +224,16 @@ static bool is_rel_or_absolute_filename(const char *name) || (name[0] != '\0' && name[1] == ':' && strchr(dir_chars,name[2])); } +static bool ignoreUserWritableDirectory() +{ + // Env var mostly for testing purposes and being independent from + // an existing installation + const char* envVarIgnoreUserWritableDirectory = + getenv("PROJ_IGNORE_USER_WRITABLE_DIRECTORY"); + return envVarIgnoreUserWritableDirectory != nullptr && + envVarIgnoreUserWritableDirectory[0] != '\0'; +} + static void* pj_open_lib_internal(projCtx ctx, const char *name, const char *mode, void* (*open_file)(projCtx, const char*, const char*), @@ -279,6 +289,17 @@ pj_open_lib_internal(projCtx ctx, const char *name, const char *mode, break; } } + + else if( !ignoreUserWritableDirectory() && + (fid = open_file(ctx, + (pj_context_get_user_writable_directory(ctx, false) + + DIR_CHAR + name).c_str(), mode)) != nullptr ) { + fname = pj_context_get_user_writable_directory(ctx, false); + fname += DIR_CHAR; + fname += name; + sysname = fname.c_str(); + } + /* if is environment PROJ_LIB defined */ else if ((sysname = getenv("PROJ_LIB")) != nullptr) { auto paths = NS_PROJ::internal::split(std::string(sysname), dirSeparator); @@ -381,6 +402,23 @@ std::unique_ptr<NS_PROJ::File> NS_PROJ::FileManager::open_resource_file( "Using %s", remote_file.c_str() ); pj_ctx_set_errno( ctx, 0 ); } + } else { + // For example for resource files like 'alaska' + auto remote_file_tif = remote_file + ".tif"; + file = open(ctx, remote_file_tif.c_str()); + if( file ) { + pj_log( ctx, PJ_LOG_DEBUG_MAJOR, + "Using %s", remote_file_tif.c_str() ); + pj_ctx_set_errno( ctx, 0 ); + } else { + // Init files + file = open(ctx, remote_file.c_str()); + if( file ) { + pj_log( ctx, PJ_LOG_DEBUG_MAJOR, + "Using %s", remote_file.c_str() ); + pj_ctx_set_errno( ctx, 0 ); + } + } } } } @@ -401,6 +401,9 @@ typedef const char* (*proj_network_get_header_value_cbk_type)( * * Read size_to_read bytes from handle, starting at offset, into * buffer. + * During this read, the implementation should make sure to store the HTTP + * headers from the server response to be able to respond to + * proj_network_get_header_value_cbk_type callback. * * error_string_max_size should be the maximum size that can be written into * the out_error_string buffer (including terminating nul character). @@ -440,6 +443,15 @@ void PROJ_DLL proj_grid_cache_set_ttl(PJ_CONTEXT* ctx, int ttl_seconds); void PROJ_DLL proj_grid_cache_clear(PJ_CONTEXT* ctx); +int PROJ_DLL proj_is_download_needed(PJ_CONTEXT* ctx, + const char* url_or_filename, + int ignore_ttl_setting); +int PROJ_DLL proj_download_file(PJ_CONTEXT *ctx, const char *url_or_filename, + int ignore_ttl_setting, + int (*progress_cbk)(PJ_CONTEXT *, double pct, + void *user_data), + void *user_data); + /*! @cond Doxygen_Suppress */ /* Manage the transformation definition object PJ */ diff --git a/src/proj_internal.h b/src/proj_internal.h index 63c53551..ce7b9d74 100644 --- a/src/proj_internal.h +++ b/src/proj_internal.h @@ -709,6 +709,7 @@ struct projCtx_t { bool iniFileLoaded = false; std::string endpoint{}; + std::string user_writable_directory{}; projGridChunkCache gridChunkCache{}; int projStringParserCreateFromPROJStringRecursionCounter = 0; // to avoid potential infinite recursion in PROJStringParser::createFromPROJString() @@ -844,8 +845,9 @@ std::string pj_context_get_url_endpoint(PJ_CONTEXT* ctx); void pj_load_ini(PJ_CONTEXT* ctx); -// For testing purposes +// Exported for testing purposes only std::string PROJ_DLL pj_context_get_grid_cache_filename(PJ_CONTEXT *ctx); +std::string PROJ_DLL pj_context_get_user_writable_directory(PJ_CONTEXT *ctx, bool create); /* classic public API */ #include "proj_api.h" diff --git a/src/sqlite3.cpp b/src/sqlite3.cpp index 0c89c0b9..90e73c2a 100644 --- a/src/sqlite3.cpp +++ b/src/sqlite3.cpp @@ -25,8 +25,17 @@ * DEALINGS IN THE SOFTWARE. *****************************************************************************/ +#ifdef __GNUC__ +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Weffc++" +#endif + #include "sqlite3.hpp" +#ifdef __GNUC__ +#pragma GCC diagnostic pop +#endif + #include <cstdlib> #include <cstring> #include <sstream> // std::ostringstream @@ -182,4 +191,4 @@ SQLiteStatement::SQLiteStatement(sqlite3_stmt *hStmtIn) : hStmt(hStmtIn) {} // --------------------------------------------------------------------------- -NS_PROJ_END
\ No newline at end of file +NS_PROJ_END diff --git a/src/transform.cpp b/src/transform.cpp index 020d62ea..811e2a6a 100644 --- a/src/transform.cpp +++ b/src/transform.cpp @@ -451,7 +451,7 @@ static int pj_apply_vgridshift( PJ *defn, if( defn->vgrids_legacy == nullptr ) { defn->vgrids_legacy = new ListOfVGrids; - auto vgrids = proj_vgrid_init(defn, "geoidgrids"); + auto vgrids = pj_vgrid_init(defn, "geoidgrids"); if( vgrids.empty() ) return 0; *static_cast<ListOfVGrids*>(defn->vgrids_legacy) = std::move(vgrids); @@ -470,7 +470,7 @@ static int pj_apply_vgridshift( PJ *defn, input.phi = y[io]; input.lam = x[io]; - value = proj_vgrid_value(defn, *static_cast<ListOfVGrids*>(defn->vgrids_legacy), input, 1.0); + value = pj_vgrid_value(defn, *static_cast<ListOfVGrids*>(defn->vgrids_legacy), input, 1.0); if( inverse ) z[io] -= value; @@ -896,7 +896,7 @@ int pj_apply_gridshift_2( PJ *defn, int inverse, if( defn->hgrids_legacy == nullptr ) { defn->hgrids_legacy = new ListOfHGrids; - auto hgrids = proj_hgrid_init(defn, "nadgrids"); + auto hgrids = pj_hgrid_init(defn, "nadgrids"); if( hgrids.empty() ) return 0; *static_cast<ListOfHGrids*>(defn->hgrids_legacy) = std::move(hgrids); @@ -914,7 +914,7 @@ int pj_apply_gridshift_2( PJ *defn, int inverse, input.phi = y[io]; input.lam = x[io]; - auto output = proj_hgrid_apply(defn, *static_cast<ListOfHGrids*>(defn->hgrids_legacy), input, inverse ? PJ_INV : PJ_FWD); + auto output = pj_hgrid_apply(defn->ctx, *static_cast<ListOfHGrids*>(defn->hgrids_legacy), input, inverse ? PJ_INV : PJ_FWD); if ( output.lam != HUGE_VAL ) { diff --git a/src/transformations/deformation.cpp b/src/transformations/deformation.cpp index eb109826..8aee50c9 100644 --- a/src/transformations/deformation.cpp +++ b/src/transformations/deformation.cpp @@ -80,20 +80,6 @@ struct deformationData { // --------------------------------------------------------------------------- -static const GenericShiftGrid* findGrid(const ListOfGenericGrids& grids, - const PJ_LP& input) -{ - for( const auto& gridset: grids ) - { - auto grid = gridset->gridAt(input.lam, input.phi); - if( grid ) - return grid; - } - return nullptr; -} - -// --------------------------------------------------------------------------- - static bool get_grid_values(PJ* P, deformationData* Q, const PJ_LP& lp, @@ -101,7 +87,8 @@ static bool get_grid_values(PJ* P, double& vy, double& vz) { - auto grid = findGrid(Q->grids, lp); + GenericShiftGridSet* gridset = nullptr; + auto grid = pj_find_generic_grid(Q->grids, lp, gridset); if( !grid ) { return false; } @@ -136,57 +123,20 @@ static bool get_grid_values(PJ* P, return false; } - const auto& extent = grid->extentAndRes(); - double grid_x = (lp.lam - extent.westLon) / extent.resLon; - double grid_y = (lp.phi - extent.southLat) / extent.resLat; - int ix = static_cast<int>(grid_x); - int iy = static_cast<int>(grid_y); - int ix2 = std::min(ix + 1, grid->width() - 1); - int iy2 = std::min(iy + 1, grid->height() - 1); - - float dx1, dy1, dz1; - if( !grid->valueAt(ix, iy, sampleE, dx1) || - !grid->valueAt(ix, iy, sampleN, dy1) || - !grid->valueAt(ix, iy, sampleU, dz1) ) { - return false; - } - - float dx2, dy2, dz2; - if( !grid->valueAt(ix2, iy, sampleE, dx2) || - !grid->valueAt(ix2, iy, sampleN, dy2) || - !grid->valueAt(ix2, iy, sampleU, dz2) ) { - return false; - } - - float dx3, dy3, dz3; - if( !grid->valueAt(ix, iy2, sampleE, dx3) || - !grid->valueAt(ix, iy2, sampleN, dy3) || - !grid->valueAt(ix, iy2, sampleU, dz3) ) { - return false; - } - - float dx4, dy4, dz4; - if( !grid->valueAt(ix2, iy2, sampleE, dx4) || - !grid->valueAt(ix2, iy2, sampleN, dy4) || - !grid->valueAt(ix2, iy2, sampleU, dz4) ) { + bool must_retry = false; + if( !pj_bilinear_interpolation_three_samples(grid, lp, + sampleE, sampleN, sampleU, + vx, vy, vz, + must_retry) ) + { + if( must_retry ) + return get_grid_values( P, Q, lp, vx, vy, vz); return false; } - - double frct_lam = grid_x - ix; - double frct_phi = grid_y - iy; - 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; // divide by 1000 to get m/year - vx = (m00 * dx1 + m10 * dx2 + m01 * dx3 + m11 * dx4) / 1000; - vy = (m00 * dy1 + m10 * dy2 + m01 * dy3 + m11 * dy4) / 1000; - vz = (m00 * dz1 + m10 * dz2 + m01 * dz3 + m11 * dz4) / 1000; + vx /= 1000; + vy /= 1000; + vz /= 1000; return true; } @@ -226,8 +176,8 @@ static PJ_XYZ get_grid_shift(PJ* P, const PJ_XYZ& cartesian) { } else { - shift.lp = proj_hgrid_value(P, Q->hgrids, geodetic.lp); - shift.enu.u = proj_vgrid_value(P, Q->vgrids, geodetic.lp, 1.0); + shift.lp = pj_hgrid_value(P, Q->hgrids, geodetic.lp); + shift.enu.u = pj_vgrid_value(P, Q->vgrids, geodetic.lp, 1.0); if (proj_errno(P) == PJD_ERR_GRID_AREA) proj_log_debug(P, "deformation: coordinate (%.3f, %.3f) outside deformation model", @@ -425,7 +375,7 @@ PJ *TRANSFORMATION(deformation,1) { if( has_grids ) { - Q->grids = proj_generic_grid_init(P, "grids"); + Q->grids = pj_generic_grid_init(P, "grids"); /* Was gridlist compiled properly? */ if ( proj_errno(P) ) { proj_log_error(P, "deformation: could not find required grid(s)."); @@ -434,13 +384,13 @@ PJ *TRANSFORMATION(deformation,1) { } else { - Q->hgrids = proj_hgrid_init(P, "xy_grids"); + Q->hgrids = pj_hgrid_init(P, "xy_grids"); if (proj_errno(P)) { proj_log_error(P, "deformation: could not find requested xy_grid(s)."); return destructor(P, PJD_ERR_FAILED_TO_LOAD_GRID); } - Q->vgrids = proj_vgrid_init(P, "z_grids"); + Q->vgrids = pj_vgrid_init(P, "z_grids"); if (proj_errno(P)) { proj_log_error(P, "deformation: could not find requested z_grid(s)."); return destructor(P, PJD_ERR_FAILED_TO_LOAD_GRID); diff --git a/src/transformations/hgridshift.cpp b/src/transformations/hgridshift.cpp index 24da4dde..122a7ab2 100644 --- a/src/transformations/hgridshift.cpp +++ b/src/transformations/hgridshift.cpp @@ -28,7 +28,7 @@ static PJ_XYZ forward_3d(PJ_LPZ lpz, PJ *P) { if ( Q->defer_grid_opening ) { Q->defer_grid_opening = false; - Q->grids = proj_hgrid_init(P, "grids"); + Q->grids = pj_hgrid_init(P, "grids"); if ( proj_errno(P) ) { return proj_coord_error().xyz; } @@ -37,7 +37,7 @@ static PJ_XYZ forward_3d(PJ_LPZ lpz, PJ *P) { if (!Q->grids.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, Q->grids, point.lp, PJ_FWD); + point.lp = pj_hgrid_apply(P->ctx, Q->grids, point.lp, PJ_FWD); } return point.xyz; @@ -51,7 +51,7 @@ static PJ_LPZ reverse_3d(PJ_XYZ xyz, PJ *P) { if ( Q->defer_grid_opening ) { Q->defer_grid_opening = false; - Q->grids = proj_hgrid_init(P, "grids"); + Q->grids = pj_hgrid_init(P, "grids"); if ( proj_errno(P) ) { return proj_coord_error().lpz; } @@ -60,7 +60,7 @@ static PJ_LPZ reverse_3d(PJ_XYZ xyz, PJ *P) { if (!Q->grids.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, Q->grids, point.lp, PJ_INV); + point.lp = pj_hgrid_apply(P->ctx, Q->grids, point.lp, PJ_INV); } return point.lpz; @@ -165,7 +165,7 @@ PJ *TRANSFORMATION(hgridshift,0) { Q->defer_grid_opening = true; } else { - Q->grids = proj_hgrid_init(P, "grids"); + Q->grids = pj_hgrid_init(P, "grids"); /* Was gridlist compiled properly? */ if ( proj_errno(P) ) { proj_log_error(P, "hgridshift: could not find required grid(s)."); diff --git a/src/transformations/vgridshift.cpp b/src/transformations/vgridshift.cpp index 3e7a015e..121b795a 100644 --- a/src/transformations/vgridshift.cpp +++ b/src/transformations/vgridshift.cpp @@ -56,7 +56,7 @@ static PJ_XYZ forward_3d(PJ_LPZ lpz, PJ *P) { if ( Q->defer_grid_opening ) { Q->defer_grid_opening = false; - Q->grids = proj_vgrid_init(P, "grids"); + Q->grids = pj_vgrid_init(P, "grids"); deal_with_vertcon_gtx_hack(P); if ( proj_errno(P) ) { return proj_coord_error().xyz; @@ -66,7 +66,7 @@ static PJ_XYZ forward_3d(PJ_LPZ lpz, PJ *P) { if (!Q->grids.empty()) { /* Only try the gridshift if at least one grid is loaded, * otherwise just pass the coordinate through unchanged. */ - point.xyz.z += proj_vgrid_value(P, Q->grids, point.lp, Q->forward_multiplier); + point.xyz.z += pj_vgrid_value(P, Q->grids, point.lp, Q->forward_multiplier); } return point.xyz; @@ -80,7 +80,7 @@ static PJ_LPZ reverse_3d(PJ_XYZ xyz, PJ *P) { if ( Q->defer_grid_opening ) { Q->defer_grid_opening = false; - Q->grids = proj_vgrid_init(P, "grids"); + Q->grids = pj_vgrid_init(P, "grids"); deal_with_vertcon_gtx_hack(P); if ( proj_errno(P) ) { return proj_coord_error().lpz; @@ -90,7 +90,7 @@ static PJ_LPZ reverse_3d(PJ_XYZ xyz, PJ *P) { if (!Q->grids.empty()) { /* Only try the gridshift if at least one grid is loaded, * otherwise just pass the coordinate through unchanged. */ - point.xyz.z -= proj_vgrid_value(P, Q->grids, point.lp, Q->forward_multiplier); + point.xyz.z -= pj_vgrid_value(P, Q->grids, point.lp, Q->forward_multiplier); } return point.lpz; @@ -193,7 +193,7 @@ PJ *TRANSFORMATION(vgridshift,0) { } else { /* Build gridlist. P->vgridlist_geoid can be empty if +grids only ask for optional grids. */ - Q->grids = proj_vgrid_init(P, "grids"); + Q->grids = pj_vgrid_init(P, "grids"); /* Was gridlist compiled properly? */ if ( proj_errno(P) ) { diff --git a/src/transformations/xyzgridshift.cpp b/src/transformations/xyzgridshift.cpp index a76f3255..e1c76518 100644 --- a/src/transformations/xyzgridshift.cpp +++ b/src/transformations/xyzgridshift.cpp @@ -51,21 +51,6 @@ struct xyzgridshiftData { }; } // anonymous namespace - -// --------------------------------------------------------------------------- - -static const GenericShiftGrid* findGrid(const ListOfGenericGrids& grids, - const PJ_LP& input) -{ - for( const auto& gridset: grids ) - { - auto grid = gridset->gridAt(input.lam, input.phi); - if( grid ) - return grid; - } - return nullptr; -} - // --------------------------------------------------------------------------- static bool get_grid_values(PJ* P, @@ -77,13 +62,14 @@ static bool get_grid_values(PJ* P, { if ( Q->defer_grid_opening ) { Q->defer_grid_opening = false; - Q->grids = proj_generic_grid_init(P, "grids"); + Q->grids = pj_generic_grid_init(P, "grids"); if ( proj_errno(P) ) { return false; } } - auto grid = findGrid(Q->grids, lp); + GenericShiftGridSet* gridset = nullptr; + auto grid = pj_find_generic_grid(Q->grids, lp, gridset); if( !grid ) { return false; } @@ -118,56 +104,20 @@ static bool get_grid_values(PJ* P, return false; } - const auto& extent = grid->extentAndRes(); - double grid_x = (lp.lam - extent.westLon) / extent.resLon; - double grid_y = (lp.phi - extent.southLat) / extent.resLat; - int ix = static_cast<int>(grid_x); - int iy = static_cast<int>(grid_y); - int ix2 = std::min(ix + 1, grid->width() - 1); - int iy2 = std::min(iy + 1, grid->height() - 1); - - float dx1, dy1, dz1; - if( !grid->valueAt(ix, iy, sampleX, dx1) || - !grid->valueAt(ix, iy, sampleY, dy1) || - !grid->valueAt(ix, iy, sampleZ, dz1) ) { - return false; - } - - float dx2, dy2, dz2; - if( !grid->valueAt(ix2, iy, sampleX, dx2) || - !grid->valueAt(ix2, iy, sampleY, dy2) || - !grid->valueAt(ix2, iy, sampleZ, dz2) ) { - return false; - } - - float dx3, dy3, dz3; - if( !grid->valueAt(ix, iy2, sampleX, dx3) || - !grid->valueAt(ix, iy2, sampleY, dy3) || - !grid->valueAt(ix, iy2, sampleZ, dz3) ) { - return false; - } - - float dx4, dy4, dz4; - if( !grid->valueAt(ix2, iy2, sampleX, dx4) || - !grid->valueAt(ix2, iy2, sampleY, dy4) || - !grid->valueAt(ix2, iy2, sampleZ, dz4) ) { + bool must_retry = false; + if( !pj_bilinear_interpolation_three_samples(grid, lp, + sampleX, sampleY, sampleZ, + dx, dy, dz, + must_retry) ) + { + if( must_retry ) + return get_grid_values( P, Q, lp, dx, dy, dz); return false; } - double frct_lam = grid_x - ix; - double frct_phi = grid_y - iy; - 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; - dx = (m00 * dx1 + m10 * dx2 + m01 * dx3 + m11 * dx4) * Q->multiplier; - dy = (m00 * dy1 + m10 * dy2 + m01 * dy3 + m11 * dy4) * Q->multiplier; - dz = (m00 * dz1 + m10 * dz2 + m01 * dz3 + m11 * dz4) * Q->multiplier; + dx *= Q->multiplier; + dy *= Q->multiplier; + dz *= Q->multiplier; return true; } @@ -341,7 +291,7 @@ PJ *TRANSFORMATION(xyzgridshift,0) { Q->defer_grid_opening = true; } else { - Q->grids = proj_generic_grid_init(P, "grids"); + Q->grids = pj_generic_grid_init(P, "grids"); /* Was gridlist compiled properly? */ if ( proj_errno(P) ) { proj_log_error(P, "xyzgridshift: could not find required grid(s)."); |