Hi All,
I am trying to take some ASCII data that is already on a regular grid and make it into a .nc grid file using xyz2grd so I can plot it with grdimage. The script below produces a plot, but if you look closely, the pixels colored by grdimage don’t line up with the source data points plotted with psxy. I’m feeling like I must be missing something. Any ideas?
The code that makes the plot is below and the source csv file is linked here. csv file (327 KB)
#!/bin/bash
#filenames
plotFile="plotPixels.eps"
#set GMT defaults
gmt set PS_MEDIA=8.5ix8.8i
gmt set MAP_FRAME_TYPE=plain
#get range of data
R=$(gmt info test.csv -Ie)
#convert the data to a nc grid file
gmt xyz2grd test.csv $R -Gtest.nc -I0.0492
#plot the data as a raster image
gmt grdimage test.nc -X0.6i -Y0.5i -R -JM7.5i -Cjet -P -K > $plotFile
#plot the source data points
gmt psxy test.csv -R -JM -Sc1.5p -Gblack -O -K >> $plotFile
#plot the axes
gmt psbasemap -R -JM -Bxa1f0.5 -Bya1f0.5 -BWeSn -O >> $plotFile
#open the plot
evince $plotFile &
The gmt info command above did return
-R-124.7187/-119.4051/35.6279/40.2035
From reading the manual, I suspected I needed pixel registration, but when I add a -rp flag to xyz2grd, I get the warning message below and an even stranger result shown below.
xyz2grd [WARNING]: 94 values gave bad indices: Pixel vs Gridline registration confusion?
It seems like there is a registration issue and/or my data isn’t evenly gridded, but when I tested the data, it appeared to be correctly gridded. So, I think I am still missing something. Could this be a bug? Seems like I must be doing something wrong. When I grid a subset of the data, things look fine.
Let me know if you need more info, and thanks for the help!
-Scott
I also played with -rp and -rg and nothing worked. Everything lines up at the top and bottom of the plot, but things get little by little off centered in the middle, and -rp when used still reports bad indices. Everything in the x-direction (longitude) is as expected, the issues are only in the y-direction (latitude).
Could this be a bug? or some sort of unexpected floating point error?
this is expected and how plotting works. see more below.
grdimagealways plots on a regular cartesian grid (this is how it works internally), but Mercator projection -JM/m... is vertically irregular (irregular longitude intervals). The result is that the top and the bottom rows of cells align perfectly while grid cells in the middle appear slightly off.
to get your grid plot exactly aligned with the map frame grid you can do one of two things
use a vertically regular projection like -JQ/q... when plotting with grdimage
use grdview -T to plot your grid. This makes grid cells to align exactly with an irregular projection. See grdview — GMT 6.5.0 documentation
Using the -R obtained with gmtinfo is the same as saying use grid registration. Using that -R and ´-rp´ corresponds to say: put the cell centers 1/2 inc away from this and naturally some points fall off their bins. There no bug here.
I didn’t try with Mercator bur the linear case worked beautifully with the limits I gave you. Did you try it?
How did I get it? Well, GDAL interprets files with .xyz extension as gridded data as long as they have the right layout (xyz2grd is more tolerant in this regard). So I renamed your file to .xyz but had also to replace the NaNs by zero because GDAL can’t recognize them. Next I did
grdconvert file.xyz -Gfile.nc
and looked at the nc grid limits.
But @mkononets is right with his explanation on why using Mercator moves the points away inside pixels.
This was enlightening and very helpful! I have used GMT for years and didn’t know that grdimage plots on a regular cartesian grid. I assumed that the -J projection option would set whatever regular or irregular grid was needed. Interesting! So many little things now make sense.
I was able to get the pixels to line up with both of your suggestions. Below is the grdview -T method.
Thanks so much @mkononets and @Joaquim!
