Monday, May 6, 2013

Tutorial: How to save data into a .txt ascii file using GrADS; Example: Gridded Temperature Data

Note: This tutorial has been updated in October 2019 to reflect the fact that the GDS has been discontinued.  The general principles of this tutorial have not changed, but now you need to download a small data file to your local computer to follow the tutorial.
Since GrADS cannot easily perform high level data analysis techniques, (barring a few complicated workarounds) it is useful to save data so it can be used with more advanced software.  While it is often logical to skip GrADS altogether, it can occasionally be easier to use GrADS as an intermediary step since it's a very intuitive scripting language for gridded datasets.

This tutorial will teach you how to save data using GrADS using the "write" function by guiding you through an example that will save latitude/longitude and monthly average surface temperature data from NCEP reanalysis into a .txt file.

Okay, so as always we will need to start by opening at data file.  For this tutorial, I am using the monthly mean 2 meter air-temperature reanalysis product.  To follow this tutorial, you will need to download the "air.mon.mean.nc" file from the ESRL PSD NCEP 1 Reanalysis website.  Once you download the file, and save it to your data directory, open the file in GrADS using "sdfopen"

  'sdfopen /home/your_datapath/air.mon.mean.nc'

Again, you are welcome to query the file to look at some of the metadata, but it's not required that you do so.  It is however good practice to plot the data and make sure it looks resonable before proceeding.  This dataset spans a long range of time, and for this tutorial, I am going to randomly choose February 2000 as the time value.  To specify the date instead of the timestep when picking a time, use the 'set time' (as opposed to 'set t') command.  Note that the time format must be %hhZ%dd%MMM%yyyy as in the example  below.

     'set time 00Z01FEB2000'
     'set gxout shaded'
     'd air'

Mean 2m air temperature from Feb 2000

Now that we have an idea what the data looks like on a map, we can move on to saving it into a .txt file.  Saving data from GrADS into an outside file is actually pretty straightforward.  The hard part is figuring out how to organize each piece of data exactly how you want it.

The first thing to do is to determine what data we would like to save, and how to organize it.  I find it a good practice to include the x/y grid indices when saving data into a text file.  So, the final text output will have five data columns, the x/y grid point indicies, the lat/lon value at each grid point, and the actual temperature data.  It's also good practice to include a header line at the top of the text file to specific which column is which.  So, the first step is to open the file and write header.  This is done by using the "write" function in GrADS.

     write('NCEP.txt', 'X    Y    LON    LAT    TEMP')

Once this command is executed in GrADS the file NCEP.txt will be created in your folder with the text "X    Y    LON    LAT    TEMP" on top.

Now that the data file has be started, it's time to start filling it up with data.  To do this, you need to set the GrADS  output to write our data out to the screen instead of plot it on a map.  This is done by setting 'gxout' to print.  In addition to that, you need to specify how GrADS should format our printed data.  This is done with the 'prnopts' command using C formating syntax.

     'set gxout print'
     'set prnopts %6.2f 1 1'    ;*%6.2f: c format, 1: values to plot on each line, ;space between values

Now that GrADS is set to print to the screen, and the print formatting is set, you are now ready to save the data.  The way to do this is to loop through all the grid-points and save the data to the file one grid-point at a time.  In this example, the nested loop to do this will first run across the x (longitude) direction, then through y (latitude), saving the information into rows as we go.  To begin, we will first need to set limits on our domain so that the loops know when to stop.  Do this by the use of the 'query' function

     'q dims'
     xline=sublin(result,2)    ;* 2nd line
     yline=sublin(result,3)    ;* 3rd line

     xmax=subwrd(xline,13)    ;*13th word on xline
     ymax=subwrd(yline,13)    ;*13th word on yline

     say 'X grid-points: 'xmax
     say 'Y grid-points: 'ymax
   
This code will print out the following to the screen:

   X grid-points: 193
   Y grid-points: 94

