Bug #921: grdcut assigns zeros to nodes, instead of NaNs - GMT - GMT — The Generic Mapping Tools

Bug #921

grdcut assigns zeros to nodes, instead of NaNs

Added by Priyesh about 1 year ago. Updated 12 months ago.

Status:

Closed

Start date:

2016-07-04

Priority:

Normal

Due date:

Assignee:

Paul

% Done:

100%

Category:

Target version:

Candidate for next bugfix release

Affected version:

5.x-svn

Platform:

Description

grdcut assigns zeros to nodes, instead of NaNs when it is outside the master grid in y direction.
Eg :
gmt grdcut -R89/92/-1/2 fag.grd -Gokfine.grd -N -V
gmt grdcut -R91/93/3/5 fag.grd -GNoNan.grd -N -V

Sample grid (fag.grd) file attached

fag.grd - sample grid for test (191 KB) Priyesh, 2016-07-04 07:45

Associated revisions

Revision 16708
Added by Paul about 1 year ago

Address issue #921

History

#1 Updated by Paul about 1 year ago

Status changed from New to Resolved
Assignee set to Paul
Target version set to Candidate for next bugfix release
% Done changed from 0 to 100

Computed bad start column when should not do so. Fixed in r16708 I believe.

#2 Updated by Joaquim about 1 year ago

Yes, fixed but running fro MEX still fails. Thought it was the same issue but apparently not.

#3 Updated by Paul about 1 year ago

Not sure about that. On command line (after the fix) I now get via grdinfo:
NoNan.grd: 7080 nodes (49.2%) set to NaN
In mex, doing the same steps I get (apart from a non-fatal message I need to check)

q = gmt('grdcut', '-R91/93/3/5 fag.grd -N -V')   
g = gmt ('grdinfo', '-M', q)
[Session Matlab (0)]: Error returned from GMT API: GMT_PTR_NOT_NULL (65)
g.text(10) 

ans = 

    ': 7080 nodes (49.2%) set to NaN'

So looks like it is working fine. Please double check.

#4 Updated by Paul about 1 year ago

Fixed the annoying fake error message as well.

#5 Updated by Joaquim about 1 year ago

No, the error is while cropping from memory.

fag = gmt('read -Tg fag.grd');
check = gmt('grdcut -R92/94/2/4', fag);
gmt('grdinfo -M', check)

': Title: Produced by grdcu'
    ': Command: grdcut -R92/94/2/4 GMTAPI-000008 -G@GMTAPI@-00000'
    ': Remark: '
    ': Pixel node registration used [Geographic grid]'
    ': Unrecognized grid file format! Probably not a GMT grid'
    ': x_min: 90 x_max: 94 x_inc: 0.0166666666667 name: longitude [degrees_east n_columns: 240'
    ': y_min: 0 y_max: 4.01338 y_inc: 0.0166530290456 name: latitude [degrees_north n_rows: 241'
    ': z_min: -1308.73657227 at x = 92.2416666667 y = 3.98840045643 z_max: 1333.40112305 at x = 90.24…'
    ': scale_factor: 1 add_offset: 0'
    'GEOGCS["unnamed ellipse",…'

#6 Updated by Paul about 1 year ago

Status changed from Resolved to In Progress

OK, I see. We set the pad accordingly but at the "writing" back to Mex the pad is brutally overridden by the default pad so all the work in grdcut is lost. Will seek a solution.

#7 Updated by Paul about 1 year ago

Status changed from In Progress to Resolved

Fixed in r16715. Works for me now with gmtmex r540.

>> fag = gmt('read -Tg fag.grd')         

fag = 

                       z: [241x240 single]
                       x: [1x240 double]
                       y: [1x241 double]
                   range: [90 94 0 4.0134 -63.9878 79.2086]
                     inc: [0.0167 0.0167]
            registration: 1
                  nodata: NaN
    projection_ref_proj4: ''
      projection_ref_wkt: 'GEOGCS["unnamed ellipse",...'
                   title: 'Produced by grdproject'
                  remark: ''
                 command: 'grdproject -R90/94/0/4.01338 -Jm1i -I @GMTAPI@-000000 -Gfag.grd --PROJ_ELLIPSOID=Spher...'
                datatype: 'float32'
                  x_unit: 'longitude [degrees_east]'
                  y_unit: 'latitude [degrees_north]'
                  z_unit: 'meters, mGal, Eotvos, micro-radians or Myr, depending on img file and -S.'

>> check = gmt('grdcut -R92/94/2/4', fag)
grdcut: Warning: (s - y_min) must equal (NY + eps) * y_inc), where NY is an integer and |eps| <= 0.0001.
grdcut: Warning: s reset from 2 to 1.99836348548
grdcut: Warning: (n - y_min) must equal (NY + eps) * y_inc), where NY is an integer and |eps| <= 0.0001.
grdcut: Warning: n reset from 4 to 3.99672697095

check = 

                       z: [120x120 single]
                       x: [1x120 double]
                       y: [1x120 double]
                   range: [92 94 1.9984 3.9967 -63.9878 79.2086]
                     inc: [0.0167 0.0167]
            registration: 1
                  nodata: 0
    projection_ref_proj4: ''
      projection_ref_wkt: 'GEOGCS["unnamed ellipse",...'
                   title: 'Produced by grdcut'
                  remark: ''
                 command: 'grdcut -R92/94/2/4 @GMTAPI@-000004 -G@GMTAPI@-000005'
                datatype: 'float32'
                  x_unit: 'longitude [degrees_east]'
                  y_unit: 'latitude [degrees_north]'
                  z_unit: 'meters, mGal, Eotvos, micro-radians or Myr, depending on img file and -S'

#8 Updated by Joaquim 12 months ago

Status changed from Resolved to Closed

Also available in: Atom PDF

GMT

Issues