diff options
| -rwxr-xr-x | scripts/reformat_cpp.sh | 4 | ||||
| -rw-r--r-- | src/4D_api.cpp | 43 | ||||
| -rw-r--r-- | src/Makefile.am | 5 | ||||
| -rw-r--r-- | src/apply_vgridshift.cpp | 387 | ||||
| -rw-r--r-- | src/gridinfo.cpp | 154 | ||||
| -rw-r--r-- | src/grids.cpp | 305 | ||||
| -rw-r--r-- | src/grids.hpp | 103 | ||||
| -rw-r--r-- | src/init.cpp | 3 | ||||
| -rw-r--r-- | src/lib_proj.cmake | 2 | ||||
| -rw-r--r-- | src/malloc.cpp | 1 | ||||
| -rw-r--r-- | src/proj_internal.h | 10 | ||||
| -rw-r--r-- | src/transform.cpp | 69 | ||||
| -rw-r--r-- | src/transformations/vgridshift.cpp | 4 | ||||
| -rw-r--r-- | test/gie/4D-API_cs2cs-style.gie | 14 |
14 files changed, 661 insertions, 443 deletions
diff --git a/scripts/reformat_cpp.sh b/scripts/reformat_cpp.sh index a8002a6c..899e665f 100755 --- a/scripts/reformat_cpp.sh +++ b/scripts/reformat_cpp.sh @@ -15,7 +15,9 @@ esac TOPDIR="$SCRIPT_DIR/.." -for i in "$TOPDIR"/include/proj/*.hpp "$TOPDIR"/include/proj/internal/*.hpp "$TOPDIR"/src/iso19111/*.cpp "$TOPDIR"/test/unit/*.cpp "$TOPDIR"/src/apps/projinfo.cpp "$TOPDIR"/src/tracing.cpp; do +for i in "$TOPDIR"/include/proj/*.hpp "$TOPDIR"/include/proj/internal/*.hpp \ + "$TOPDIR"/src/iso19111/*.cpp "$TOPDIR"/test/unit/*.cpp "$TOPDIR"/src/apps/projinfo.cpp \ + "$TOPDIR"/src/tracing.cpp "$TOPDIR"/src/grids.hpp "$TOPDIR"/src/grids.cpp; do if ! echo "$i" | grep -q "lru_cache.hpp"; then "$SCRIPT_DIR"/reformat.sh "$i"; fi diff --git a/src/4D_api.cpp b/src/4D_api.cpp index f37594b5..ed6a62d6 100644 --- a/src/4D_api.cpp +++ b/src/4D_api.cpp @@ -1558,6 +1558,45 @@ PJ_GRID_INFO proj_grid_info(const char *gridname) { /* in case the grid wasn't found */ if (gridinfo->filename == nullptr || gridinfo->ct == nullptr) { + + const auto gridSet = NS_PROJ::VerticalShiftGridSet::open(ctx, gridname); + if( gridSet ) + { + const auto& grids = gridSet->grids(); + if( !grids.empty() ) + { + const auto& grid = grids.front(); + const auto& extent = grid->extentAndRes(); + + /* name of grid */ + strncpy (grinfo.gridname, gridname, sizeof(grinfo.gridname) - 1); + + /* full path of grid */ + pj_find_file(ctx, gridname, grinfo.filename, sizeof(grinfo.filename) - 1); + + /* grid format */ + strncpy (grinfo.format, gridSet->format().c_str(), sizeof(grinfo.format) - 1); + + /* grid size */ + grinfo.n_lon = grid->width(); + grinfo.n_lat = grid->height(); + + /* cell size */ + grinfo.cs_lon = extent.resLon; + grinfo.cs_lat = extent.resLat; + + /* bounds of grid */ + grinfo.lowerleft.lam = extent.westLon; + grinfo.lowerleft.phi = extent.southLat; + grinfo.upperright.lam = extent.eastLon; + grinfo.upperright.phi = extent.northLat; + + pj_gridinfo_free(ctx, gridinfo); + + return grinfo; + } + } + pj_gridinfo_free(ctx, gridinfo); strcpy(grinfo.format, "missing"); return grinfo; @@ -1585,8 +1624,8 @@ PJ_GRID_INFO proj_grid_info(const char *gridname) { /* bounds of grid */ grinfo.lowerleft = gridinfo->ct->ll; - grinfo.upperright.lam = grinfo.lowerleft.lam + grinfo.n_lon*grinfo.cs_lon; - grinfo.upperright.phi = grinfo.lowerleft.phi + grinfo.n_lat*grinfo.cs_lat; + grinfo.upperright.lam = grinfo.lowerleft.lam + (grinfo.n_lon-1)*grinfo.cs_lon; + grinfo.upperright.phi = grinfo.lowerleft.phi + (grinfo.n_lat-1)*grinfo.cs_lat; pj_gridinfo_free(ctx, gridinfo); diff --git a/src/Makefile.am b/src/Makefile.am index 80596ef8..3d0868a7 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -212,7 +212,10 @@ libproj_la_SOURCES = \ proj_json_streaming_writer.hpp \ proj_json_streaming_writer.cpp \ \ - tracing.cpp + tracing.cpp \ + \ + grids.hpp \ + grids.cpp # The sed hack is to please MSVC diff --git a/src/apply_vgridshift.cpp b/src/apply_vgridshift.cpp index daa44858..0d18b5dc 100644 --- a/src/apply_vgridshift.cpp +++ b/src/apply_vgridshift.cpp @@ -26,7 +26,9 @@ * DEALINGS IN THE SOFTWARE. *****************************************************************************/ -#define PJ_LIB__ +#ifndef FROM_PROJ_CPP +#define FROM_PROJ_CPP +#endif #include <assert.h> #include <stdio.h> @@ -34,273 +36,116 @@ #include <math.h> #include "proj_internal.h" +#include "proj/internal/internal.hpp" -static int is_nodata(float value, double vmultiplier) -{ - /* nodata? */ - /* GTX official nodata value if -88.88880f, but some grids also */ - /* use other big values for nodata (e.g naptrans2008.gtx has */ - /* nodata values like -2147479936), so test them too */ - return value * vmultiplier > 1000 || value * vmultiplier < -1000 || value == -88.88880f; -} +using namespace NS_PROJ; + +static double read_vgrid_value( PJ *defn, PJ_LP input, double vmultiplier) { -static double read_vgrid_value( PJ *defn, PJ_LP input, double vmultiplier, int *gridlist_count_p, PJ_GRIDINFO **tables, struct CTABLE *ct) { - int itable = 0; - double value = HUGE_VAL; - double grid_x, grid_y; - long grid_ix, grid_iy; - long grid_ix2, grid_iy2; - float *cvs; /* do not deal with NaN coordinates */ /* cppcheck-suppress duplicateExpression */ if( isnan(input.phi) || isnan(input.lam) ) - itable = *gridlist_count_p; - - /* keep trying till we find a table that works */ - for ( ; itable < *gridlist_count_p; itable++ ) { - PJ_GRIDINFO *gi = tables[itable]; - - ct = gi->ct; - - /* skip tables that don't match our point at all (latitude check). */ - if( ct->ll.phi > input.phi - || ct->ll.phi + (ct->lim.phi-1) * ct->del.phi < input.phi ) - continue; - - bool fullWorldLongExtent = false; - if( fabs(ct->lim.lam * ct->del.lam - 2 * M_PI) < 1e-10 ) - { - fullWorldLongExtent = true; - } - - /* skip tables that don't match our point at all (longitude check). */ - else if( ct->ll.lam > input.lam - || ct->ll.lam + (ct->lim.lam-1) * ct->del.lam < input.lam ) - continue; - - /* If we have child nodes, check to see if any of them apply. */ - while( gi->child != nullptr ) - { - PJ_GRIDINFO *child; - - for( child = gi->child; child != nullptr; child = child->next ) - { - struct CTABLE *ct1 = child->ct; - - fullWorldLongExtent = false; - - if( ct1->ll.phi > input.phi - || ct1->ll.phi+(ct1->lim.phi-1)*ct1->del.phi < input.phi) - continue; - - if( fabs(ct1->lim.lam * ct1->del.lam - 2 * M_PI) < 1e-10 ) - { - fullWorldLongExtent = true; - } - else if( ct1->ll.lam > input.lam - || ct1->ll.lam+(ct1->lim.lam-1)*ct1->del.lam < input.lam) - continue; - - break; - } - - /* we didn't find a more refined child node to use, so go with current grid */ - if( child == nullptr ) - { - break; - } - - /* Otherwise let's try for childrens children .. */ - gi = child; - ct = child->ct; - } - - /* load the grid shift info if we don't have it. */ - if( ct->cvs == nullptr ) - { - pj_gridinfo_load( pj_get_ctx(defn), gi ); - } - if( ct->cvs == nullptr ) - { - pj_ctx_set_errno( defn->ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); - return PJD_ERR_FAILED_TO_LOAD_GRID; - } + return HUGE_VAL; + } - /* Interpolation a location within the grid */ - grid_x = (input.lam - ct->ll.lam) / ct->del.lam; - if( fullWorldLongExtent ) { - // 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 + ct->lim.lam, ct->lim.lam) + ct->lim.lam, - ct->lim.lam); - } - grid_y = (input.phi - ct->ll.phi) / ct->del.phi; - grid_ix = lround(floor(grid_x)); - assert(grid_ix >= 0 && grid_ix < ct->lim.lam); - grid_iy = lround(floor(grid_y)); - assert(grid_iy >= 0 && grid_iy < ct->lim.phi); - grid_x -= grid_ix; - grid_y -= grid_iy; + const VerticalShiftGrid* grid = nullptr; + for( const auto& gridset: defn->vgrids ) + { + grid = gridset->gridAt(input.lam, input.phi); + if( grid ) + break; + } + if( !grid ) + { + return HUGE_VAL; + } - grid_ix2 = grid_ix + 1; - if( grid_ix2 >= ct->lim.lam ) { - if( fullWorldLongExtent ) { - grid_ix2 = 0; - } else { - grid_ix2 = ct->lim.lam - 1; - } - } - grid_iy2 = grid_iy + 1; - if( grid_iy2 >= ct->lim.phi ) - grid_iy2 = ct->lim.phi - 1; + const auto& extent = grid->extentAndRes(); - cvs = (float *) ct->cvs; - { - float value_a = cvs[grid_ix + grid_iy * ct->lim.lam]; - float value_b = cvs[grid_ix2 + grid_iy * ct->lim.lam]; - float value_c = cvs[grid_ix + grid_iy2 * ct->lim.lam]; - float value_d = cvs[grid_ix2 + grid_iy2 * ct->lim.lam]; - double total_weight = 0.0; - int n_weights = 0; - value = 0.0f; - if( !is_nodata(value_a, vmultiplier) ) - { - double weight = (1.0-grid_x) * (1.0-grid_y); - value += value_a * weight; - total_weight += weight; - n_weights ++; - } - if( !is_nodata(value_b, vmultiplier) ) - { - double weight = (grid_x) * (1.0-grid_y); - value += value_b * weight; - total_weight += weight; - n_weights ++; - } - if( !is_nodata(value_c, vmultiplier) ) - { - double weight = (1.0-grid_x) * (grid_y); - value += value_c * weight; - total_weight += weight; - n_weights ++; - } - if( !is_nodata(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; + /* 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; } - return value * vmultiplier; -} - -/************************************************************************/ -/* 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. */ -/************************************************************************/ -int pj_apply_vgridshift( PJ *defn, const char *listname, - PJ_GRIDINFO ***gridlist_p, - int *gridlist_count_p, - int inverse, - long point_count, int point_offset, - double *x, double *y, double *z ) - -{ - int i; - static int debug_count = 0; - PJ_GRIDINFO **tables; - struct CTABLE ct; + double total_weight = 0.0; + int n_weights = 0; + double value = 0.0f; - if( *gridlist_p == nullptr ) + if( !grid->isNodata(value_a, vmultiplier) ) { - *gridlist_p = - pj_gridlist_from_nadgrids( pj_get_ctx(defn), - pj_param(defn->ctx,defn->params,listname).s, - gridlist_count_p ); - - if( *gridlist_p == nullptr || *gridlist_count_p == 0 ) - return defn->ctx->last_errno; + double weight = (1.0-grid_x) * (1.0-grid_y); + value += value_a * weight; + total_weight += weight; + n_weights ++; } - - if( *gridlist_count_p == 0 ) + if( !grid->isNodata(value_b, vmultiplier) ) { - pj_ctx_set_errno( defn->ctx, PJD_ERR_FAILED_TO_LOAD_GRID); - return PJD_ERR_FAILED_TO_LOAD_GRID; + double weight = (grid_x) * (1.0-grid_y); + value += value_b * weight; + total_weight += weight; + n_weights ++; } - - tables = *gridlist_p; - defn->ctx->last_errno = 0; - - for( i = 0; i < point_count; i++ ) + if( !grid->isNodata(value_c, vmultiplier) ) { - double value; - long io = i * point_offset; - PJ_LP input; - - input.phi = y[io]; - input.lam = x[io]; - - value = read_vgrid_value(defn, input, 1.0, gridlist_count_p, tables, &ct); - - if( inverse ) - z[io] -= value; - else - z[io] += value; - if( value != HUGE_VAL ) - { - if( debug_count++ < 20 ) { - proj_log_trace(defn, "pj_apply_gridshift(): used %s", ct.id); - break; - } - } - - if( value == HUGE_VAL ) - { - int itable; - 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( itable = 0; itable < *gridlist_count_p; itable++ ) - { - PJ_GRIDINFO *gi = tables[itable]; - if( itable == 0 ) - gridlist += " tried: "; - else - gridlist += ','; - gridlist += gi->gridname; - } - - proj_log_debug(defn, "%s", gridlist.c_str()); - pj_ctx_set_errno( defn->ctx, PJD_ERR_GRID_AREA ); - - return PJD_ERR_GRID_AREA; - } + 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 0; + return value * vmultiplier; } /**********************************************/ -int proj_vgrid_init(PJ* P, const char *grids) { +int proj_vgrid_init(PJ* P, const char *gridkey) { /********************************************** Initizalize and populate gridlist. @@ -314,29 +159,39 @@ int proj_vgrid_init(PJ* P, const char *grids) { ***********************************************/ - /* prepend "s" to the "grids" string to allow usage with pj_param */ - char *sgrids = (char *) pj_malloc( (strlen(grids)+1+1) *sizeof(char) ); - sprintf(sgrids, "%s%s", "s", grids); - - if (P->vgridlist_geoid == nullptr) { - P->vgridlist_geoid = pj_gridlist_from_nadgrids( - P->ctx, - pj_param(P->ctx, P->params, sgrids).s, - &(P->vgridlist_geoid_count) - ); + std::string key("s"); + key += gridkey; + const char* grids = pj_param(P->ctx, P->params, key.c_str()).s; + if( grids == nullptr ) + return 0; - if( P->vgridlist_geoid == nullptr || P->vgridlist_geoid_count == 0 ) { - pj_dealloc(sgrids); - return 0; + 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 = VerticalShiftGridSet::open(P->ctx, gridname); + if( !gridSet ) + { + if( !canFail ) + { + pj_ctx_set_errno( P->ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); + P->vgrids.clear(); + return 0; + } + } + else + { + P->vgrids.emplace_back(std::move(gridSet)); } } - if (P->vgridlist_geoid_count == 0) { - proj_errno_set(P, PJD_ERR_FAILED_TO_LOAD_GRID); - } - - pj_dealloc(sgrids); - return P->vgridlist_geoid_count; + return static_cast<int>(P->vgrids.size()); } /***********************************************/ @@ -350,11 +205,9 @@ double proj_vgrid_value(PJ *P, PJ_LP lp, double vmultiplier){ ************************************************/ - struct CTABLE used_grid; double value; - memset(&used_grid, 0, sizeof(struct CTABLE)); - value = read_vgrid_value(P, lp, vmultiplier, &(P->vgridlist_geoid_count), P->vgridlist_geoid, &used_grid); + value = read_vgrid_value(P, 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; diff --git a/src/gridinfo.cpp b/src/gridinfo.cpp index 7f7d930f..01d392b3 100644 --- a/src/gridinfo.cpp +++ b/src/gridinfo.cpp @@ -352,51 +352,6 @@ int pj_gridinfo_load( projCtx_t* ctx, PJ_GRIDINFO *gi ) return 1; } -/* -------------------------------------------------------------------- */ -/* GTX format. */ -/* -------------------------------------------------------------------- */ - else if( strcmp(gi->format,"gtx") == 0 ) - { - int words = gi->ct->lim.lam * gi->ct->lim.phi; - PAFile fid; - - fid = pj_open_lib( ctx, gi->filename, "rb" ); - - if( fid == nullptr ) - { - pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); - pj_release_lock(); - return 0; - } - - pj_ctx_fseek( ctx, fid, gi->grid_offset, SEEK_SET ); - - ct_tmp.cvs = (FLP *) pj_malloc(words*sizeof(float)); - if( ct_tmp.cvs == nullptr ) - { - pj_ctx_set_errno( ctx, ENOMEM ); - pj_release_lock(); - return 0; - } - - if( pj_ctx_fread( ctx, ct_tmp.cvs, sizeof(float), words, fid ) - != (size_t)words ) - { - pj_dalloc( ct_tmp.cvs ); - pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); - pj_release_lock(); - return 0; - } - - if( IS_LSB ) - swap_words( (unsigned char *) ct_tmp.cvs, 4, words ); - - pj_ctx_fclose( ctx, fid ); - gi->ct->cvs = ct_tmp.cvs; - pj_release_lock(); - return 1; - } - else { pj_release_lock(); @@ -737,108 +692,6 @@ static int pj_gridinfo_init_ntv1( projCtx ctx, PAFile fid, PJ_GRIDINFO *gi ) } /************************************************************************/ -/* pj_gridinfo_init_gtx() */ -/* */ -/* Load a NOAA .gtx vertical datum shift file. */ -/************************************************************************/ - -static int pj_gridinfo_init_gtx( projCtx ctx, PAFile fid, PJ_GRIDINFO *gi ) - -{ - unsigned char header[40]; - struct CTABLE *ct; - double xorigin,yorigin,xstep,ystep; - int rows, columns; - - /* cppcheck-suppress sizeofCalculation */ - STATIC_ASSERT( sizeof(pj_int32) == 4 ); - /* cppcheck-suppress sizeofCalculation */ - STATIC_ASSERT( sizeof(double) == 8 ); - -/* -------------------------------------------------------------------- */ -/* Read the header. */ -/* -------------------------------------------------------------------- */ - if( pj_ctx_fread( ctx, header, sizeof(header), 1, fid ) != 1 ) - { - pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); - return 0; - } - -/* -------------------------------------------------------------------- */ -/* Regularize fields of interest and extract. */ -/* -------------------------------------------------------------------- */ - if( IS_LSB ) - { - swap_words( header+0, 8, 4 ); - swap_words( header+32, 4, 2 ); - } - - memcpy( &yorigin, header+0, 8 ); - memcpy( &xorigin, header+8, 8 ); - memcpy( &ystep, header+16, 8 ); - memcpy( &xstep, header+24, 8 ); - - memcpy( &rows, header+32, 4 ); - memcpy( &columns, header+36, 4 ); - - if( xorigin < -360 || xorigin > 360 - || yorigin < -90 || yorigin > 90 ) - { - pj_log( ctx, PJ_LOG_ERROR, - "gtx file header has invalid extents, corrupt?"); - pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); - return 0; - } - -/* -------------------------------------------------------------------- */ -/* Fill in CTABLE structure. */ -/* -------------------------------------------------------------------- */ - ct = (struct CTABLE *) pj_malloc(sizeof(struct CTABLE)); - if (!ct) { - pj_ctx_set_errno(ctx, ENOMEM); - return 0; - } - strcpy( ct->id, "GTX Vertical Grid Shift File" ); - - ct->ll.lam = xorigin; - ct->ll.phi = yorigin; - ct->del.lam = xstep; - ct->del.phi = ystep; - ct->lim.lam = columns; - ct->lim.phi = rows; - - /* some GTX files come in 0-360 and we shift them back into the - expected -180 to 180 range if possible. This does not solve - problems with grids spanning the dateline. */ - if( ct->ll.lam >= 180.0 ) - ct->ll.lam -= 360.0; - - if( ct->ll.lam >= 0.0 && ct->ll.lam + ct->del.lam * ct->lim.lam > 180.0 ) - { - pj_log( ctx, PJ_LOG_DEBUG_MAJOR, - "This GTX spans the dateline! This will cause problems." ); - } - - pj_log( ctx, PJ_LOG_DEBUG_MINOR, - "GTX %dx%d: LL=(%.9g,%.9g) UR=(%.9g,%.9g)", - ct->lim.lam, ct->lim.phi, - ct->ll.lam, ct->ll.phi, - ct->ll.lam + (columns-1)*xstep, ct->ll.phi + (rows-1)*ystep); - - ct->ll.lam *= DEG_TO_RAD; - ct->ll.phi *= DEG_TO_RAD; - ct->del.lam *= DEG_TO_RAD; - ct->del.phi *= DEG_TO_RAD; - ct->cvs = nullptr; - - gi->ct = ct; - gi->grid_offset = 40; - gi->format = "gtx"; - - return 1; -} - -/************************************************************************/ /* pj_gridinfo_init() */ /* */ /* Open and parse header details from a datum gridshift file */ @@ -929,13 +782,6 @@ PJ_GRIDINFO *pj_gridinfo_init( projCtx ctx, const char *gridname ) pj_gridinfo_init_ntv2( ctx, fp, gilist ); } - else if( strlen(gridname) > 4 - && (strcmp(gridname+strlen(gridname)-3,"gtx") == 0 - || strcmp(gridname+strlen(gridname)-3,"GTX") == 0) ) - { - pj_gridinfo_init_gtx( ctx, fp, gilist ); - } - else if( header_size >= 9 && strncmp(header + 0,"CTABLE V2",9) == 0 ) { struct CTABLE *ct = nad_ctable2_init( ctx, (struct projFileAPI_t*)fp ); diff --git a/src/grids.cpp b/src/grids.cpp new file mode 100644 index 00000000..1a18f975 --- /dev/null +++ b/src/grids.cpp @@ -0,0 +1,305 @@ +/****************************************************************************** + * Project: PROJ + * Purpose: Grid management + * Author: Even Rouault, <even.rouault at spatialys.com> + * + ****************************************************************************** + * Copyright (c) 2000, Frank Warmerdam <warmerdam@pobox.com> + * Copyright (c) 2019, Even Rouault, <even.rouault at spatialys.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 "grids.hpp" +#include "proj/internal/internal.hpp" +#include "proj_internal.h" + +NS_PROJ_START + +using namespace internal; + +/************************************************************************/ +/* swap_words() */ +/* */ +/* Convert the byte order of the given word(s) in place. */ +/************************************************************************/ + +static const int byte_order_test = 1; +#define IS_LSB (1 == ((const unsigned char *)(&byte_order_test))[0]) + +static void swap_words(void *dataIn, size_t word_size, size_t word_count) + +{ + unsigned char *data = static_cast<unsigned char *>(dataIn); + for (size_t word = 0; word < word_count; word++) { + for (size_t i = 0; i < word_size / 2; i++) { + unsigned char t; + + t = data[i]; + data[i] = data[word_size - i - 1]; + data[word_size - i - 1] = t; + } + + data += word_size; + } +} + +bool ExtentAndRes::fullWorldLongitude() const { + return eastLon - westLon + resLon >= 2 * M_PI - 1e-10; +} + +// --------------------------------------------------------------------------- + +VerticalShiftGrid::VerticalShiftGrid(int widthIn, int heightIn, + const ExtentAndRes &extentIn) + : m_width(widthIn), m_height(heightIn), m_extent(extentIn) {} + +// --------------------------------------------------------------------------- + +VerticalShiftGrid::~VerticalShiftGrid() = default; + +// --------------------------------------------------------------------------- + +static ExtentAndRes globalExtent() { + ExtentAndRes extent; + extent.westLon = -M_PI; + extent.southLat = -M_PI / 2; + extent.eastLon = M_PI; + extent.northLat = M_PI / 2; + extent.resLon = M_PI; + extent.resLat = M_PI / 2; + return extent; +} + +// --------------------------------------------------------------------------- + +class NullVerticalShiftGrid : public VerticalShiftGrid { + + public: + NullVerticalShiftGrid() : VerticalShiftGrid(3, 3, globalExtent()) {} + + bool valueAt(int, int, float &out) const override; + bool isNodata(float, double) const override { return false; } +}; + +// --------------------------------------------------------------------------- + +bool NullVerticalShiftGrid::valueAt(int, int, float &out) const { + out = 0.0f; + return true; +} + +// --------------------------------------------------------------------------- + +class GTXVerticalShiftGrid : public VerticalShiftGrid { + PJ_CONTEXT *m_ctx; + PAFile m_fp; + + GTXVerticalShiftGrid(const GTXVerticalShiftGrid &) = delete; + GTXVerticalShiftGrid &operator=(const GTXVerticalShiftGrid &) = delete; + + public: + GTXVerticalShiftGrid(PJ_CONTEXT *ctx, PAFile fp, int widthIn, int heightIn, + const ExtentAndRes &extentIn) + : VerticalShiftGrid(widthIn, heightIn, extentIn), m_ctx(ctx), m_fp(fp) { + } + + ~GTXVerticalShiftGrid() override; + + bool valueAt(int x, int y, float &out) const override; + bool isNodata(float val, double multiplier) const override; + + static GTXVerticalShiftGrid *open(PJ_CONTEXT *ctx, PAFile fp); +}; + +// --------------------------------------------------------------------------- + +GTXVerticalShiftGrid::~GTXVerticalShiftGrid() { pj_ctx_fclose(m_ctx, m_fp); } + +// --------------------------------------------------------------------------- + +GTXVerticalShiftGrid *GTXVerticalShiftGrid::open(PJ_CONTEXT *ctx, PAFile fp) { + unsigned char header[40]; + + /* -------------------------------------------------------------------- */ + /* Read the header. */ + /* -------------------------------------------------------------------- */ + if (pj_ctx_fread(ctx, header, sizeof(header), 1, fp) != 1) { + pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID); + return nullptr; + } + + /* -------------------------------------------------------------------- */ + /* Regularize fields of interest and extract. */ + /* -------------------------------------------------------------------- */ + if (IS_LSB) { + swap_words(header + 0, 8, 4); + swap_words(header + 32, 4, 2); + } + + double xorigin, yorigin, xstep, ystep; + int rows, columns; + + memcpy(&yorigin, header + 0, 8); + memcpy(&xorigin, header + 8, 8); + memcpy(&ystep, header + 16, 8); + memcpy(&xstep, header + 24, 8); + + memcpy(&rows, header + 32, 4); + memcpy(&columns, header + 36, 4); + + if (xorigin < -360 || xorigin > 360 || yorigin < -90 || yorigin > 90) { + pj_log(ctx, PJ_LOG_ERROR, + "gtx file header has invalid extents, corrupt?"); + pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID); + return nullptr; + } + + /* some GTX files come in 0-360 and we shift them back into the + expected -180 to 180 range if possible. This does not solve + problems with grids spanning the dateline. */ + if (xorigin >= 180.0) + xorigin -= 360.0; + + if (xorigin >= 0.0 && xorigin + xstep * columns > 180.0) { + pj_log(ctx, PJ_LOG_DEBUG_MAJOR, + "This GTX spans the dateline! This will cause problems."); + } + + ExtentAndRes extent; + extent.westLon = xorigin * DEG_TO_RAD; + extent.southLat = yorigin * DEG_TO_RAD; + extent.resLon = xstep * DEG_TO_RAD; + extent.resLat = ystep * DEG_TO_RAD; + extent.eastLon = (xorigin + xstep * (columns - 1)) * DEG_TO_RAD; + extent.northLat = (yorigin + ystep * (rows - 1)) * DEG_TO_RAD; + + return new GTXVerticalShiftGrid(ctx, fp, columns, rows, extent); +} + +// --------------------------------------------------------------------------- + +bool GTXVerticalShiftGrid::valueAt(int x, int y, float &out) const { + assert(x >= 0 && y >= 0 && x < m_width && y < m_height); + + pj_ctx_fseek(m_ctx, m_fp, 40 + sizeof(float) * (y * m_width + x), SEEK_SET); + if (pj_ctx_fread(m_ctx, &out, sizeof(out), 1, m_fp) != 1) { + pj_ctx_set_errno(m_ctx, PJD_ERR_FAILED_TO_LOAD_GRID); + return false; + } + if (IS_LSB) { + swap_words(&out, sizeof(float), 1); + } + return true; +} + +// --------------------------------------------------------------------------- + +bool GTXVerticalShiftGrid::isNodata(float val, double multiplier) const { + /* nodata? */ + /* GTX official nodata value if -88.88880f, but some grids also */ + /* use other big values for nodata (e.g naptrans2008.gtx has */ + /* nodata values like -2147479936), so test them too */ + return val * multiplier > 1000 || val * multiplier < -1000 || + val == -88.88880f; +} + +// --------------------------------------------------------------------------- + +VerticalShiftGridSet::VerticalShiftGridSet() = default; + +// --------------------------------------------------------------------------- + +VerticalShiftGridSet::~VerticalShiftGridSet() = default; + +// --------------------------------------------------------------------------- + +std::unique_ptr<VerticalShiftGridSet> +VerticalShiftGridSet::open(PJ_CONTEXT *ctx, const std::string &filename) { + if (filename == "null") { + auto set = + std::unique_ptr<VerticalShiftGridSet>(new VerticalShiftGridSet()); + set->m_name = filename; + set->m_format = "null"; + set->m_grids.push_back(std::unique_ptr<NullVerticalShiftGrid>( + new NullVerticalShiftGrid())); + return set; + } + + PAFile fp; + if (!(fp = pj_open_lib(ctx, filename.c_str(), "rb"))) { + ctx->last_errno = 0; /* don't treat as a persistent error */ + return nullptr; + } + if (ends_with(filename, "gtx") || ends_with(filename, "GTX")) { + auto grid = GTXVerticalShiftGrid::open(ctx, fp); + if (!grid) { + pj_ctx_fclose(ctx, fp); + return nullptr; + } + auto set = + std::unique_ptr<VerticalShiftGridSet>(new VerticalShiftGridSet()); + set->m_name = filename; + set->m_format = "gtx"; + set->m_grids.push_back(std::unique_ptr<VerticalShiftGrid>(grid)); + return set; + } + + pj_log(ctx, PJ_LOG_DEBUG_MAJOR, "Unrecognized vertical grid format"); + pj_ctx_fclose(ctx, fp); + return nullptr; +} + +// --------------------------------------------------------------------------- + +const VerticalShiftGrid *VerticalShiftGrid::gridAt(double lon, + double lat) const { + for (const auto &child : m_children) { + const auto &extentChild = child->extentAndRes(); + if ((extentChild.fullWorldLongitude() || + (lon >= extentChild.westLon && lon <= extentChild.eastLon)) && + lat >= extentChild.southLat && lat <= extentChild.northLat) { + return child->gridAt(lon, lat); + } + } + return this; +} +// --------------------------------------------------------------------------- + +const VerticalShiftGrid *VerticalShiftGridSet::gridAt(double lon, + double lat) const { + for (const auto &grid : m_grids) { + if (dynamic_cast<NullVerticalShiftGrid *>(grid.get())) { + return grid.get(); + } + const auto &extent = grid->extentAndRes(); + if ((extent.fullWorldLongitude() || + (lon >= extent.westLon && lon <= extent.eastLon)) && + lat >= extent.southLat && lat <= extent.northLat) { + return grid->gridAt(lon, lat); + } + } + return nullptr; +} + +NS_PROJ_END diff --git a/src/grids.hpp b/src/grids.hpp new file mode 100644 index 00000000..0e3ccffb --- /dev/null +++ b/src/grids.hpp @@ -0,0 +1,103 @@ +/****************************************************************************** + * Project: PROJ + * Purpose: Grid management + * Author: Even Rouault, <even.rouault at spatialys.com> + * + ****************************************************************************** + * Copyright (c) 2019, Even Rouault, <even.rouault at spatialys.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 GRIDS_HPP_INCLUDED +#define GRIDS_HPP_INCLUDED + +#include <memory> +#include <vector> + +#include "proj/util.hpp" + +typedef struct projCtx_t PJ_CONTEXT; + +NS_PROJ_START + +struct ExtentAndRes { + double westLon; // in radian + double southLat; // in radian + double eastLon; // in radian + double northLat; // in radian + double resLon; // in radian + double resLat; // in radian + + bool fullWorldLongitude() const; +}; + +// --------------------------------------------------------------------------- + +class VerticalShiftGrid { + protected: + int m_width; + int m_height; + ExtentAndRes m_extent; + std::vector<std::unique_ptr<VerticalShiftGrid>> m_children{}; + + public: + VerticalShiftGrid(int widthIn, int heightIn, const ExtentAndRes &extentIn); + virtual ~VerticalShiftGrid(); + + int width() const { return m_width; } + int height() const { return m_height; } + const ExtentAndRes &extentAndRes() const { return m_extent; } + + const VerticalShiftGrid *gridAt(double lon, double lat) const; + + virtual bool isNodata(float /*val*/, double /* multiplier */) const { + return false; + } + + // x = 0 is western-most column, y = 0 is southern-most line + virtual bool valueAt(int x, int y, float &out) const = 0; +}; + +// --------------------------------------------------------------------------- + +class VerticalShiftGridSet { + std::string m_name{}; + std::string m_format{}; + std::vector<std::unique_ptr<VerticalShiftGrid>> m_grids{}; + + VerticalShiftGridSet(); + + public: + virtual ~VerticalShiftGridSet(); + + static std::unique_ptr<VerticalShiftGridSet> + open(PJ_CONTEXT *ctx, const std::string &filename); + + const std::string &name() const { return m_name; } + const std::string &format() const { return m_format; } + const std::vector<std::unique_ptr<VerticalShiftGrid>> &grids() const { + return m_grids; + } + const VerticalShiftGrid *gridAt(double lon, double lat) const; +}; + +NS_PROJ_END + +#endif // GRIDS_HPP_INCLUDED diff --git a/src/init.cpp b/src/init.cpp index 1ed46e5a..69b4ec5f 100644 --- a/src/init.cpp +++ b/src/init.cpp @@ -649,9 +649,6 @@ pj_init_ctx_with_allow_init_epsg(projCtx ctx, int argc, char **argv, int allow_i PIN->gridlist = nullptr; PIN->gridlist_count = 0; - PIN->vgridlist_geoid = nullptr; - PIN->vgridlist_geoid_count = 0; - /* Set datum parameters. Similarly to +init parameters we want to expand */ /* +datum parameters as late as possible when dealing with pipelines. */ /* otherwise only the first occurrence of +datum will be expanded and that */ diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake index 1b0db531..1cc3b90c 100644 --- a/src/lib_proj.cmake +++ b/src/lib_proj.cmake @@ -284,6 +284,8 @@ set(SRC_LIBPROJ_CORE proj_json_streaming_writer.hpp proj_json_streaming_writer.cpp tracing.cpp + grids.hpp + grids.cpp ${CMAKE_CURRENT_BINARY_DIR}/proj_config.h ) diff --git a/src/malloc.cpp b/src/malloc.cpp index d2e2d87a..050721e6 100644 --- a/src/malloc.cpp +++ b/src/malloc.cpp @@ -227,7 +227,6 @@ PJ *pj_default_destructor (PJ *P, int errlev) { /* Destructor */ /* free grid lists */ pj_dealloc( P->gridlist ); - pj_dealloc( P->vgridlist_geoid ); /* We used to call pj_dalloc( P->catalog ), but this will leak */ /* memory. The safe way to clear catalog and grid is to call */ diff --git a/src/proj_internal.h b/src/proj_internal.h index ffa533e7..1e43e066 100644 --- a/src/proj_internal.h +++ b/src/proj_internal.h @@ -65,6 +65,7 @@ #include <vector> #include "proj.h" +#include "grids.hpp" #ifdef PROJ_API_H #error proj_internal.h must be included before proj_api.h @@ -482,8 +483,7 @@ struct PJconsts { int gridlist_count = 0; int has_geoid_vgrids = 0; /* TODO: Description needed */ - struct _pj_gi **vgridlist_geoid = nullptr; /* TODO: Description needed */ - int vgridlist_geoid_count = 0; + std::vector<std::unique_ptr<NS_PROJ::VerticalShiftGridSet>> vgrids{}; double from_greenwich = 0.0; /* prime meridian offset (in radians) */ double long_wrap_center = 0.0; /* 0.0 for -180 to 180, actually in radians*/ @@ -839,12 +839,6 @@ void nad_free(struct CTABLE *); /* higher level handling of datum grid shift files */ -int pj_apply_vgridshift( PJ *defn, const char *listname, - PJ_GRIDINFO ***gridlist_p, - int *gridlist_count_p, - int inverse, - long point_count, int point_offset, - double *x, double *y, double *z ); int pj_apply_gridshift_2( PJ *defn, int inverse, long point_count, int point_offset, double *x, double *y, double *z ); 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; diff --git a/src/transformations/vgridshift.cpp b/src/transformations/vgridshift.cpp index de0cdd8c..51f129ab 100644 --- a/src/transformations/vgridshift.cpp +++ b/src/transformations/vgridshift.cpp @@ -22,7 +22,7 @@ static PJ_XYZ forward_3d(PJ_LPZ lpz, PJ *P) { PJ_COORD point = {{0,0,0,0}}; point.lpz = lpz; - if (P->vgridlist_geoid != nullptr) { + if (!P->vgrids.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, point.lp, Q->forward_multiplier); @@ -37,7 +37,7 @@ static PJ_LPZ reverse_3d(PJ_XYZ xyz, PJ *P) { PJ_COORD point = {{0,0,0,0}}; point.xyz = xyz; - if (P->vgridlist_geoid != nullptr) { + if (!P->vgrids.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, point.lp, Q->forward_multiplier); diff --git a/test/gie/4D-API_cs2cs-style.gie b/test/gie/4D-API_cs2cs-style.gie index e5722b5e..8d541823 100644 --- a/test/gie/4D-API_cs2cs-style.gie +++ b/test/gie/4D-API_cs2cs-style.gie @@ -448,6 +448,20 @@ expect 4.05 52.1 -10 ------------------------------------------------------------------------------- ------------------------------------------------------------------------------- +Test null grid with vgridshift +------------------------------------------------------------------------------- +operation proj=vgridshift grids=tests/test_nodata.gtx,null ellps=GRS80 +------------------------------------------------------------------------------- +ignore pjd_err_failed_to_load_grid +accept 4.05 52.1 0 +expect 4.05 52.1 -10 + +# Outside validity area of test_nodata.gtx. Fallback on null +accept 4.05 -52.1 0 +expect 4.05 -52.1 0 +------------------------------------------------------------------------------- + +------------------------------------------------------------------------------- Test bug fix of https://github.com/OSGeo/proj.4/issues/1025. Using geocent in the new API with a custom ellipsoid should return coordinates that correspond to that particular ellipsoid and not WGS84 as demonstrated in |