Now that domain dimensions are known, you can initialize the loops.  There is very little code within the nested loop.  First the 'display' command (which because of the 'gxout' and 'prnopts' configuration described above) is used to print out the temperature data at each grid-point to the "result" variable.  To save the data to the file, the only thingsleft to do are a) pull the temperature from the results data, b) pull the lat/lon values using the 'q dims' command, and c) write the data to a new row in the NCEP.txt file.  The code is as follows:

     y=1
     while(y<=ymax)
        x=1
        while(x<=xmax)
           'set x 'x
           'set y 'y
           'd air'

     * NOTE: It may be useful to test this to find out where the data is contained with in the result
     * It just so happens that in this case, the data is the 1st word of the 2nd line, this is not always true

           tmp=sublin(result,2)
           tmp=subwrd(tmp,1)

    *Get Lat/Lon Data
           'q dims'
           lons=sublin(result,2)
           lats=sublin(result,3)
           lon=subwrd(lons,6)
           lat=subwrd(lats,6)

    *Save data to file

    *Note the "append", so to add to the file instead of overwriting it
           write('NCEP.txt', x'    'y'    'lon'    'lat'    'tmp, append)
           x=x+1
      endwhile
      y=y+1
   endwhile

That's actually all there is to it.  This script will take a few moments to loop through all of the grid-points, and as it does so you will likely see the words "unknown command 0" printed a bunch of times on your screen.  I have often found it helpful to include a "say" command inside the loop, so to know where I am in the saving process.  When this script finishes you will have a lot of data saved into a .txt file that can be read into any plotting or statistical program.  For example, to confirm GrADS saved the data correctly to text, I created the image below by loading NCEP.txt into a Python array and plotted it on a map using Matplotlib and Cartopy.

NCEP.txt data plotted using Cartopy

As always, I included a downloadable script for this tutorial: Download this script

Note: The downloadable script has not been updated since it was originally created, and will therefore not work when trying to access the GDS.


