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 _rotate_
latitudes/longitudes (as our model uses a rotated lat/lon coordinate
system) and this will give erroneous (if only slightly) results.

Kind regards,
Oli

-----Original Message-----
From: Dennis Shea [mailto:shea@ucar.edu]
Sent: Tue 5/4/2010 8:55 PM
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