From eafeb61ce59aeb34dabf38f55f70ba9a3b779c4b Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Mon, 6 Jan 2020 23:03:29 +0100 Subject: Grid refactoring: address review comments of https://github.com/OSGeo/PROJ/pull/1777 - Move content of legacy apply_gridshift.cpp and apply_vgridshift.cpp in grids.cpp - Rename nad_ functions to pj_hgrid_ - Rename internal proj_hgrid_/proj_vgrid_ functions to pj_ --- src/transformations/deformation.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) (limited to 'src/transformations/deformation.cpp') diff --git a/src/transformations/deformation.cpp b/src/transformations/deformation.cpp index eb109826..993647fc 100644 --- a/src/transformations/deformation.cpp +++ b/src/transformations/deformation.cpp @@ -226,8 +226,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 +425,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 +434,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); -- cgit v1.2.3 From 237296b7e84a8bb270e3be06a690737a601d73e7 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 7 Jan 2020 03:37:23 +0100 Subject: Remote grid: add mechanism to re-open a grid if it has changed while being opened --- src/transformations/deformation.cpp | 48 +++++++++++++++++++------------------ 1 file changed, 25 insertions(+), 23 deletions(-) (limited to 'src/transformations/deformation.cpp') diff --git a/src/transformations/deformation.cpp b/src/transformations/deformation.cpp index 993647fc..778c1a2e 100644 --- a/src/transformations/deformation.cpp +++ b/src/transformations/deformation.cpp @@ -81,13 +81,16 @@ struct deformationData { // --------------------------------------------------------------------------- static const GenericShiftGrid* findGrid(const ListOfGenericGrids& grids, - const PJ_LP& input) + const PJ_LP& input, + GenericShiftGridSet *&gridSetOut) { for( const auto& gridset: grids ) { auto grid = gridset->gridAt(input.lam, input.phi); - if( grid ) + if( grid ) { + gridSetOut = gridset.get(); return grid; + } } return nullptr; } @@ -101,7 +104,8 @@ static bool get_grid_values(PJ* P, double& vy, double& vz) { - auto grid = findGrid(Q->grids, lp); + GenericShiftGridSet* gridset = nullptr; + auto grid = findGrid(Q->grids, lp, gridset); if( !grid ) { return false; } @@ -145,30 +149,28 @@ static bool get_grid_values(PJ* P, 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) || + bool error =( !grid->valueAt(ix, iy, sampleE, dx1) || + !grid->valueAt(ix, iy, sampleN, dy1) || + !grid->valueAt(ix, iy, sampleU, dz1) || + !grid->valueAt(ix2, iy, sampleE, dx2) || + !grid->valueAt(ix2, iy, sampleN, dy2) || + !grid->valueAt(ix2, iy, sampleU, dz2) || + !grid->valueAt(ix, iy2, sampleE, dx3) || + !grid->valueAt(ix, iy2, sampleN, dy3) || + !grid->valueAt(ix, iy2, sampleU, dz3) || + !grid->valueAt(ix2, iy2, sampleE, dx4) || !grid->valueAt(ix2, iy2, sampleN, dy4) || - !grid->valueAt(ix2, iy2, sampleU, dz4) ) { + !grid->valueAt(ix2, iy2, sampleU, dz4) ); + if( grid->hasChanged() ) { + if( gridset->reopen(P->ctx) ) { + return get_grid_values( P, Q, lp, vx, vy, vz); + } + error = true; + } + if( error ) { return false; } -- cgit v1.2.3 From 408cb1e4279786e3729e942c2f8b2ee2ad435c94 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 7 Jan 2020 15:01:36 +0100 Subject: deformation.cpp and xyzgridshift.cpp: factor out common code (addresses https://github.com/OSGeo/PROJ/pull/1790#discussion_r363748615) --- src/transformations/deformation.cpp | 74 ++++++------------------------------- 1 file changed, 11 insertions(+), 63 deletions(-) (limited to 'src/transformations/deformation.cpp') diff --git a/src/transformations/deformation.cpp b/src/transformations/deformation.cpp index 778c1a2e..8aee50c9 100644 --- a/src/transformations/deformation.cpp +++ b/src/transformations/deformation.cpp @@ -80,23 +80,6 @@ struct deformationData { // --------------------------------------------------------------------------- -static const GenericShiftGrid* findGrid(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; -} - -// --------------------------------------------------------------------------- - static bool get_grid_values(PJ* P, deformationData* Q, const PJ_LP& lp, @@ -105,7 +88,7 @@ static bool get_grid_values(PJ* P, double& vz) { GenericShiftGridSet* gridset = nullptr; - auto grid = findGrid(Q->grids, lp, gridset); + auto grid = pj_find_generic_grid(Q->grids, lp, gridset); if( !grid ) { return false; } @@ -140,55 +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(grid_x); - int iy = static_cast(grid_y); - int ix2 = std::min(ix + 1, grid->width() - 1); - int iy2 = std::min(iy + 1, grid->height() - 1); - - float dx1, dy1, dz1; - float dx2, dy2, dz2; - float dx3, dy3, dz3; - float dx4, dy4, dz4; - bool error =( !grid->valueAt(ix, iy, sampleE, dx1) || - !grid->valueAt(ix, iy, sampleN, dy1) || - !grid->valueAt(ix, iy, sampleU, dz1) || - !grid->valueAt(ix2, iy, sampleE, dx2) || - !grid->valueAt(ix2, iy, sampleN, dy2) || - !grid->valueAt(ix2, iy, sampleU, dz2) || - !grid->valueAt(ix, iy2, sampleE, dx3) || - !grid->valueAt(ix, iy2, sampleN, dy3) || - !grid->valueAt(ix, iy2, sampleU, dz3) || - !grid->valueAt(ix2, iy2, sampleE, dx4) || - !grid->valueAt(ix2, iy2, sampleN, dy4) || - !grid->valueAt(ix2, iy2, sampleU, dz4) ); - if( grid->hasChanged() ) { - if( gridset->reopen(P->ctx) ) { + 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); - } - error = true; - } - 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; // 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; } -- cgit v1.2.3