Grdtrack gives negative values from an input grid with only positive values

Q&A
1

I’m using the earth_age grid, which contain positive numbers only.

For some reason, grdtrack gives me negative values.
Have anyone had any issues with this?

Example:

# create grid
$ gmt grdcut @earth_age_01m -R-7/2/71/75 -Gage_wgs84.nc

# check grid - looks good
$ gmt grdinfo age_wgs84.nc
age_wgs84.nc: Title: 
age_wgs84.nc: Command: gmt grdcut @earth_age_01m_g/ -R-7/2/71/75 -Gage_wgs84.nc
age_wgs84.nc: Remark: 
age_wgs84.nc: Gridline node registration used [Geographic grid]
age_wgs84.nc: Grid file format: nf = GMT netCDF format (32-bit float), CF-1.7
age_wgs84.nc: x_min: -7 x_max: 2 x_inc: 0.0166666666667 (1 min) name: longitude n_columns: 541
age_wgs84.nc: y_min: 71 y_max: 75 y_inc: 0.0166666666667 (1 min) name: latitude n_rows: 241
age_wgs84.nc: v_min: 0.0100021362305 v_max: 50.7999992371 name: z
age_wgs84.nc: scale_factor: 1 add_offset: 0
age_wgs84.nc: format: netCDF-4 chunk_size: 136,241 shuffle: on deflation_level: 3
age_wgs84.nc: Default CPT: @earth_age.cpt

# convert to utm31
$ gdalwarp -s_srs EPSG:4326 -t_srs EPSG:32631 age_wgs84.nc age_utm31.nc
Processing age_wgs84.nc [1/1] : 0Using internal nodata values (e.g. nan) for image age_wgs84.nc.
...10...20...30...40...50...60...70...80...90...100 - done.

# check again - looks good
$ gmt grdinfo age_utm31.nc
age_utm31.nc: Title: z
age_utm31.nc: Command: gmt grdcut @earth_age_01m_g/ -R-7/2/71/75 -Gage_wgs84.nc
age_utm31.nc: Remark: 
age_utm31.nc: Pixel node registration used [Cartesian grid]
age_utm31.nc: Grid file format: nf = GMT netCDF format (32-bit float), CF-1.7
age_utm31.nc: x_min: 137638.77365 x_max: 470972.377092 x_inc: 903.34309876 name: x coordinate of projection [m] n_columns: 369
age_utm31.nc: y_min: 7876431.83181 y_max: 8348880.27247 y_inc: 903.34309876 name: y coordinate of projection [m] n_rows: 523
age_utm31.nc: v_min: 0.0100021362305 v_max: 50.7999992371 name: z
age_utm31.nc: scale_factor: 1 add_offset: 0
age_utm31.nc: format: classic
age_utm31.nc: Default CPT: @earth_age.cpt
PROJCS["WGS 84 / UTM zone 31N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",3],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32631"]]

# negative values???
$ gmt grdtrack -Gage_utm31.nc 20.cells | grep - | head
421950	8010775	1	-0.0115028944039
421950	8010750	1	-0.00968893837643
421975	8010750	1	-0.00932320722593
421950	8010725	1	-0.00759451048117
421950	8010700	1	-0.00521364665683
421950	8010675	1	-0.00254038284217           <------ negative numbers
417275	8009275	1	-8.65589427942e-05
417275	8009250	1	-0.000713623031298
417275	8009225	1	-0.00115224396813
417300	8009225	1	-0.00075609916143
[...]

In total 971 values are negative.
There should be none!

2

I suspect that is due to the resampling of your grid. By default is bicubic interpolation. So try using the adding the -nn option in grdtrack.

1 Like
3

You seem to have fixed my problem.
Thanks @Esteban82!

I now only get positive values - as expected. I didn’t think this operation required resampling. From the manpage, -A option:

For track resampling (if -C or -E are set) we can select how this is to be performed.

No -C or -E used in my code.

4

You are welcome! By the way, I think a linear interpolation (-nl) could be better for your grid.

1 Like
5

Except if the points fall over the grid nodes, there is always the need for interpolation from nearby nodes and is doing cubic by default.

1 Like
6

That makes sense.

I thought that the points indeed fell over the grid nodes, but of course they won’t; I’m translating a 1’ grid to UTM, and give grdtrack input points derived from a 25mx25m grid with a “clean” region (multiples of 25m).

Thanks alot all.