Replacing only negative Z values in GeoTIFF file, grdmath?

Q&A
1

Dear GMT masters,

I am trying to modify only negative Z values in a GeoTIFF file but I am still scratching my head about how to achieve this. My first thought was to use the neat grdmath, but after checking the operators I found the IFELSE and I think I would need to use an IF operator instead. Is there another GMT tool that may solve this issue? I also considered to use awk to do the job:

grd2xyz test.tif | awk '{for (i=1; i<=NF; i++) if ($3 < 0) $3 = ($3-180)*-1; print}' | tail -200 | head -10

so, re-gridding the result but I would prefer to avoid this way because this will be a time consuming process considering the amount of TIF files I have to analyze. I have uploaded a test file here as a MWE test.tif (314.5 KB), so I will be very much appreciated in case you can take a look at it and suggest me some GMT tools and/or ideas.

Best regards,

Gery

2

Hello gery

Have you try with grdclip?

3

Hola Esteban!

Thanks for the suggestion, it could work but not sure how to change the values below 0 -Sb0/???.

Thanks again,

PS. I found that gdal_calc.py may also work but it applies the changes to all values (ie, long, lat, z):

gdal_calc.py -A test.tif --outfile=out.tif --calc="abs(A-180)*(A<0)"
4

Yes, exactly. Check this example from the docs:

To set all values ≤ 0 to NaN in the file errors.nc:

gmt grdclip errors.nc -Gpos_errors.nc -Sb0/NaN+e -V

5

Unfortunately no luck so far, it seems that operations cannot be included, eg. -Sb0/((z-180)*-1) so everything below 0 should be positive after adding -180 to the Z value, perhaps I am misunderstanding something here…

6

Ok, I see. Do you want to use a formula to change the values? Ok, I think that grdclip only accepts a constant value.

7

You could try this:

  1. Keep the negative values.
  2. Changes values
  3. Combine new grid with previos tif file.
gmt grdclip test.tif -Sa0/NaN -Gtest_negative.nc
gmt grdmath test_negative.nc 180 SUB -1 MUL = new.nc
gmt grdmath new.nc test.tif AND = new_test.tif

I got this image from your test file

newtest
newtest1933×1406 275 KB

1 Like
8

My code most likely could be written in one single command using grdmath LE instead of grdclip.

9

I tried grdmath out of curiosity and came up with the following code:

gmt grdmath test.tif 0 LT 180 test.tif SUB MUL test.tif 0 GE test.tif MUL ADD = out.tif

my conclusion is that your code remains very much readable for humans, not only for that RPN calculator of grdmath. A big advantage for so many of us!

Funny thing is that the topic starter has actually got the result from gdal_calc.py. It could happily end up there but gdal_calc.py has not been able to recognize coordinate system from test.tif and silently threw it away. The result was a Cartesian grid without any CRS definition and swapped coordinates (gmt reads tiffs using gdal and gdal obviously assumes latitude comes first in a raster without any CRS). gmt appeared to be more forgiving and thoroughly kept coordinate system definition from test.tif.

1 Like
10

Both solutions worked nicely! thanks Esteban82 and mkononets for the nice support.