force grads to interpolate

[Bernd Becker]

Patrick,

you may like to try this:


function interpolate(argument_list)

locale=subwrd(argument_list,1)                     * name of the output
 part1=subwrd(argument_list,2)                     * variable
  llat=subwrd(argument_list,3)                     * latitude
  llon=subwrd(argument_list,4)                     * longitude

*say locale' 'part1' 'llat' 'llon

* learn current dimension setting

'q dims'
line2=sublin(result,2)
line3=sublin(result,3)
*say line2
origlon1=subwrd(line2,6)
origlon2=subwrd(line2,8)

*say line3
origlat1=subwrd(line3,6)
origlat2=subwrd(line3,8)

'q file'
line5=sublin(result,5)
xmax=subwrd(line5,3)
xmax=0.+xmax

* convert lat  lon into x y space

*say 'q w2gr 'llon' 'llat
    'q w2gr 'llon' 'llat
say result
xdim=subwrd(result,3)
xdim=0.+xdim
if ( xdim <  1.) ; xdim=xmax+xdim ; endif
ydim=subwrd(result,6)

* nearest x dimensions :

x1= math_int(xdim)
x2 = x1 + 1
if ( x2 > xmax ) ; x2  = x2 - xmax ; endif
*say "xdim = " xdim' 'x1' 'x2

* the weights are

xw1= xdim-x1
xw2= 1.-xw1

* nearest y dimensions :

y1= math_int(ydim)
y2 = y1 + 1

* the weights are

yw1= ydim-y1
yw2= 1.-yw1

*part4=       part1%'(x='x1',y='y1')*'xw2'*'yw2'+'
*part4= part4%part1%'(x='x2',y='y1')*'xw1'*'yw2'+'
*part4= part4%part1%'(x='x1',y='y2')*'xw2'*'yw1'+'
*part4= part4%part1%'(x='x2',y='y2')*'xw1'*'yw1

_part2= '(x='x1',y='y1')*'xw2'*'yw2'+'
_part3= '(x='x2',y='y1')*'xw1'*'yw2'+'
_part5= '(x='x1',y='y2')*'xw2'*'yw1'+'
_part6= '(x='x2',y='y2')*'xw1'*'yw1

*  station =           f(x1,y1)*xw2*yw2+
*                      f(x2,y1)*xw1*yw2+
*                      f(x1,y2)*xw2*yw1+
*                      f(x2,y2)*xw1*yw1

      part4=       part1%_part2
      part4= part4%part1%_part3
      part4= part4%part1%_part5
      part4= part4%part1%_part6

'set x 1'
'set y 1'
*say 'define 'locale' = 'part4
'define 'locale' = 'part4
*'d 'locale
*say subwrd(result,4)

* return to dimension setting as  before interpolation

'set lat 'origlat1' 'origlat2
'set lon 'origlon1' 'origlon2

return




call with:

rc = interpolate( local, par ,lati, loni)


On Sun, 2007-12-02 at 19:23 +0100, Patrick Reuter wrote:
> Dear Brian,
>
> thanks for your answer. Actually, I was just making this simple
> example to get easily understood.
>
> My actual problem is, that I have floating point latitude/logitude
> data for a city, and a WRF model step of 0.0791 latitudes/longitudes.
>
> That is why is was talking about linear interpolation : Taking the 4
> gridpoints that surround the city, and weight the influence of the
> gridpoints wrt the distance.
>
> Does this linear interpolation exist ? Or do I have to find the 4
> gridpoints on my own and do the interpolation ?
>
> I hope this was clear ..
>
> Thanks in advance
>
>     Patrick
>
>
>
> Brian Doty <doty at COLA.IGES.ORG> a écrit :
>
> > If you have one degree data, and want to "interpolate" to lon 90.5,
> > you can use the expression:
> >
> >     display (var(lon=90)+var(lon=91))/2
> >
> > ...Brian
> >
> > On Dec 1, 2007, at 12:39 PM, Patrick Reuter wrote:
> >
> >> Hi everybody,
> >>
> >> maybe this question is pretty basic - but I just can't find it.
> >>
> >> ga-> set lat 90.5
> >> LAT set to 91 91
> >>
> >> How could I force grads to use linear interpolation between the values
> >> ? I mean, that it takes the values of grid point 90 and 91, and
> >> linearly interpolates, so that I can display the (one dimensional)
> >> data ?
> >>
> >> Thanks in advance !
> >>
> >>  Patrick
--
Bernd Becker   The Monthly Outlook
Met Office  FitzRoy Road  Exeter   Devon EX1 3PB  United Kingdom
Tel.: +44 (0) 1392 884511 Fax: +44 (0)870 900 5050
E-mail:bernd.becker at metoffice.com  -  http://www.metoffice.com