[gradsusr] Displaying data from NetCDF; outputting text files; using comparitive statistics

Muhammad Rahiz muhammad.rahiz at ouce.ox.ac.uk
Thu Dec 16 09:33:23 EST 2010



-- 
Muhammad Rahiz
Researcher & DPhil Candidate (Climate Systems & Policy)
School of Geography & the Environment
University of Oxford

On Thu, 16 Dec 2010, James Ciarlo` wrote:

> Dear Muhammad,
> The problem with the nc file is now solved thank you. I will contact the RegCNET to understand
> better the process that produces the nc files so that I can understand the problem with that
> file.
> The processing that I was doing was from the ESRL web site, to obtain a select region and time
> frame. But that I can do with GrADS so it is not a problem.
> 
> Regarding the fprintf. I need to 
> sdfopen input.nc
> set lon x
> set lon y
> set lev z   #so one value for each, one point in the map
> set t 1 150 #so basically I am setting a time line
> d prmsl # this gives me a graph of pressure v time
> #Then from fprintf I need to get the values for both pressure and time
> I use 'frpintf prmsl data.txt' to get the pressure values
> But then I need to get the corresponding time values, even in a separate text document, it
> doesn't matter, as long as I can get it.

It seems that you want to extract the timeseries of a particular grid 
point. The attached script would do that (items marked with ** needs 
to be changed to suit your data).

What it does is to make a loop, and prints out two text files - one for 
data, and one for time stamp.

Once you get the two text files out, you can join them together using the 
command

  $ paste x.time.txt x.data.txt > x.output.txt

The timestamp that I have for my file is something like 12Z28JAN2061. You 
may want to remove unwanted characters like 12Z so that it becomes 
28JAN2061. To do so, run

  $ sed -i 's/12Z//g' x.time.txt

> DO you think it's possible?
> 
> James
> 
> On 15 December 2010 16:16, Muhammad Rahiz <muhammad.rahiz at ouce.ox.ac.uk> wrote:
> 
>
>       --
>       Muhammad Rahiz
>       Researcher & DPhil Candidate (Climate Systems & Policy)
>       School of Geography & the Environment
>       University of Oxford
>
>       On Wed, 15 Dec 2010, James Ciarlo` wrote:
>
>       Dear Muhammad,
>       A) I have sent you the nc file for file3.jpg, but I cannot send you the nc
>       file for file1.jpg
>       and file2.jpg as it is 35 GB large. Perhaps we can figure it out through
>       this one?
> 
> 
> The attached ESRL_SLP.nc file has undergone some processing. You can look at the history
> of commands on the file by
> 
>  $ ncdump -h ESRL_SLP.nc
> 
> I went the source of the data at
> www.esrl.noaa.gov/psd/data/gridded/data.20thC_ReanV2.monolevel.mm.html and downloaded
> the file prmsl.mon.mean.nc (it is the same file you have judging from the file history).
> I then made a test plot in GrADs (see prmsl.png). Plots fine.
> 
> As I said, some processing was done on the file you have. I don't know exactly what
> language (I suspect NCL) and what it's trying to do. Perhaps you can tell what is it
> that you're trying to do with the file, I could possibly help.
> 
>
>       B) I have set up the script fprintf.gs and have managed to get it
>       operational. It works well
>       and I thank you for it. I haven't been able to produce the text for time and
>       elevation though.
>       I tried the same procedure as you mentioned for lon.txt and lat.txt for
>       time, but it didn't
>       work? does it require another script or maybe a different way to define
>       time?
> 
> 
> I'm not sure what you mean by wanting to get the ascii grid for time. Each ascii that
> you generate corresponds to a particular time. Hence if you want to get ascii grids for
> all time steps, just do something like
> 
> 'sdfopen input.nc'
> 'set t 1'
> 'd prmsl'
> 'fprintf prmsl prmsl.t1.txt'
> 
> 'set t 2'
> 'd prmsl'
> 'fprintf prmsl prmsl.t2.txt'
> 
> Of course, you can do this in a loop for all time steps.
> 
> I'm not sure if elevation (topography) is a variable of climate data (I've not come
> across any). Do you mean pressure levels like 500mb, 100mb etc? If yes, then you can do;
> 
> 'set lev 500'
> 'set t 1'
> 'd prmsl'
> 'fprintf prmsl prmsl.time1.lev500.txt'
> 
> This gives you the data for 500mb pressure level at time, t=1.
> 
> 
>
>       Regards,
>
>       James
>
>       On 15 December 2010 10:46, Muhammad Rahiz <muhammad.rahiz at ouce.ox.ac.uk>
>       wrote:
>            Hi James,
>
>            I suspect a corrupted data file. Could you attach the original
>       datafile, if it's
>            not too big? And attach the commands that you use to plot as well.
>
>            --
>            Muhammad Rahiz
>            Researcher & DPhil Candidate (Climate Systems & Policy)
>            School of Geography & the Environment
>            University of Oxford
>
>       On Wed, 15 Dec 2010, James Ciarlo` wrote:
>
>       Dear Muhammad,
>       I am sorry it took me so long to answer. Thank you for your
>       answer on point 2. 
>
>       Regarding point 1, opening .nc files. I have attached some
>       files. 
>       The images from file1.jpg and file2.jpg are grads images
>       (printim) from nc files produced with RegCM4, one is gxout
>       contour and the other is gxout shaded.
>       The image from file3.jpg is a grads image (printim) gxout shaded
>       for an nc file from the ESRL database. 
>       The last image was produced with the PANOPLY data viewer from
>       the same file that produced file3.jpg
>
>       As you can see the data is not being plotted with grads as with
>       normal .ctl files (which are working properly).
>       I cannot seem to understand why it is not reading it properly.
>       Maybe there are some definitions that I have to set? How should
>       I check if all the settings are correct?
>
>       Regarding point 3, I think I can do them as a normal
>       mathematical expression, but I'll be trying that later on. 
>
>       Regards,
>
>       James
>
>       On 11 December 2010 12:24, James Ciarlo`
>       <james.ciarlo at physics.org> wrote:
>            Dear Muhammad,
>       I just had a look at it this morning. Unfortunately I
>       haven't yet had the time to put it to use. Thank you very
>       much for your help. I will get back to you about it for
>       sure. Thanks again
>
>       Regards,
>
>       James
> 
>
>       On 11 December 2010 09:59, Muhammad Rahiz
>       <muhammad.rahiz at ouce.ox.ac.uk> wrote:
>            Hi James,
>            I was wondering if you got my reply to the
>            query you posted because I'm not sure if it
>            gets delivered successfully.
>
>            Anyway, if you did not, these are what I
>            suggested
>
>            1. display error
>            - Try set gxout to something other than the
>            default e.g. set gxout grfill
>            - set clevs
>            - check units
>
>            2. print to ascii
>            what you need is the fprintf.gs script
>            available
>       athttp://cookbooks.opengrads.org/index.php?title=Recipe-002:_Savi
>
>            ng_GrADS_variable_data_to_a_text_file
>
>            Copy the code and name the file as fprintf.gs.
>            Then do the following,
>
>            sdfopen input.nc
>            d var                   # where var = variable
>            fprintf var dat.txt      # prints data of var
>            fprintf lon lon.txt     # prints corresponding
>            lon
>            fprintf lat lat.txt     # prints corresponding
>            lat
>
>            3. statistical analysis
>            I know how to do so in R. If you're interested
>            in using R for this purpose, I can help.
>
>            HTH,
>
>            --
>            Muhammad Rahiz
>            Researcher & DPhil Candidate (Climate Systems
>            & Policy)
>            School of Geography & the Environment
>            University of Oxford
>
>            On Fri, 10 Dec 2010, James Ciarlo` wrote:
>
>            Dear all,
>            I am using GrADS 2.0.a9 on LINUX. I am
>            trying to compare measured and modelled
>            data by using GrADS. However since I am
>            new with the program, there are
>            some issues that I cannot seem to figure
>            out, or find a solution to understand.
>            I am experiencing:
>
>            1) Difficulty with displaying data from
>            NetCDF files
>            I have opened the file with the command
>            'sdfopen'
>            When I tried 'd prmsl' I got the
>            following error:
>            Operation error:  Invalid dimension
>            environment
>              Min longitude > max longitude: 320 144
>            So I used ' set lon -40 145' and tried
>            the display command again
>            I got the following line and the image
>            attached (ncimage.jpg):
>            Contouring: 99000 to 102300 interval 300
>            And I am getting similar results with
>            other nc files
>
>            2) Difficulty with finding a command to
>            output text files, or dat, or anything
>            similar.
>            I have used set lat, set lon, and set
>            lev to obtain one value, and set t to
>            get a long time interval, I have used d
>            psa to display the time series, I
>            can use printim file.jpg to get an
>            image, but I cannot find a command to
>            produce a file.txt so that I can make
>            use of the actual numbers for
>            statistical purposes.
>
>            3) Difficulty with finding a command to
>            work out comparative statistics,
>            particularly correlation and
>            significance.
>
>            Regards,
>
>            James Ciarlo`
> 
> 
> 
> 
> 
> 
> 
> 
>
-------------- next part --------------
'sdfopen sresb1.ncar_pcm1.tas.mon.nc'

** Set Lon
'set lon 77'

** Set Lat
'set lat 28'

** Set Level
*'set lev 500'

** Enter number of timesteps
ntot = 11 


n = 1
while(n < 11)

* Prints timestamp to file
'set t 'n''
'q dims'
r1 = sublin(result,5)
r2 = subwrd(r1,6)

of1 = 'x.time.txt'
res = write(of1,r2,append)
IO= sublin(res,1)
if(IO !=0)
 say 'write error #'IO
 'quit'
endif

* Prints data to file
'set gxout print'
'd tas'
r3 = sublin(result,2)
r4 = subwrd(r3,1)

of2 = 'x.data.txt'
res = write(of2,r4,append)
IO= sublin(res,1)
if(IO !=0)
 say 'write error #'IO
 'quit'
endif

n = n +1
endwhile

'quit'


More information about the gradsusr mailing list