[gradsusr] Help on calculating Area Maximum from NetCDF file using Grads

Wed Jul 16 11:45:55 EDT 2014

Dear Arlindo and All Grads Users,

As directed, I have developed  below Script, defined a global grid u4x4 and calculated Amax and used ‘set defval’ command. When I exported the defined grid to netcdf file. The output was vertical straight lines. Can anyone please point out what is wrong while defining ‘set defval’ in script? 

Any help please.

'reinit'

'open D:\WeatherGribFiles\2014\06\11\rtofs_glo_2ds_f048_daily_prog.ctl'

* original grid dimensions

xmin = 74.16

xmax = 434.0

ymin = -90

ymax = 90

* defining "u4x4" global grid for u

'define u4x4 = re(u_velocity,2)'

* defining var for further use

'define var = u_velocity'

xinc = 2

yinc = 1

xnow = xmin

ynow = ymin

******* jj and jj_max are the latitude of Global grid u4x4 ****************

jj = 1

jj_max = 90

count = 1

while (xnow < xmax)

xtempMax = xnow+xinc;ytempMax = ynow+yinc

'd amax(var,lon='xnow',lon='xtempMax',lat='ynow',lat='ytempMax')'

maxval=subwrd(result,4)

  jj = 1

   while (jj <= jj_max)

   'set defval u4x4 'count' 'jj' 'maxval

   jj = jj + 1

  endwhile

xnow = xnow+xinc;ynow = ynow+yinc

count=count+1

endwhile

*************** Write a netcdf file ****************

'set sdfwrite a_grid.nc'

'sdfwrite u4x4'

Thanks & Regards

Anil Kumar 

----------------------------------------

---------- Forwarded message ----------
From: Arlindo da Silva <dasilva at alum.mit.edu <mailto:dasilva at alum.mit.edu> >
Date: Thu, Jul 10, 2014 at 7:48 PM
Subject: Re: [gradsusr] Help on calculating Area Maximum from NetCDF file using Grads
To: Anil Kumar <anilku at oceanmanager.com <mailto:anilku at oceanmanager.com> >, GrADS Users Forum <gradsusr at gradsusr.org <mailto:gradsusr at gradsusr.org> >

On Thu, Jul 10, 2014 at 5:40 AM, Anil Kumar <anilku at oceanmanager.com <mailto:anilku at oceanmanager.com> > wrote:

Dear Arlindo,

Thanks for your detailed reply. I have a few queries regarding the task I need to accomplish.

1)      If Grads is not the best tool to accomplish my task, can you please suggest me the another tool, if you know, to accomplish the job ?

Probably C or Fortran. You could also accomplish some of that with python with pygrads. All depends on how much programming experience you have.

2)      Secondly, the approach which you have already shared with me as below

1) define a 2x2 global grid (say, use lterp with linear interpolation),

2) loop over every 2x2 region in the 0.25x0.5 deg grid, compute the 2 x2 area max with amax(), and use "set defval" to replace the corresponding value in the previous defined variable in 1). As I said, this is very inefficient.

Should I need to create 2x2 global grid first physically first and then need to replace its values. If so, I don’t find any hint on google,  <http://gradsuser.org> gradsuser.org to create a new grid. What is the procedure to create grid ?

What I said was to define a 2x2 grid, for example

define u2x2 = re(u,2)

   I really don't have the time to walk you thru the details as it would require me to actually write the code for you. If you can program in C or fortran I'd go that route. If you follow the UDF approach  that I outlined you can still do grads for your plotting. And please do not write directly to me, I cannot afford to provide individualized grads support.

   Good Luck,

        Arlindo

Thanks again for your help.

Thanks & Regards

Anil Kumar 

OceanManager Inc. |  <http://www.oceanmanager.com/> www.oceanmanager.com

OceanManager - Marine Software

From:  <mailto:gradsusr-bounces at gradsusr.org> gradsusr-bounces at gradsusr.org [mailto: <mailto:gradsusr-bounces at gradsusr.org> gradsusr-bounces at gradsusr.org] On Behalf Of Arlindo da Silva
Sent: Tuesday, July 01, 2014 7:12 AM
To: GrADS Users Forum

Subject: Re: [gradsusr] Help on calculating Area Maximum from NetCDF file using Grads

On Mon, Jun 30, 2014 at 9:37 AM, Anil Kumar <anilku at oceanmanager.com <mailto:anilku at oceanmanager.com> > wrote:

Hi Grads users,

I am awaiting an useful comments or direction on the below case. Please lead me, to get all the maximum values from the small – small areas of every 2 degree in Lat and lon say for example. Below is my complete query.

My Requirement :  I have a larger dataset (.NC file) 1440*720 points. I want a new destination grid with linear interpolation having all the all area maximum values using amax() of small-small regions. And ultimately, need to query the destination grid for output.

