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.
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
ReplyDeleteThank you
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.
DeleteI 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 ?
DeleteI'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.
DeleteSorry, what we want to do is:
DeleteSpecify 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!
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.
DeleteHopefully, this is kind of what you were looking for?
Wow thanks for the quick answer!
ReplyDeleteI'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");
Hi, I need some help. I make this:
ReplyDelete'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...
Hello,
ReplyDeleteI have a problem, in the file.txt in the part of the result say "printing"
What does it means?
Thanks!
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.
DeletePrinting 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!
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
ReplyDeleteI'm having the same problem.
DeleteI 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
ReplyDeleteHi 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
ReplyDeletewrite only works when embedded in a script
Deletesame here. it always says unknown command
ReplyDeletewrite only works when embedded in a script, not in interactive command window
DeleteThanks 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.
ReplyDeleteI think one way you could do this is setting the extent you want in grads:
Delete'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!
Thanks Anonymous but in which line xmax and ymax values be replaced.
ReplyDeleteThis comment has been removed by the author.
ReplyDeleteThis comment has been removed by the author.
ReplyDeleteHi Admin
ReplyDeleteI 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
Can someone specify how to run this code in version 2.0.2? As others have pointed out, the `write` function is not available.
ReplyDeleteHi
ReplyDeleteThank 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.
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:
ReplyDeletetxt = '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.
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.
ReplyDeleteI 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