28 comments:

  1. Is it also possible to print the single value for specified lat/lon coordinates ? Example: what is the 2m temperature for the coordinates: 51.445449, 5.581055. And put the result in kelvin or celcius in a variable or text file

    Thank you

    ReplyDelete
    Replies
    1. Yes, you can do this using a very similar method to the one described here. Just don't use the loop. Use the 'set lat' and 'set lon' command to fix your desired latitude and longitude coordinates, then use the 'set gxout print' command, display your variable (you can easily convert temperature to Kelvin, or Celcius in GrADS). Then use the write() function to write out the data to a file.

      Delete
    2. I got the output on commandline, works fine.. but now tje final problem is how to get this in a variable in PHP so we can use the values in calculations... Did you know this ?

      Delete
    3. I'm not exactly sure what you are trying to do, are you trying to take your output from GrADS and run the ascii file into a php script? In that case it isn't too hard to load up.

      Delete
    4. Sorry, what we want to do is:
      Specify a location in our script to get the value for that grid point in a variable in php. So for example, if we want the capesfc for the coordinates 56.0 -5.4 we retrieve one value of the amount of cape for that location. as lets say $cape

      We want to use the $cape in other calculations so the cape of the location specified above should be connected to $cape.

      Is that possible ? Thanks!

      Delete
    5. Presuming, you are running GrADS from a php program (i'm not sure this is possible as I am not an expert in php), you can probably feed arguments to GrADS when you run the program: e.g., "run grads -bx script.gs lat lon" or whatever command runs a program from php. Then I recommend writing the CAPE variable out from GrADS into a .txt file. One GrADS has done its thing, you can then use php to open up the .txt file and read the cape variable as $cape. At that point you should be able to do with it as you wish using php.

      Hopefully, this is kind of what you were looking for?

      Delete
  2. Wow thanks for the quick answer!

    I'm use grads with php, do you know how to use this in php?

    I have this so far:

    fwrite('NCEP.txt', 'X Y LON LAT TEMP');
    $ga->set("gxout","print");
    $ga->set("prnopts"," %6.2f 1 1");

    ReplyDelete
  3. Hi, I need some help. I make this:

    'q pos'

    say result
    px=subwrd(result,3)
    py=subwrd(result,4)

    say 'posicao x ='px
    say 'posicao y ='py

    'q xy2w 'px' 'py''
    say result

    plon=subwrd(result,3)
    plat=subwrd(result,6)


    **********************************************************
    * Guardar puntos geograficos

    longitud='colonlat.csv'
    fmt= '%8.3'
    rc1=math_format(fmt,plon)
    rc2=math_format(fmt,plat)
    varf=rc1","rc2
    rec= write(longitud,varf,closed)

    I get the file but I cannot get any value, so I dont now what to do...

    ReplyDelete
  4. Hello,

    I have a problem, in the file.txt in the part of the result say "printing"

    What does it means?

    Thanks!

    ReplyDelete
    Replies
    1. Hi, sorry I didn't get back to you sooner. This is a common problem, and it is pretty simple to fix! Okay, so when you set your graphics to print ('set gxout print') and you display your variable, you get output to the screen. What I did to get the variable "tmp" to write to the .txt file, I needed to grab the part of that output that is the actual value. To do this, I use the "sublin" and "subwrd" commands. Example below.

      Printing Grid -- 1 Values -- Undef = -9.99e+08
      278.05

      tmp=sublin(result,2)
      tmp=subwrd(tmp,1)

      In this example, the temperature is on the 2nd line, so the 2nd argument in the sublin function is 2, if it were on the 3rd line, it would be 3, and so on. Since it is the 1st word on the 2nd line, my argument for the subwrd function is 1.

      It looks like you are grabbing wrong line in your sublin command: which is why you get "printing" as a value. See what line your data is on, and change the sublin command to match.

      Good luck, hope this helps!

      Delete
  5. I was able to output an image and values using the script though reading a different data set from NCEP sst daily mean values. I just want to ask if there was an error somewhere when in my terminal screen it says there "unknown command : 0. How about if i want to output using a specified domain. How can i get through it...Thanks

    ReplyDelete
  6. I was able to output an image and values using the script though reading a different data set from NCEP sst daily mean values. I just want to ask if there was an error somewhere when in my terminal screen it says there "unknown command : 0. How about if i want to output using a specified domain. How can i get through it...Thanks

    ReplyDelete
  7. Hi I was wondering what version of Grads do I need to enable the write function. When I try running a script, a get an error message about the write function

    ReplyDelete
    Replies
    1. write only works when embedded in a script

      Delete
  8. same here. it always says unknown command

    ReplyDelete
    Replies
    1. write only works when embedded in a script, not in interactive command window

      Delete
  9. Thanks for this script.But i wonde how can i define a extent and extract the value for defines extent only. Because i have a dat for entire world when executing this script it is extracting the value for whoole world.

    ReplyDelete
    Replies
    1. I think one way you could do this is setting the extent you want in grads:

      'set lat 'lat1' 'lat2 ... etc

      Then doing a 'q dims', and I think one of the lines has the x/y dimensions for your subsetted domain. Then you can just replace the xmax and ymax values in the example provided here. That should get you started!

      Delete
  10. Thanks Anonymous but in which line xmax and ymax values be replaced.

    ReplyDelete
  11. This comment has been removed by the author.

    ReplyDelete
  12. This comment has been removed by the author.

    ReplyDelete
  13. Hi Admin

    I need a help.
    I just need to print my lon, lat and var,
    But, I could not get my var in the printed file.

    How, I know my var line?
    I still confused with subline and subwrd

    Thank you very much


    ReplyDelete
  14. Can someone specify how to run this code in version 2.0.2? As others have pointed out, the `write` function is not available.

    ReplyDelete
  15. Hi

    Thank you for the tutorial. I'm very new to GrADS. I'm following your tutorial to learn how to extract ascii data from netcdf file. Unfortunately, I failed to display the file. When I used http://monsoondata.org:9090/dods/rean2d.info, GrADS indicate 'couldn't ingest SDF meta data. I'm using the 2.0.a7 version. Could you help please.

    ReplyDelete
  16. How about in reverse? I can read a txt file in GrADS, but it's useless to me in any calculation. My script to read is as follows:

    txt = 'KYMN_Soil_Sites_FCHV.txt'
    i = 1
    while i < 10000
    input = read(txt)
    chk = sublin(input,1)
    if chk = 0
    nfo = sublin(input,2)
    date = subwrd(nfo,1)
    time = subwrd(nfo,2)
    st2 = subwrd(nfo,3)
    sm2 = subwrd(nfo,4)
    say date
    say time
    say st2
    say sm2
    else
    break
    endif
    i = i + 1
    endwhile
    input = close(txt)

    But, what can I do to make it usable for calculations from another script. Those variables I name as I read the data are gone as soon as it finishes.

    ReplyDelete
  17. what i wanted is to extract spatial average and is of different times and what i got is temporal average but of different spaces. what can i do in this situation.

    ReplyDelete
  18. I want to write time series data of a specific variable for a single location and save it as a .txt file which I plan to transfer to a spreadsheet. But I've tried using the write() command and it does not work at all; GrADS simply returns with 'unknown command'. How do I solve this problem?

    ReplyDelete