aboutsummaryrefslogtreecommitdiff
path: root/src/transformations
diff options
context:
space:
mode:
Diffstat (limited to 'src/transformations')
-rw-r--r--src/transformations/deformation.cpp86
-rw-r--r--src/transformations/hgridshift.cpp10
-rw-r--r--src/transformations/vgridshift.cpp10
-rw-r--r--src/transformations/xyzgridshift.cpp80
4 files changed, 43 insertions, 143 deletions
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).");
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-15 22:26:15 +0000