Earthdata Forum

Yes, the echo command needs to interpret '\n' as a linefeed character.
Another way to force this behavior is with the '-e' option.  (I have updated
my other post to show this.)  Another possibility would be to just use
convert to do the conversion to PGM (or TIFF or other 16-bit format).
From my example:
convert -size 8862x8730 -depth 16 GRAY:nak.dat $i.16bit.log.pgm

Note that you must use a projection suitable to your  Laguna de Términos
region.  In my example, I used an Albers equal area conic projection for
the South Atlantic.  If you just copy those projection parameters, none of
the level-2 pixels from your input file will map to that region.

Norman

Hi Norman,

Not sure if doing the right thing, but I'm still gettin blak images, here are my steps:

perl -e 'open X,"lon" or die;open Y,"lat" or die;open Z,"'rhos_469_q'" or die;while(read(Z,$z,4)){read(X,$x,4);read(Y,$y,4);print $x,$y,$z}' > rhos_469_q.bin
gmt mapproject rhos_469_q.bin -bi3f -bo3f -Jl-201/12/17.5/29.5/1:800000 -R-99.59/-79.1928/1.469/21.8741 > rhos_469_q_conic.bin
gmt nearneighbor rhos_469_q_conic.bin -Grhos_469_q_conic.grd -I0.01 -N4 -R-99.59/-79.1928/1.469/21.8741 -S0.04 -r -bi3f -V
gmt grdmath rhos_469_q_conic.grd 0.008 MAX 0.008 DIV LOG 2 0.008 DIV LOG DIV 1 MIN 65535 MUL = nak.grd
gmt grd2xyz nak.grd -ZTLHw > nak.dat
echo $'P5\n7732 7733\n65535'|cat - nak.dat>rhos_469_q.16bit.log.pgm

I don't want to reproject (if not neccesary) but to leave on WGS84 with 250 m resolution like original data (something like 0.002599 degress p/pixel).

There are a few problems here.  You are specifying a Lambert conformal conic
projection centered at (201 W, 12 N) with standard parallels at 17.5. and 29.5 N.
With that projection and the boundaries you give in the -R option, you get a map
like the one I attach below (only 20 times larger).  I'm guessing that this is not what
you had in mind.

The next problem is that the -R option you pass to nearneighbor is identical to
the one you pass to mapproject.  Since nearneighbor is operating on projected
Cartesian coordinates, you need to give it boundaries in that same space.

Try using mapproject with the -W option to get the boundaries you need.
For example, to get the width and height of the map I included below, I did this.

gmt mapproject -W -Jl-201/12/17.5/29.5/1:16000000 -R-99.59/-79.1928/1.469/21.8741
7.91705791492  7.91804080942

Note that I changed your map scale to get a much smaller map.  With the above
mapproject command, you would pass the following -R option to near neighbor.

-R0/7.92/0/7.92

In my examples, I also assume that PROJ_LENGTH_UNIT is set to inch.

See the manual pages for GMT for more information.
https://gmt.soest.hawaii.edu/doc/5.4.2/

Regards,
Norman

P.S.  It is necessary to reproject if you wish to eliminate the bow-tie effect.

Hi Norman, thank you for your help it was hard to understand how GMT works. Here my preliminary results:



I will share a complete script asap.

I'm glad that that is working for you.

Note that GMT has its own forum where you can ask questions
about topics that you cannot find answers to on their web site.

https://gmt.soest.hawaii.edu/projects/gmt/boards/1

Regards,
Norman