diff options
| author | Kristian Evers <kristianevers@gmail.com> | 2018-05-19 15:21:03 +0200 |
|---|---|---|
| committer | GitHub <noreply@github.com> | 2018-05-19 15:21:03 +0200 |
| commit | d9a0cefbf6bf84257f7f8da2b5ab32fecc9bfef3 (patch) | |
| tree | 706ca897c6dd4616092f587a63bfebcfe98da8c8 /src | |
| parent | 905b15d4682876e306f0eb97014baa09eb1a3402 (diff) | |
| parent | e9e9e72c9d36eac00e65bb137242d039891c78f5 (diff) | |
| download | PROJ-d9a0cefbf6bf84257f7f8da2b5ab32fecc9bfef3.tar.gz PROJ-d9a0cefbf6bf84257f7f8da2b5ab32fecc9bfef3.zip | |
Merge pull request #1004 from rouault/fix_1002
Vertical grid shift: do not interpolate node values at nodata value (fixes #1002)
Diffstat (limited to 'src')
| -rw-r--r-- | src/pj_apply_vgridshift.c | 65 |
1 files changed, 50 insertions, 15 deletions
diff --git a/src/pj_apply_vgridshift.c b/src/pj_apply_vgridshift.c index 58605df8..ca932674 100644 --- a/src/pj_apply_vgridshift.c +++ b/src/pj_apply_vgridshift.c @@ -35,6 +35,15 @@ #include "proj_internal.h" #include "projects.h" +static int is_nodata(float value) +{ + /* 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 > 1000 || value < -1000 || value == -88.88880f; +} + static double read_vgrid_value( PJ *defn, LP input, int *gridlist_count_p, PJ_GRIDINFO **tables, struct CTABLE *ct) { int itable = 0; double value = HUGE_VAL; @@ -112,23 +121,49 @@ static double read_vgrid_value( PJ *defn, LP input, int *gridlist_count_p, PJ_GR grid_iy2 = ct->lim.phi - 1; cvs = (float *) ct->cvs; - value = cvs[grid_ix + grid_iy * ct->lim.lam] - * (1.0-grid_x) * (1.0-grid_y) - + cvs[grid_ix2 + grid_iy * ct->lim.lam] - * (grid_x) * (1.0-grid_y) - + cvs[grid_ix + grid_iy2 * ct->lim.lam] - * (1.0-grid_x) * (grid_y) - + cvs[grid_ix2 + grid_iy2 * ct->lim.lam] - * (grid_x) * (grid_y); + { + 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) ) + { + double weight = (1.0-grid_x) * (1.0-grid_y); + value += value_a * weight; + total_weight += weight; + n_weights ++; + } + if( !is_nodata(value_b) ) + { + double weight = (grid_x) * (1.0-grid_y); + value += value_b * weight; + total_weight += weight; + n_weights ++; + } + if( !is_nodata(value_c) ) + { + double weight = (1.0-grid_x) * (grid_y); + value += value_c * weight; + total_weight += weight; + n_weights ++; + } + if( !is_nodata(value_d) ) + { + 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; + } } - /* 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 */ - if( value > 1000 || value < -1000 || value == -88.88880f ) - value = HUGE_VAL; - return value; } |