fwrite wrap data

[Eric Altshuler]

Javier,

There is another situation where unwanted behavior can occur with fwrite. If you fwrite a global grid, be sure to set the dimension environment using 'set x' and 'set y' before doing the fwrite. For example, with a 0.5 degree (720x361) global grid, when you open the file in grads (and if it's the first file you open in a grads session), the initial dimension environment will be set to the following by default:

lon: 0 to 360 (X from 1 to 721)
lat: -90 to 90 (Y from 1 to 361)

If you fwrite without changing the dimension environment, the resulting file will have 721x361 elements because the data for lon=0 (X=1) will be duplicated as an extra column (lon=360, X=721). To prevent this behavior and obtain a 720x361 array in the fwrite file, do this before the fwrite:

set x 1 720

This will (in this case) set lon from 0 to 359.5 and prevent the wraparound behavior.

Best regards,

Eric L. Altshuler
Assistant Research Scientist
Center for Ocean-Land-Atmosphere Studies
4041 Powder Mill Road, Suite 302
Calverton, MD 20705-3106
USA

E-mail: ela at cola.iges.org
Phone: (301) 902-1257
Fax: (301) 595-9793

----- Original Message -----
From: "Javier G. Corripio" <Javier.Corripio at UIBK.AC.AT>
To: GRADSUSR at LIST.CINECA.IT
Sent: Monday, September 28, 2009 5:05:31 PM GMT -05:00 US/Canada Eastern
Subject: Re: fwrite wrap data

Thanks a lot!

Spot on, Erik, it was that extra column multiplied by every row that made all the trouble.

And thanks Michael, when I checked by q dim on a higher resolution set, set lat and set lon gives
fractional values for X and Y which may cause trouble later.

I owe you two a pint of beer (or many more) if we ever met.

Cheers.

Javier


On Mon, 28 Sep 2009 16:37:36 -0400, Eric Altshuler <ela at COLA.IGES.ORG> wrote:

>From looking at the plots, it seems that this is more than just a wraparound problem. The fields
show significant differences aside from the wraparound. One possible problem might be the
number of grid points in the xdef and ydef lines in the ctl file for the subset. The dimensions
should be 41 and 21 rather than 40 and 20.
>
>Best regards,
>
>Eric L. Altshuler
>Assistant Research Scientist
>Center for Ocean-Land-Atmosphere Studies
>4041 Powder Mill Road, Suite 302
>Calverton, MD 20705-3106
>USA
>
>E-mail: ela at cola.iges.org
>Phone: (301) 902-1257
>Fax: (301) 595-9793
>
>----- Original Message -----
>From: "Michael L CIV 63134 Sestak" <michael.l.sestak at NAVY.MIL>
>To: GRADSUSR at LIST.CINECA.IT
>Sent: Monday, September 28, 2009 4:17:31 PM GMT -05:00 US/Canada Eastern
>Subject: Re: fwrite wrap data
>
>
>
>When using fwrite you should use 'set x' and 'set y' instead of 'set lat' and 'set lon' so the number
of grid points in each dimension is certain. With set lat and set lon, what a person thinks is the last
x grid position may not be written due to roundoff error, especially for fine resolution grids.
>
>Michael Sestak
>Fleet Numerical Meteorology and Oceanography Center
>Monterey, CA
>
>
>
>From: Javier G. Corripio
>Sent: Mon 9/28/2009 12:03
>To: GRADSUSR at LIST.CINECA.IT
>Subject: fwrite wrap data
>
>
>Dear all,
>
>I would appreciate some advice with fwrite.
>
>I am trying to write some data which is a subset of a larger grib.  It seems that the data is
wrapped
>around in a peculiar way when using fwrite.
>
>Here I show and example for clarity, the topography of the Alps from the GFS (hgtsfc). When
using
>the data written by fwrite it will display with the eastern end on the northwestern corner
>
>The ctl files seems to be OK, so I wonder if I failed to specify any settings for fwrite.
>
>These are the commands:
>
>open gfs.ctl
>set lon 0 20
>set lat 40 50
>d hgtsfc;************** first map correct
>set gxout fwrite
>set fwrite alpsset.dat
>d hgtsfc
>disaple fwrite
>reinit
>open alpsset.ctl
>d topo;**************** second map wrapped around
>
>
>ctl file for subset:
>
>ga-> q ctlinfo
>dset alpsset.dat
>title
>undef 9.999e+20
>xdef 40 linear 0 0.5
>ydef 20 linear 40 0.5
>zdef 1 linear 1013 1
>tdef 1 linear 00Z04SEP2009 60mn
>vars 1
>topo  0  1,1,0,0  surface height (m) []
>endvars
>
>
>
>original ctl file:
>
>ga-> q ctlinfo 2
>dset gfs
>title gfs
>undef 9.999e+20
>xdef 720 linear 0 0.5
>ydef 361 linear -90 0.5
>zdef 26 levels 1000 975 950 925 900 850 800 750
> 700 650 600 550 500 450 400 350 300 250
> 200 150 100 70 50 30 20 10
>tdef 51 linear 00Z28SEP2009 180mn
>vars 52
>.....
>endvars
>
>Image of what it should be and what I obtained:
>
>http://tinyurl.com/y8f2xcq
>
>
>or in a schematic way,
>
>THIS:
>
>*
>****
>********
>************
>*****************
>****************************
>*****************
>************
>********
>****
>*
>
>
>BECOMES THIS:
>
>***
>**************
>***
>*
>****
>********
>************
>**************
>**************
>**************
>************
>********
>****
>*
>
>
>Thanks a lot for any help.
>
>Javier