Monday, May 6, 2013

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

Since GrADS cannot easily perform high level data analysis techniques, (barring a few complicated workarounds) it may be useful to save data so it can be used with more advanced software.  It might seem logical at times to simply skip GrADS in the all together, which is fine.  However, I often find it easier to use GrADS as an intermediary step as it allows me to use only the data I need, and arrange the data as I see fit.  This tutorial will teach you how to save data using GrADS 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 the data file.  For this tutorial, I will use the monthly mean 2d reanalysis product.

     'sdfopen http://monsoondata.org:9090/dods/rean2d.info'

Before we do anything else, lets see what our variable looks like.  Since this dataset spans a large range of time, you can pick a time to look at; I am going to choose February,2000 as my time value.  I do this by using the 'set time' (as opposed to 'set t') command.  Once the time is set, simply display your variable.  The resultant image is shown below.

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

Average 2 meter temperature for Feb 2000
Now, that we see what our 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 very easy, we make use of a couple of simple commands and functions to do this.  The hard part is figuring out how to organize each piece of data exactly how you want it.

The first thing we want to do is determine what data we would like to save, and how we would like to organize it.  I have found that when writing 2d data to a text file, I like to include the x/y grid location, the lat/lon value at each grid point, and the actual temperature data.  So, our data file will have five columns to contain these data.  I like to know which column is which, so the first thing we will do is write header to start off our data file.  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 we have our data file started, it's time to start filling it up with data.  To do this, the first thing that we need to do is tell GrADS to write our data out to the screen instead of plot it on a map.  This is done by setting the output to print.  In addition to that, we need to tell GrADS how to format our printed data.  This is done with the 'prnopts' command.

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

So, GrADS is set to print to the screen, and our printing options are set, we are ready to save the data.  The way we are going to this is one grid-point at a time using a nested loop, looping first though x, then through y grid points, saving the information into rows as we go.  To begin, we will first need to set limits on our domain.  We 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 brick of code should print out the following to the screen:

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

Now we know our domain dimensions.  We can set up our loop.  Again, we will use a simple nested loop configuration to go through the data.  All we do in the loop, is print out the temperature data at each grid-point, get the lat/lon values using the 'q dims' command, and 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 t2m'

     * 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.

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


24 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
  8. same here. it always says unknown command

    ReplyDelete
  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