Can't get grdcontour to plot correlation contours

smauget · June 21, 2022, 5:29pm

I’m trying to use GMT ver 6.3.0 on an iMac (12.3.1) to plot contours of correlations between a principal component (PC) time series and the data time series used to generate the PC. The result should be a map of correlations over a tropical pacific region from 90W to 230 W. and ~ -10 S to 10 N.

The lon - lat - correlation values on a 1.0 deg lon X 0.3333 deg lat grid are in a text file called eof.xyz:

-229.50   -9.83     0.8244
-228.50   -9.83     0.7973
-227.50   -9.83     0.7303
-226.50   -9.83     0.6468
-225.50   -9.83 -9999.0000
-224.50   -9.83 -9999.0000
-223.50   -9.83 -9999.0000
-222.50   -9.83 -9999.0000
-221.50   -9.83 -9999.0000
.........................
.........................
 -98.50   10.17    -0.0890
 -97.50   10.17    -0.0969
 -96.50   10.17    -0.1081
 -95.50   10.17    -0.1353
 -94.50   10.17    -0.1644
 -93.50   10.17    -0.2029
 -92.50   10.17    -0.2444
 -91.50   10.17    -0.3059
 -90.50   10.17    -0.3697
#
#
# In a c-shell script I first plot the  equatorial Pacific map...  
gmt pscoast -R-230.0/-90.0/-10.0000/10.333333  -Ba15/a10WSne  -Dc -P  -G120/120/120 -K  -X1.5 -Y8.0   -JQ-160.0/6.0 >! map.ps
#
# ... then convert eof.xyz to temp.grd
gmt xyz2grd eof.xyz -Gtemp.grd -J -I1.0/0.333333333 -Ddegree/degree/sshg/=/=/=/= -R  -V
#
# Next, all null values (-9999.0000 ) in temp.grd are reassigned  as NaN in a temp0.grd file
gmt grdclip temp.grd -Gtemp0.grd -Sa1.0/NaN -Sb-1.0/NaN -V
#
# A first grdcontour call contours all positive correlation values @ 0.10 increments in red
gmt grdcontour  temp0.grd -R -B -J -P  -A0.10    -LP  -W0.75p,red  -K -O  -V     >> map.ps
#
# A second call contours all negative correlation values @ 0.10 increments in blue
gmt grdcontour  temp0.grd -R -B -J -P  -A0.10    -LN  -W0.75p,blue    -O  -V     >> map.ps

Problem is, no contours are plotted. When I try to dump the contours to a txt file using -Dtemplate.txt, nothing ends up in template.txt. I don’t see anything in the following verbose output of the xyz2grd-grdclip-grdcontour calls that looks out of place. Both grdcontour calls say they are tracing the contours but no contours are generated in template.txt or plotted in map.ps. Any idea what I might be doing wrong?

xyz2grd [INFORMATION]: Processing input table data
xyz2grd [INFORMATION]: Cartesian input grid
xyz2grd [INFORMATION]: Cartesian input grid
xyz2grd [INFORMATION]: n_columns = 141  n_rows = 62  nm = 8742  size = 8742
xyz2grd [INFORMATION]: Reading Data Table from file eof.xyz
xyz2grd [INFORMATION]: Data records read: 8540  used: 8540  nodes filled: 2911 nodes empty: 5831 [set to NaN]
xyz2grd [INFORMATION]: Extend grid via copy onto larger memory-aligned grid
xyz2grd [INFORMATION]: Writing grid to file temp.grd
xyz2grd [INFORMATION]: Proj4 string to be converted to WKT:
	+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=-160 +x_0=0 +y_0=0 +units=m +a=6378137.000 +b=6356752.314 +ellps=WGS84 +datum=WGS84 +units=m +no_defs
xyz2grd [INFORMATION]: WKT converted from proj4 string:
PROJCS["unknown",
    GEOGCS["unknown",
        DATUM["unknown",
            SPHEROID["unknown",6378137,298.257220143428]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]]],
    PROJECTION["Equirectangular"],
    PARAMETER["standard_parallel_1",0],
    PARAMETER["central_meridian",-160],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH]]
