Monday, May 13, 2013

Tutorial: Use the tcorr function to create a map of regression/correlation coefficients in GrADS

Making a correlation map in GrADS is actually quite simple.  A regression/correlation map, is basically a map of regression and correlation coefficients plotted against a time-series of some variable.  The classic example is the monthly average time-series of El-Nino vs. surface temperature on the globe.  The example used in this tutorial will be 500mb height correlated with the 500mb height over the northeast U.S.  For this tutorial we will use the NCEP reanalysis data.

GrADS makes correlation maps very simple.  All you need is a data set that varies in both space and time.  The first step is to define your time-series variable.  In the interest of time and simplicity, I will just take the 500mb height at one point rather than average over a small spatial area, as is often done when doing these maps.  So lets open the data file and set our time-series variable!

    'set t 1 636'
    'set lev 500'
    'set lat 44'
    'set lon 282'

The variable "ts_var" is our time-series variable, and basically we saved the time-series of geopotential height at 500mb at the lat/lon point of 44/282.  Now, we are going to use two intrinsic GrADS functions to plot the regression coefficients (shaded) and the correlation coefficients (contoured above 0.5).

So, to plot our coefficients, we use the 'tregr()' and the 'tcorr()' functions.  Both of these functions do essentially the same thing; you feed them a time-series variable, and a spatial variable, as well as time-constraints, and the function plots your coefficients.

Example:  tregr(x, y, t=1, t=50)

This variable takes the regression coefficients of spatial variable y and time-series x through the first 50 time-steps.  In our example of 500mb height our code looks like:
    'set gxout shaded'
    'set t 1'
    'set lat 0 90'
    'set lon 180 360'
    'd tregr(ts_var,z,t=1,t=636)'

    'set gxout contour'
    'set cstyle 2'
    'set cthick 6'
    'set ccolor 1'
    'set clevs 0.5 0.6 0.7 0.8 0.9 1.0'
    'd tcorr(ts_var,z,t=1,t=636)'

Most of the above code is for display of the contour variable; style, thickness, levels, but you can see the tregr() and the tcorr() functions used correctly.  If you are following along at home, the resultant map is plotted.

Shaded:Regression Coefficients, Contoured: Correlation Coefficients

So, what you can gather from this map is that the 500mb geopotential height is generally well correlated above 20 degrees N (r>0.8).  Furthermore, you can pick out the Rossby wave pattern in the regression coefficients.  In anycase, this example show give you an idea of how to make regression/correlation maps in GrADS.

Download Example Script


  1. Is there a way to create a correlation map between two spatially varying variables? Say I want to, for example, compare two Land Surface Models over a certain period, and I want to show a correlation map between the output from these models. It seems here that the first argument in tcorr must only vary in time but not in space.

    The correlation is therefore between the time series of each grid square among the two models and the result is a map of these correlations.

  2. Unfortunately, I do not know of a great way to do this in GrADS. As I understand your question, you basically want to grab the correlation coefficient between two time-varying variables at each lat/lon coordinate and then plot the resultant map of correlation coefficients on the map. Is this correct?

    If this is correct, the best way to do this with GrADS would be to grab the correlation coefficient at each grid point and save it to a .txt file along with the lat/lon coordinate and then plot the map in matlab, or idl or something.

    This option is very possible, I have another tutorial on how to save data to a .txt file:

    The only issue with doing it this way, is that I don't think that GrADS has an intrinsic function to get a simple correlation coefficient between to variables, so you would have to work around it and calculate it yourself. This is a little bit tricky, but basically, you would have to run a triple-nested loop in lon, lat, and then time. The time loop you would grab the value at each time-step and build two arrays, one for each model. The run the arrays through a function to calculate the correlation coefficient. You can follow the model used in my script to do this:

    But yeah, GrADS isn't great for advanced statistical analysis, and while there are work arounds to it, they can by somewhat complicated.

    Hope this helps!, Cheers!

  3. I'm probably in the wrong tutorial to ask this, but I didn't found a better place to post.
    Is there a way to change the size of my map's GRID ?
    Let me try to explain better: I have a template that reads and plots a map with 0.25 degrees of spacing, but I need to plot a map with 0.40 to compare with another map that I have, which is with 0.40 of spacing. I already try to use the lterp function, but didn't works... I probably made something wrong.

    My .ctl have this: XDEF 210 LINEAR -82.00 0.25
    But I need this: XDEF 126 LINEAR -82.00 0.40

    Thank you!

    1. The lterp function should work for you, simply take the hi-resolution file as the first argument in the function, and the course resolution file as the second.

      Example 4 from this website, Shows you exactly what you want to get.

      Only in your case, you would set x to 1 126 as that is the size of your course grid.

      Hope this helps!

    2. It's that one that I tryed, and I did what you told me to... I think I'm doing something wrong.

      I have a lot of .bin and .ctl files, so, I created a template to read them all, but the resolution is 0.25 and I need a 0.40 resolution to compare with another template.
      So, I tryed this lterp function, but I only have one template to put in the arguments, so, I created another template, just like the first one, but with 0.40 resolution to put in the argument of the function. But this template didn't work, so, the final result of my script with the lterp function didn't work as well.

      I know I'm probably doing something wrong, but I'm new in GrADS and I don't know exatcly what to do.

      I don't know if you understand me or if you can help me, but thanks anyway!

  4. please explain- 'ts_var=z'. It shows an error when I type it interactively. I understand using the set function you have created a time series. But this 'ts_var=z" isnt clear. Thanks

  5. dear sir please how to write script in between geopotential hgight and slp,at 500mb. At every grid point correlation how to write.
    pls send me sir,
    thank you sir