I am afraid the regridding method you are looking for is not implemented by lterp() or re().  This is a good example of when GrADS is not the best tool for the job at hand. But if you insist, a very inefficient approach would be

1) define a 2x2 global grid (say, use lterp with linear interpolation),

2) loop over every 2x2 region in the 0.25x0.5 deg grid, compute the 2 x2 area max with amax(), and use "set defval" to replace the corresponding value in the previous defined variable in 1). As I said, this is very inefficient.

If you are up for some C/Fortran programming, you can use the opengrads IPC extension (which is available with your Windows build):

    http://opengrads.org/doc/udxt/libipc/

to write your variables to disk,

ga-> ipc_save u_velocity u.bin

ga-> ipc_save v_velocity v.bin

Then write a C/Fortran program that read these files (u.bin,v.bin) and compute the amax at the 2x2 grid you require, writing output files (u2x2.bin.v2x2.bin). 

ga-> ! my_amaxregrid 

Consult the IPC documentation for the description of the transfer file format (very similar to the old GrADS 1.x UDF.) You can finally read these back with ipc_load() extension:

ga-> d ipc_load("u2x2.bin");ipc_load("v2x2.bin")

Notice that this is a poor-man version of the file-based UDF that existed in GrADS 1.x. Of course, you can accomplish the same more efficiently by writing a genuine OpenGrADS extension with no transfer files involved.

  Good Luck,

      Arlindo

What Am doing so far : 

I have a NetCDF file having  below CTL Info.

*****************************************************************

dset D:\WeatherGribFiles\2014\06\22\ <http://rtofs_glo_2ds_f024_daily_prog.nc> rtofs_glo_2ds_f024_daily_prog.nc

title 2D Sfc Daily Prognostic 00Z22jun2014: Forecast, downloaded Jun 22 15:20 UTC

undef 1.26765e+30

dtype netcdf

xdef 1440 linear 74.16 0.24999

ydef 720 linear -90 0.24999

zdef 1 linear 1 1

tdef 1 linear 00Z23JUN2014 1mn

vars 2

u_velocity=>u_velocity  1  t,z,y,x  eastward_sea_water_velocity (m/s)

v_velocity=>v_velocity  1  t,z,y,x  northward_sea_water_velocity (m/s)

endvars

*****************************************************************

I want a new destination grid ; which would contain the area maximum values as points. Destination grid should contain maximum value from every 8 points in x-axis and y axis. For this, I have created below, CTL file named as  rtofs_glo_2ds_f024_daily_prog.ctl:

*****************************************************************

dset D:\WeatherGribFiles\2014\06\22\ <http://rtofs_glo_2ds_f024_daily_prog.nc> rtofs_glo_2ds_f024_daily_prog.nc

index  D:\WeatherGribFiles\2014\06\22\rtofs_glo_2ds_f024_daily_prog_DUMMY.nc 

title 2D Sfc Daily Prognostic 00Z22jun2014: Forecast, downloaded Jun 22 15:20 UTC 

undef 1.26765e+30 

dtype netcdf 

xdef 180 linear 74.16 2.0  

ydef 90 linear -90 2.0  

zdef 1 linear 1 1 

tdef 1 linear 00Z23JUN2014 1mn

options template 

vars 2 

u_velocity=>u_velocity  1  t,z,y,x  eastward_sea_water_velocity (m/s)  

v_velocity=>v_velocity  1  t,z,y,x  northward_sea_water_velocity (m/s)  

endvars

*****************************************************************

Below is the command I use to display the data from destination grid rtofs_glo_2ds_f024_daily_prog_DUMMY.nc defined in above CTL file. Below is the command set using.

*****************************************************************

reinit

sdfopen D:/WeatherGribFiles/2014/06/22/ <http://rtofs_glo_2ds_f024_daily_prog.nc> rtofs_glo_2ds_f024_daily_prog.nc

open D:/WeatherGribFiles/2014/06/22/rtofs_glo_2ds_f024_daily_prog2.ctl

set parea 0 11 0 8.5

set vpage 0 11 0 8.5

set lon -74.16 434

set lat -90 90

set cthick 6

set arrscl 0.10 

define val = mag(lterp(u_velocity.1*1.95,lat.2,amax),lterp(v_velocity.1*1.95,lat.2,amax))

d lterp(u_velocity.1*1.95,lat.2,amax);lterp(v_velocity.1*1.95,lat.2,amax);val

*****************************************************************

When I change the destination grid dimensions slightly  in rtofs_glo_2ds_f024_daily_prog.ctl file

xdef 135 linear 74.16 2.5  

ydef 68 linear -90 2.5  

and checks for the output, It varies at the same point i.e. if I look at certain point before changing and after changing grid dimensions, the result comes of different magnitude.

What I want is, When I increase or decrease the destination grid points in “rtofs_glo_2ds_f024_daily_prog.ctl”, and the maximum value should be at its place. 

How can I achieve this ????? 

Thanks & Regards

Anil Kumar – Sr. Software Developer