grdclip [INFORMATION]: Processing input grid
grdclip [INFORMATION]: Round-off patrol changed grid limit for ymin from -9.999999836065573 to -9.999999836065573
grdclip [INFORMATION]: Round-off patrol changed grid limit for ymax from 10.33333316393443 to 10.33333316393443
grdclip [INFORMATION]: Cartesian input grid
grdclip [INFORMATION]: Reading grid from file temp.grd
grdclip [INFORMATION]: Cartesian input grid
grdclip [INFORMATION]: gmt_grd_BC_set: Set boundary condition for all edges: natural
grdclip [INFORMATION]: gmt_grd_BC_set: Set boundary condition for left   edge: natural
grdclip [INFORMATION]: gmt_grd_BC_set: Set boundary condition for right  edge: natural
grdclip [INFORMATION]: gmt_grd_BC_set: Set boundary condition for bottom edge: natural
grdclip [INFORMATION]: gmt_grd_BC_set: Set boundary condition for top    edge: natural
grdclip [INFORMATION]: Writing grid to file temp0.grd
grdclip [INFORMATION]: 119 values < -1 set to nan
grdclip [INFORMATION]: 0 values grdcontour [INFORMATION]: Processing input grid
grdcontour [INFORMATION]: Spherical approximation used
grdcontour [INFORMATION]: Map scale is 2594.55 km per inch or 1:1.02148e+08.
grdcontour [INFORMATION]: Reading grid from file temp0.grd
grdcontour [INFORMATION]: gmt_grd_BC_set: Set boundary condition for all edges: natural
grdcontour [INFORMATION]: gmt_grd_BC_set: Set boundary condition for left   edge: natural
grdcontour [INFORMATION]: gmt_grd_BC_set: Set boundary condition for right  edge: natural
grdcontour [INFORMATION]: gmt_grd_BC_set: Set boundary condition for bottom edge: natural
grdcontour [INFORMATION]: gmt_grd_BC_set: Set boundary condition for top    edge: natural
grdcontour [INFORMATION]: Tracing the 0 contour
grdcontour [INFORMATION]: Tracing the 0.1 contour
grdcontour [INFORMATION]: Tracing the 0.2 contour
grdcontour [INFORMATION]: Tracing the 0.3 contour
grdcontour [INFORMATION]: Tracing the 0.4 contour
grdcontour [INFORMATION]: Tracing the 0.5 contour
grdcontour [INFORMATION]: Tracing the 0.6 contour
grdcontour [INFORMATION]: Tracing the 0.7 contour
grdcontour [INFORMATION]: Tracing the 0.8 contour
grdcontour [INFORMATION]: Processing input grid
grdcontour [INFORMATION]: Spherical approximation used
grdcontour [INFORMATION]: Map scale is 2594.55 km per inch or 1:1.02148e+08.
grdcontour [INFORMATION]: Reading grid from file temp0.grd
grdcontour [INFORMATION]: gmt_grd_BC_set: Set boundary condition for all edges: natural
grdcontour [INFORMATION]: gmt_grd_BC_set: Set boundary condition for left   edge: natural
grdcontour [INFORMATION]: gmt_grd_BC_set: Set boundary condition for right  edge: natural
grdcontour [INFORMATION]: gmt_grd_BC_set: Set boundary condition for bottom edge: natural
grdcontour [INFORMATION]: gmt_grd_BC_set: Set boundary condition for top    edge: natural
grdcontour [INFORMATION]: Tracing the -0.8 contour
grdcontour [INFORMATION]: Tracing the -0.7 contour
grdcontour [INFORMATION]: Tracing the -0.6 contour
grdcontour [INFORMATION]: Tracing the -0.5 contour
grdcontour [INFORMATION]: Tracing the -0.4 contour
grdcontour [INFORMATION]: Tracing the -0.3 contour
grdcontour [INFORMATION]: Tracing the -0.2 contour
grdcontour [INFORMATION]: Tracing the -0.1 contour
grdcontour [INFORMATION]: Tracing the 0 contour

Joaquim · June 21, 2022, 5:54pm

Please wrap your script/numbers under triple back ticks as I just did.

John_Geodesy · June 21, 2022, 6:46pm

As a start, did you try running a grdinfo -M on temp.grd and temp0.grd? This will help to make sure that the .grd files look as you expect them to look.

Also, if the .grd files actually were made, then you could use grd2xyz just to double check the values and the coordinates. Again, this is more of a check.

I’m not that familiar with the option -LP in grdcontour. Especially without specifying a low/high value, which I don’t think you can do. That might be something else to check. I’m not sure that you need that -L option at all.

pwessel · June 21, 2022, 7:11pm

Hard to tell without seeing the grid, but I spot a few problem areas:

Looking at your input data, I see latitudes such as 10.17. You say the interval in y is 1/3 degree (20m) but 10.17 is not a multiple of that. It is 10 degrees 10.2 minutes. 10 minute is not a multiple of 20 from -10 degrees so you probably have a pixel-registered file and you need to add -r to your xyz2grd command. Otherwise I expect every other node is NaN and no contours can be drawn.
You are drawing the zero contour twice. Probably should use either -Lp or -Ln to avoid that.
0.33333 etc is not very good, I suggest you use 20m as the interval and 10:20 as the north border to avoid roundoff.
Of course, the process that made this data should afford a bit more decimals to not be off by so much in the coordinates.

Joaquim · June 21, 2022, 7:23pm

I suspect the problem is that -J that inherits from history and creates one grid saying it’s projected. Don’t use that -J

pwessel · June 21, 2022, 7:35pm

Probably adds some bad metadata, but GMT is not actively using those things. I think my point 1 is the issue.