GrADS to Ascii with header

[Stefan Fronzek]

Hi,

an alternative method to python is implemented in the grads script below
that writes out grads expressions to ESRI's ascii grid format. The script is
modified from an earlier post on this list.

Best regards,
Stefan



***********
* Name:    print.asciigrid
* Args:    <expression> <preout>
* Purpose: prints out the given expression given by the argument to an
*          ASCII grid in ESRI format that can be imported by ArcView/ArcInfo.
*          It will write a file called <preout>ddMMMyyyy.asc for each time step.
*          ESRI's ASCII grids are defined only on quadratic grid cells; the
*          result will be distorted if the lat- and lon-resolution of the GrADS
*          expression are not equal and if they are not linear.
*
*          This script is modified from Daniel de Castro Victoria's posting
*          on the gradsusr list.
* Author:  Stefan Fronzek, e-mail: Stefan(dot)Fronzek(at)ymparisto(dot)fi
* Version: 1.2 (Tue Jan 3 16:34:16 FLEDT 2006)
***********

function main(args)

var=subwrd(args,1)
pre=subwrd(args,2)
timefmt=subwrd(args,3)

if(var = '')
say "print.asciigrid - exporting grads grids to  ASCII grid (ESRI)"
say "                  output file names: <preout>ddMMMyyyy.asc"
say "Usage: run print.asciigrid <expr> <preout> [mm]"
say "       where the optional mm-flag changes the default file names"
say "       to <preout>01.asc for January, <preout>02.asc for Feb ..."
return
endif

say "print.asciigrid - exporting ("var") with the following dimensions:"
say
'q dims'
say result

* get the X, Y, and T dimensions
'q dims'
lx=sublin(result,2)
ly=sublin(result,3)
lt=sublin(result,5)

minlon=subwrd(lx,6)
maxlon=subwrd(lx,8)
minx=subwrd(lx,11)
maxx=subwrd(lx,13)
* don't allow fixed X or Y dimensions as the resolution would be undefined
if(subwrd(lx,3) = 'fixed')
say 'X dimension may not be fixed'
return
endif
minlat=subwrd(ly,6)
maxlat=subwrd(ly,8)
miny=subwrd(ly,11)
maxy=subwrd(ly,13)
if(subwrd(ly,3) = 'fixed')
say 'Y dimension may not be fixed'
return
endif
initime=subwrd(lt,11)
fimtime=subwrd(lt,13)
if(subwrd(lt,3) = 'fixed')
initime=subwrd(lt,9)
fimtime=initime
endif

* calculate the resolution
res=(maxlat-minlat)/(maxy-miny)

** time-step loop
ti=initime
while (ti <= fimtime)
   'set t 'ti
*  set the filename to the time: <ddMMMyyyy>.asc
   'q time'
   nomel=subwrd(result,3)
   nome=substr(nomel,4,9)
   if(timefmt='mm')
     month=substr(nomel,6,3)
     mm='01'
     if(month='FEB') ; mm='02' ; endif
     if(month='MAR') ; mm='03' ; endif
     if(month='APR') ; mm='04' ; endif
     if(month='MAY') ; mm='05' ; endif
     if(month='JUN') ; mm='06' ; endif
     if(month='JUL') ; mm='07' ; endif
     if(month='AUG') ; mm='08' ; endif
     if(month='SEP') ; mm='09' ; endif
     if(month='OCT') ; mm='10' ; endif
     if(month='NOV') ; mm='11' ; endif
     if(month='DEC') ; mm='12' ; endif
     nome=pre''mm'.asc'
  else
     nome=pre''nome'.asc'
   endif

   say "Exporting "nome" ..."

*  write the header information
   lixo=write(nome,"ncols       "maxx-minx+1)
   lixo=write(nome,"nrows       "maxy-miny+1)
   lixo=write(nome,"xllcenter   "minlon)
   lixo=write(nome,"yllcenter   "minlat)
   lixo=write(nome,"cellsize    "res)
   lixo=write(nome,"NODATA_value        -9999")

*  loop Y - line or lat (variable j)
   j=maxy
   while (j>=miny)
*      *isolate a line
      'set y 'j
*      *loop X - column or lon (variable i)
      i=minx
      linha=""
      while (i<=maxx)
*         *isolate a columns
         'set x 'i
*         *get the date and write to the end of a string
         'd 'var
         temp=subwrd(result,4)
         linha=linha''temp' '
         i=i+1
      endwhile
*      *write one line of data
      lixo=write(nome,linha)
      j=j-1
   endwhile

   lixo=close(nome)
   ti=ti+1
endwhile
* reset the dimensions
'set x 'minx' 'maxx
'set y 'miny' 'maxy
'set t 'initime' 'fimtime





On Fri, 19 Dec 2008 22:40:18 -0200, Arlindo da Silva
<arlindo.dasilva at GMAIL.COM> wrote:

>On Fri, Dec 19, 2008 at 3:17 AM, Ilukpitiya Mudi S P Jayawardena <
>imspj at hawaii.edu> wrote:
>
>> Dear,
>>
>>
>>
>> Last time, I asked how to convert GrADS to Ascii. So you have replied me.
>> Thank you very much for replying.
>>
>> Your method is successfully working. Here ascii file is generating without
>> a header.
>>
>> My ultimate purpose is utilizing this ascii file in ArcGIS. Therefore, I
>> want the ascii file with a header.
>>
>>
>>
>> Pls tell me how to add a header to the ascii file in the process of
>> converting GrADS to Ascii.
>>
>> Pls reply me.
>>
>
>If you are willing to use python, the method described in Recipe-003 is
>quite flexible:
>
>http://cookbooks.opengrads.org/index.php?title=Table_of_Contents#SEQUENTIAL_INDEX_OF_RECIPES
>
>As illustrated in this recipe, you can even do this using data on an OPeNDAP
>server without the need for downloading the data first.
>
>    Arlindo
>
>
>
>>
>>
>> With best rgds,
>>
>> Shiromani
>>
>
>
>
>--
>Arlindo da Silva
>dasilva at alum.mit.edu
>