# Can't get grdcontour to plot correlation contours

**URL:** <https://forum.generic-mapping-tools.org/t/cant-get-grdcontour-to-plot-correlation-contours/3041>\
**Category:** Q&A\
**Created:** [June 21, 2022, 5:29pm UTC](https://forum.generic-mapping-tools.org/t/cant-get-grdcontour-to-plot-correlation-contours/3041 "2022-06-21T17:29:12Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![smauget](https://avatars.discourse-cdn.com/v4/letter/s/edb3f5/32.png) [@smauget](https://forum.generic-mapping-tools.org/u/smauget)\
**Post date:** [June 21, 2022, 5:29pm UTC](https://forum.generic-mapping-tools.org/t/cant-get-grdcontour-to-plot-correlation-contours/3041/1 "2022-06-21T17:29:12Z")

</div>

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.

 ![Screen Shot 2022-06-21 at 11.46.21 AM](https://canada1.discourse-cdn.com/flex047/uploads/gmt/original/2X/0/043eb413a105557e27656972e57f36f408b4fe64.png)

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:

```auto
-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?

```auto
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

```

---

<div class="post-metadata">

**Author:** ![Joaquim](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/joaquim/32/16_2.png) [@Joaquim](https://forum.generic-mapping-tools.org/u/Joaquim)\
**Post date:** [June 21, 2022, 5:54pm UTC](https://forum.generic-mapping-tools.org/t/cant-get-grdcontour-to-plot-correlation-contours/3041/2 "2022-06-21T17:54:24Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![John\_Geodesy](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/john_geodesy/32/116_2.png) [@John\_Geodesy](https://forum.generic-mapping-tools.org/u/John_Geodesy)\
**Post date:** [June 21, 2022, 6:46pm UTC](https://forum.generic-mapping-tools.org/t/cant-get-grdcontour-to-plot-correlation-contours/3041/3 "2022-06-21T18:46:52Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![pwessel](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/pwessel/32/30_2.png) [@pwessel](https://forum.generic-mapping-tools.org/u/pwessel)\
**Post date:** [June 21, 2022, 7:11pm UTC](https://forum.generic-mapping-tools.org/t/cant-get-grdcontour-to-plot-correlation-contours/3041/4 "2022-06-21T19:11:39Z")

</div>

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

1. 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.
2. You are drawing the zero contour twice. Probably should use either -Lp or -Ln to avoid that.
3. 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.
4. Of course, the process that made this data should afford a bit more decimals to not be off by so much in the coordinates.

---

<div class="post-metadata">

**Author:** ![Joaquim](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/joaquim/32/16_2.png) [@Joaquim](https://forum.generic-mapping-tools.org/u/Joaquim)\
**Post date:** [June 21, 2022, 7:23pm UTC](https://forum.generic-mapping-tools.org/t/cant-get-grdcontour-to-plot-correlation-contours/3041/5 "2022-06-21T19:23:43Z")

</div>

> [@smauget](#):
>
> `gmt xyz2grd eof.xyz -Gtemp.grd -J`

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

---

<div class="post-metadata">

**Author:** ![pwessel](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/pwessel/32/30_2.png) [@pwessel](https://forum.generic-mapping-tools.org/u/pwessel)\
**Post date:** [June 21, 2022, 7:35pm UTC](https://forum.generic-mapping-tools.org/t/cant-get-grdcontour-to-plot-correlation-contours/3041/6 "2022-06-21T19:35:17Z")

</div>

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