From 126384854bf8e1b7461aebcc28966a6559971de1 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 4 Dec 2019 18:56:15 +0100 Subject: vertical grid shift: rework to no longer load whole grid into memory --- src/transform.cpp | 69 +++++++++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 65 insertions(+), 4 deletions(-) (limited to 'src/transform.cpp') diff --git a/src/transform.cpp b/src/transform.cpp index d111d835..863da605 100644 --- a/src/transform.cpp +++ b/src/transform.cpp @@ -431,6 +431,70 @@ static int height_unit (PJ *P, PJ_DIRECTION dir, long n, int dist, double *z) { } +/************************************************************************/ +/* pj_apply_vgridshift() */ +/* */ +/* This implementation takes uses the gridlist from a coordinate */ +/* system definition. If the gridlist has not yet been */ +/* populated in the coordinate system definition we set it up */ +/* now. */ +/************************************************************************/ +static int pj_apply_vgridshift( PJ *defn, + int inverse, + long point_count, int point_offset, + double *x, double *y, double *z ) + +{ + if( defn->vgrids.empty() ) + { + proj_vgrid_init(defn, "geoidgrids"); + } + + for( int i = 0; i < point_count; i++ ) + { + double value; + long io = i * point_offset; + PJ_LP input; + + input.phi = y[io]; + input.lam = x[io]; + + value = proj_vgrid_value(defn, input, 1.0); + + if( inverse ) + z[io] -= value; + else + z[io] += value; + + if( value == HUGE_VAL ) + { + std::string gridlist; + + proj_log_debug(defn, + "pj_apply_vgridshift(): failed to find a grid shift table for\n" + " location (%.7fdW,%.7fdN)", + x[io] * RAD_TO_DEG, + y[io] * RAD_TO_DEG ); + + for( const auto& gridset: defn->vgrids ) + { + if( gridlist.empty() ) + gridlist += " tried: "; + else + gridlist += ','; + gridlist += gridset->name(); + } + + proj_log_debug(defn, "%s", gridlist.c_str()); + pj_ctx_set_errno( defn->ctx, PJD_ERR_GRID_AREA ); + + return PJD_ERR_GRID_AREA; + } + } + + return 0; +} + /* -------------------------------------------------------------------- */ /* Transform to ellipsoidal heights if needed */ @@ -441,10 +505,7 @@ static int geometric_to_orthometric (PJ *P, PJ_DIRECTION dir, long n, int dist, return 0; if (z==nullptr) return PJD_ERR_GEOCENTRIC; - err = pj_apply_vgridshift (P, "sgeoidgrids", - &(P->vgridlist_geoid), - &(P->vgridlist_geoid_count), - dir==PJ_FWD ? 1 : 0, n, dist, x, y, z ); + err = pj_apply_vgridshift (P, dir==PJ_FWD ? 1 : 0, n, dist, x, y, z ); if (err) return pj_ctx_get_errno(P->ctx); return 0; -- cgit v1.2.3 From 6875ef7116b9dab4021afeb06e2b79cd5679743b Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Thu, 5 Dec 2019 01:03:46 +0100 Subject: horizontal grid shift: rework to no longer load whole grid into memory --- src/transform.cpp | 57 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 57 insertions(+) (limited to 'src/transform.cpp') diff --git a/src/transform.cpp b/src/transform.cpp index 863da605..f6002b70 100644 --- a/src/transform.cpp +++ b/src/transform.cpp @@ -867,6 +867,58 @@ int pj_geocentric_from_wgs84( PJ *defn, return 0; } + +/************************************************************************/ +/* pj_apply_gridshift_2() */ +/* */ +/* This implementation uses the gridlist from a coordinate */ +/* system definition. If the gridlist has not yet been */ +/* populated in the coordinate system definition we set it up */ +/* now. */ +/************************************************************************/ +static +int pj_apply_gridshift_2( PJ *defn, int inverse, + long point_count, int point_offset, + double *x, double *y, double * /*z*/ ) + +{ + if( defn->hgrids.empty() ) + { + proj_hgrid_init(defn, "nadgrids"); + } + + for( long i = 0; i < point_count; i++ ) + { + PJ_LP input; + + long io = i * point_offset; + input.phi = y[io]; + input.lam = x[io]; + + auto output = proj_hgrid_apply(defn, input, inverse ? PJ_INV : PJ_FWD); + + if ( output.lam != HUGE_VAL ) + { + y[io] = output.phi; + x[io] = output.lam; + } + else + { + if( defn->ctx->debug_level >= PJ_LOG_DEBUG_MAJOR ) + { + pj_log( defn->ctx, PJ_LOG_DEBUG_MAJOR, + "pj_apply_gridshift(): failed to find a grid shift table for\n" + " location (%.7fdW,%.7fdN)", + x[io] * RAD_TO_DEG, + y[io] * RAD_TO_DEG ); + } + } + } + + return 0; +} + + /************************************************************************/ /* pj_datum_transform() */ /* */ @@ -1107,3 +1159,8 @@ static int adjust_axis( projCtx ctx, return 0; } +// --------------------------------------------------------------------------- + +void pj_deallocate_grids() +{ +} -- cgit v1.2.3 From 916a92062ffa2f2b59007047fae2176bbb463ca3 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Fri, 6 Dec 2019 14:07:27 +0100 Subject: Remove hgrids and vgrids member from PJ structure, and store them in hgridshift/vgridshift/deformation structures --- src/transform.cpp | 33 ++++++++++++++++++++++++++------- 1 file changed, 26 insertions(+), 7 deletions(-) (limited to 'src/transform.cpp') diff --git a/src/transform.cpp b/src/transform.cpp index f6002b70..020d62ea 100644 --- a/src/transform.cpp +++ b/src/transform.cpp @@ -34,6 +34,9 @@ #include "proj.h" #include "proj_internal.h" #include "geocent.h" +#include "grids.hpp" + +using namespace NS_PROJ; static int adjust_axis( projCtx ctx, const char *axis, int denormalize_flag, long point_count, int point_offset, @@ -445,9 +448,17 @@ static int pj_apply_vgridshift( PJ *defn, double *x, double *y, double *z ) { - if( defn->vgrids.empty() ) + if( defn->vgrids_legacy == nullptr ) + { + defn->vgrids_legacy = new ListOfVGrids; + auto vgrids = proj_vgrid_init(defn, "geoidgrids"); + if( vgrids.empty() ) + return 0; + *static_cast(defn->vgrids_legacy) = std::move(vgrids); + } + if( static_cast(defn->vgrids_legacy)->empty() ) { - proj_vgrid_init(defn, "geoidgrids"); + return 0; } for( int i = 0; i < point_count; i++ ) @@ -459,7 +470,7 @@ static int pj_apply_vgridshift( PJ *defn, input.phi = y[io]; input.lam = x[io]; - value = proj_vgrid_value(defn, input, 1.0); + value = proj_vgrid_value(defn, *static_cast(defn->vgrids_legacy), input, 1.0); if( inverse ) z[io] -= value; @@ -476,7 +487,7 @@ static int pj_apply_vgridshift( PJ *defn, x[io] * RAD_TO_DEG, y[io] * RAD_TO_DEG ); - for( const auto& gridset: defn->vgrids ) + for( const auto& gridset: *static_cast(defn->vgrids_legacy) ) { if( gridlist.empty() ) gridlist += " tried: "; @@ -882,9 +893,17 @@ int pj_apply_gridshift_2( PJ *defn, int inverse, double *x, double *y, double * /*z*/ ) { - if( defn->hgrids.empty() ) + if( defn->hgrids_legacy == nullptr ) { - proj_hgrid_init(defn, "nadgrids"); + defn->hgrids_legacy = new ListOfHGrids; + auto hgrids = proj_hgrid_init(defn, "nadgrids"); + if( hgrids.empty() ) + return 0; + *static_cast(defn->hgrids_legacy) = std::move(hgrids); + } + if( static_cast(defn->hgrids_legacy)->empty() ) + { + return 0; } for( long i = 0; i < point_count; i++ ) @@ -895,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, input, inverse ? PJ_INV : PJ_FWD); + auto output = proj_hgrid_apply(defn, *static_cast(defn->hgrids_legacy), input, inverse ? PJ_INV : PJ_FWD); if ( output.lam != HUGE_VAL ) { -- cgit v1.2.3