Re: rcm2points

Hi Dennis,

Having looked at the Fortran source code for the rcm2points I've
learned two things...

1) It's the inverse distance weighted interpolation which is at the
source of the problem. The interpolation routine uses a parameter
(K=2) which determines the size of the spacing of the points to use
for interpolation ( [i,j] [i+K,j] [i,j+K] [i+K,j+K] ). This leads to
ugly jumps when switching from one box to the next when moving along
a line. This problem disappears by setting K=1 in the drcm2points.f
Fortran routine (see attached graphics invdist_k=1,2,3.png). Of
course, the inverse distance weighted interpolation still has a
tendency to be closest to the closest gridpoint and thus deviates
from the line when close to a source gridpoint (e.g.
invdist_k=1.png). I've replaced the inverse distance weighting by a
simple bi-linear interpolation (see code snippet from drcm2points.f
below) and the line is even better represented (e.g. bilin_k=1.png).

                        WX = (XO(NXY)-XI(IX,IY))/
      + (XI(IX+K,IY)-XI(IX,IY))
                        WY = (YO(NXY)-YI(IX,IY))/
      + (YI(IX,IY+K)-YI(IX,IY))
                        W(1,1) = (1.D0-WX)*(1.D0-WY)
                        W(2,1) = WX*(1.D0-WY)
                        W(1,2) = (1.D0-WX)*WY
                        W(2,2) = WX*WY

! W(1,1) = (1.D0/DGCDIST(YO(NXY),XO(NXY),
! + YI(IX,IY),XI(IX,IY),2))**2
! W(2,1) = (1.D0/DGCDIST(YO(NXY),XO(NXY),
! + YI(IX+K,IY),XI(IX+K,IY),2))**2
! W(1,2) = (1.D0/DGCDIST(YO(NXY),XO(NXY),
! + YI(IX,IY+K),XI(IX,IY+K),2))**2
! W(2,2) = (1.D0/DGCDIST(YO(NXY),XO(NXY),
! + YI(IX+K,IY+K),XI(IX+K,IY+K),2))**2

Conclusion: It would be nice to have a more versatile interpolation
function in NCL for regridding 2d fields in curvilinear coordinates.
As far as I understand linint2 does not do the trick for 2d
curvilinear grids. Maybe something else might work?

Note: The attached graphics have been generated without NCL by using
a small Fortran test program around the drcm2points.f routine in
order to emulate the small NCL test script I sent around in my
original post.

2) The routine should really only be applied for full latitudes/
longitudes as coordinates, as the the distance between two points is
computed as great circles. I've been using the routine with _rotated_
latitudes/longitudes (as our model uses a rotated lat/lon coordinate
system) and this will give erroneous (if only slightly) results.

Kind regards,
Oli

 
________________________________________

Oliver Fuhrer
Numerical Models

Federal Departement of Home Affairs FDHA
Federal Office of Meteorology and Climatology MeteoSwiss

Kraehbuehlstrasse 58, P.O. Box 514, CH-8044 Zurich, Switzerland

Tel. +41 44 256 93 59
Fax +41 44 256 92 78
oliver.fuhrer@meteoswiss.ch
www.meteoswiss.ch - First-hand information
  
 

> -----Original Message-----
> From: Dennis Shea [mailto:shea@ucar.edu]
> Sent: Dienstag, 4. Mai 2010 20:55
> To: Fuhrer Oliver
> Cc: ncl-talk@ucar.edu
> Subject: Re: rcm2points
>
> Hi Oli
>
> I'll look in the near future.
>
> I have added type of interpolation. Nothing fancy.
> It is a simple inverse distance squared algorithm.
> It was created for very dense grids like WRF.
>
> D
>
>
> On 5/4/10 8:12 AM, Oliver.Fuhrer@meteoswiss.ch wrote:
> > Hi NCL list,
> >
> > I am using the rcm2points interpolation routine to
> interpolate fields
> > from a curvilinear grid (lat/lon) to a list of specific
> locations given
> > in (lat/lon). Doing that, I've stumbled over a strange
> behaviour of the
> > rcm2points function. The NCL documentation page of
> rcm2points does not
> > specify what type of interpolation (bi-linear, nearest
> neighbour, ...)
> > rcm2points is supposed to do. I've attached an example
> script showing
> > the strange behaviour (see attached figure for output of
> script). In my
> > humble opinion, the output should be two almost identical
> lines (black
> > original, red interpolated). Even if rcm2points would do a nearest
> > neighbour interpolation, the resulting figure cannot really be
> > explained. The strange behaviour changes by increasing the number of
> > points in the mesh and disappears as a very large number of
> points is
> > used (e.g. s=(/100,100/) ).
> >
> > Any ideas?
> >
> > Cheers,
> > Oli
> >
> > ; *****************************************************************
> > ; test script showing off rcm2points problem
> >
> > begin
> >
> > ; define dimensions
> > s = (/4,4/) ; # points in mesh
> > np = 23 ; # points in line
> >
> > ; create regular mesh of unit square. the number of points
> > ; used in each direction is specified by s
> > x = new(s,float)
> > do i=0,s(1)-1
> > x(:,i) = fspan(0.0,1.0,s(0))
> > end do
> > y = new(s,float)
> > do i=0,s(1)-1
> > y(i,:) = fspan(0.0,1.0,s(1))
> > end do
> >
> > ; create straight line of np points from (x,y) =
> (0.1,0.1) to (0.8,0.9)
> > p = new((/np,2/),float)
> > tmp = fspan(0.,1.,np)
> > p(:,0) = 0.1 + tmp*0.7
> > p(:,1) = 0.1 + tmp*0.8
> > delete(tmp)
> >
> > ; interpolate mesh coordinates to line position
> > ; note: this should give back line
> > ; note: in normal usage, the third argument of rcm2points
> would be any
> > field to interpolate
> > pi = p
> > pi(:,0) = rcm2points(x,y,x,p(:,0),p(:,1),1)
> > pi(:,1) = rcm2points(x,y,y,p(:,0),p(:,1),1)
> >
> > ; open graphic port
> > wks = gsn_open_wks("ps","int_problem")
> >
> > ; plot original/interpolated line
> > r = True
> > r@gsnDraw = False
> > r@gsnFrame = False
> > pl1 = gsn_xy(wks,p(:,0),p(:,1),r)
> > r@xyLineColor = "red"
> > pl2 = gsn_xy(wks,pi(:,0),pi(:,1),r)
> >
> > ; overlay
> > overlay(pl1,pl2)
> > draw(pl1)
> > frame(wks)
> >
> > end
> >
> > ________________________________________
> >
> > Oliver Fuhrer
> > Numerical Models
> >
> > Federal Departement of Home Affairs FDHA
> > Federal Office of Meteorology and Climatology MeteoSwiss
> >
> > Kraehbuehlstrasse 58, P.O. Box 514, CH-8044 Zurich, Switzerland
> >
> > Tel. +41 44 256 93 59
> > Fax +41 44 256 92 78
> > oliver.fuhrer@meteoswiss.ch
> > www.meteoswiss.ch - First-hand information
> >
> >
> >
> >
> > _______________________________________________
> > ncl-talk mailing list
> > List instructions, subscriber options, unsubscribe:
> > http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk