[ncl-talk] NCL spline question

Wed Dec 23 16:01:58 MST 2015

Greetings NCL-Talk,

I'm trying to generate an interpolatory spline* for a lat/lon tropical 
cyclone trajectory from a set of lat/lon points at semi-regular time 
intervals. The input points are regularly spaced (for the most part) 
6-hourly Best Track points, however there are occasional irregular time 
points in between, such as the exact time of landfall. The desired 
output is a semi-regular grid of hourly points, but with the same few 
irregular time points as the for the input times.

(*An interpolatory spline is one in which the interpolated spline passes 
through the points that are being interpolated from, as opposed to an 
approximating spline which does not have to pass through the points.)

It was suggested that NCL's ftkurv function (from the FitGrid library) 
might do what I'd like, since this does an interpolation for parametric 
curves. If I understand my problem correctly, lat and lon can be thought 
of as parameters of time.

The example plot for ftkurv looks promising:
http://www.ncarg.ucar.edu/ngmath/fitgrid/plot4.html

ftkurv (
                 xi [*] : numeric,
                 yi [*] : numeric,
                 t  [*] : numeric,
                 xo [*] : float,    ; or double
                 yo [*] : float     ; or double
         )

In the example given for this function, the parameter array passed in 
for 't' is a regularly-spaced normalized parameter array that ranges 
from 0 to 1. I don't understand where this function gets any information 
about the input time values however. I built a test case for a real 
hurricane track and validated the interpolated spline points at the 
input points. I get zero "error" at the end points of the track, with 
maximum error in the middle. The error increases and decreases smoothly, 
suggesting an offset issue (probably related to time drift).

- Does anybody know if this function allow for irregular parameter 
values for 't'?
- Is there a way to use ftkurv in a way that makes it aware of the times 
of the input lat/lon?

I've created a fairly detailed stand-alone test code, which is attached. 
Feel free to run to see the error characteristics of the ftkurv spline.

Ultimately, I'd like to be able to do piecewise cubic spline 
interpolation along the lines of what is done for gps tracks: 
http://topofusion.com/spline.php (piecewise interpolation, matching 
derivatives to ensure continuity of slope, Bessel interpolation used to 
compute the first derivative by fitting a 3-point parabola for each 
point). If someone has already implemented something like this in NCL, 
that would of course be awesome. But I can make do with something more 
basic provided it can handle the irregular input/output times.

If anyone has experience with this type of spline problem, or has 
suggestions on other NCL routines to try, I'd appreciate any input you 
can provide.

Jonathan
-------------- next part --------------
;********************************************************
; module_spline_helper_routines.ncl
;********************************************************
; This module contains procedures to add polymarkers, 
; polylines, and polygons to a plot, one at a time. This 
; module is quite helpful in that it solves the issue of
; ensuring that each poly-item have a unique identifier,
; saving the user from having to worry about the hassle of
; creating arrays, etc. There is also a time conversion
; routine and a function which creates a variety of
; new polymarkers.
;
;
; The first group are essentially wrappers to the gsn_add_* functions 
; and accept limited input as arguments (instead of resources; 
; add_text is an exception to this).
;
;   procedure  add_text
;   procedure  add_polyline
;   procedure  add_polymarker
;   procedure  add_polygon
;
;
; The next functions are used to create new polymarker styles and 
; set some things specific to certain plots.
;
;   function   create_new_polymarkers
;
;
; The next functions do time conversion from timeoffset since
; base time to yyyymmddhh format:
;
;   function   convert_timeoffset_to_yyyymmddhh
;   function   convert_timeoffset_to_yyyymmddhhmm
;
; I believe Mary Haley provided the essential form of add_*
; procedures. I packaged them into this module and added 
; a few bells and whistles so I can use them with greater 
; ease in various codes.
;
; NOTICE of COPYRIGHT: This script was developed by 
; Jonathan Vigh while he was at Colorado State University 
; (with subsequent development at the National Center for 
; Atmospheric Research). Mary Haley (University Corporation 
; for Atmospheric Research provided some assistance). 
;
; All Rights Reserved by author. Please do not redistribute this 
; script without the author's permission. Use of this script
; for commercial expressly prohibited.
;
; ~ Jonathan Vigh, National Center for Atmospheric Research
;********************************************************

;********************************************************
; procedure add_text
;
; This procedure adds text to a plot (in plot coordinates), making sure that each
; set is returned to a unique variable name, and that this
; variable is retained even outside this procedure call. 
;********************************************************
undef("add_text")
procedure add_text(wks[1]:graphic,plot[1]:graphic,text[1]:string,x[1]:numeric,y[1]:numeric,txres[1]:logical)

local txres

begin

  textid = unique_string("text")  ; "unique_string" will return a unique
                                  ; string every time it is called from
                                  ; within a single NCL session.