# [ncl-talk] NCL spline question

Jonathan Vigh jvigh at ucar.edu
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).
;
;
;
; 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).
;
; script without the author's permission. Use of this script
; for commercial expressly prohibited.
;
; ~ Jonathan Vigh, National Center for Atmospheric Research
;********************************************************

;********************************************************
;
; 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.
;********************************************************

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.

; You can then use this unique string as an attribute variable name
; that gets attached to the plot variable. This ensures that this
; value will live for the duration of the script.

end

;********************************************************
;
; This procedure adds markers to a plot, making sure that each
; set is returned to a unique variable name, and that this
; variable is retained even outside this procedure call.
;********************************************************

local pmres, str

begin

pmres = True

pmres at gsMarkerColor = color
pmres at gsMarkerIndex = type
pmres at gsMarkerSizeF = size

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

; You can then use this unique string as an attribute variable name
; that gets attached to the plot variable. This ensures that this
; value will live for the duration of the script.

end

;********************************************************
;
; This procedure adds a line to a plot, making sure that each
; set is returned to a unique variable name, and that this
; variable is retained even outside this procedure call.
;********************************************************

local plres,polyline

begin

plres = True

plres at gsLineColor 	  = color
plres at gsLineThicknessF  = thickness
plres at gsLineDashPattern = dash_pattern
plres at gsLineDashSegLenF = dash_segment_length

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

; You can then use this unique string as an attribute variable name
; that gets attached to the plot variable. This ensures that this
; value will live for the duration of the script.

end

;********************************************************
;
; This procedure adds a polygon to a plot, making sure that each
; set is returned to a unique variable name, and that this
; variable is retained even outside this procedure call.
;********************************************************

local pyres,polygon

begin

pyres = True

pyres at gsFillColor = color

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

; You can then use this unique string as an attribute variable name
; that gets attached to the plot variable. This ensures that this
; value will live for the duration of the script.

end

;********************************************************

; This procedure adds text labels to a plot (in data coordinates),
; making sure that each set is returned to a unique variable name, and that this
; variable is retained even outside this procedure call.
;
; Unlike the add_text procedure above, this routine can add multiple
; text items.
;********************************************************

local txres

begin

if ( (dimsizes(x) .eq. dimsizes(y)) .and. (dimsizes(x) .eq. dimsizes(text)) .and. (dimsizes(x) .ge. 1) ) then
ntext = dimsizes(x)
else
print("ERROR: Jgsn_add_text needs at least one x and y points to draw the text. The number of x and y points, and text strings must all be the same size.")
exit
end if

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

; You can then use this unique string as an attribute variable name
; that gets attached to the plot variable. This ensures that this
; value will live for the duration of the script.

end do

end

;********************************************************
;
; This procedure adds markers to a plot, making sure that each
; set is returned to a unique variable name, and that this
; variable is retained even outside this procedure call.
;
; Unlike the add_marker procedure above, this routine can add multiple
; polymarkers.
;********************************************************

local pmres, str

begin

if ( (dimsizes(x) .eq. dimsizes(y)) .and. (dimsizes(x) .ge. 1) ) then
nmarkers = dimsizes(x)
else
print("ERROR: Jgsn_add_polymarker needs at least one x and y points to draw a line. The x and y inputs must be the same size.")
exit
end if

do im = 0,nmarkers-1
polymarker = unique_string("polymarker")  ; "unique_string" will return a unique
; string every time it is called from
; within a single NCL session.

; You can then use this unique string as an attribute variable name
; that gets attached to the plot variable. This ensures that this
; value will live for the duration of the script.

end do

end

;********************************************************
;
; This procedure adds a line to a plot, making sure that each
; set is returned to a unique variable name, and that this
; variable is retained even outside this procedure call.
;
; Unlike the add_polyline procedure above, this routine can add multiple
; polylines.
;********************************************************

local plres,polyline

begin

plres = True

; You can then use this unique string as an attribute variable name
; that gets attached to the plot variable. This ensures that this
; value will live for the duration of the script.

if ( (dimsizes(x) .eq. dimsizes(y)) .and. (dimsizes(x) .ge. 2) ) then
npts = dimsizes(x)
else
print("ERROR: Jgsn_add_polyline needs at least two x and y points to draw a line. The x and y inputs must be the same size.")
exit
end if

do ip = 1, npts - 1

if (.not.ismissing(x(ip - 1)) .and.     \
.not.ismissing(x(ip)) .and.         \
.not.ismissing(y(ip - 1)) .and.     \
.not.ismissing(y(ip))) then

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

plot@\$polyline\$ = gsn_add_polyline(wks,plot,x(ip - 1:ip),y(ip - 1:ip),plres)

end if

end do

end

;********************************************************
;
; This procedure adds a polygon to a plot, making sure that each
; set is returned to a unique variable name, and that this
; variable is retained even outside this procedure call.
;********************************************************

local pyres,polygon

begin

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

; You can then use this unique string as an attribute variable name
; that gets attached to the plot variable. This ensures that this
; value will live for the duration of the script.

end

;********************************************************
; Function: create_new_polymarkers
;
; This function creates some new marker styles and returns
; a logical variable with the marker identifiers attached.
;
; It has been adapted from an example from the NCL web page.
;********************************************************
undef("create_new_polymarkers")
function create_new_polymarkers(wks)

begin

MARKERS = True

;************************************************
; Create our own markers using NhlNewMarker
;************************************************
; the arguments for this function are:
; wks
; marker_string[*]
; font_table_number
; x-offset
; y-offset
; aspect_ratio
; size
; angle

; this example will create filled squares. You will have to play with
; the numbers a but to get the size and shape you desire. On the
; documentation page for NhlNewMarker, there is a table of values for
; the current marker set, to give you an idea of where to start.

;index 	description 			character 	font table 	X offset 	Y offset 	aspect ratio 	size 	angle
;1 	dot (small filled circle) 	Z 		37 		0.0 		0.0 		1.3125 		0.175 	0.0
;2 	plus sign 			+ 		18 		0.0 		0.075 		1.3125 		0.95 	0.0
;3 	asterisk 			* 		1 		0.0 		0.0 		1.3125 		1.0 	0.0
;4 	hollow circle 			x 		19 		0.0 		0.075 		1.3125 		1.2 	0.0
;5 	cross (x) 			U 		18 		0.0 		0.075 		1.3125 		1.1 	0.0
;6 	hollow square 			Z 		19 		0.0 		0.083 		1.3125 		1.45 	0.0
;7 	up pointing triangle 		[ 		19 		0.0 		0.03 		1.5 		1.25 	0.0
;8 	down pointing triangle 		X 		19 		0.0 		0.87 		2.15 		0.67 	0.0
;9 	diamond 			\\ 		19 		0.0 		0.075 		1.0 		1.15 	0.0
;10 	left pointing filled triangle 	` 		19 		0.0 		0.08 		1.5 		1.55 	0.0
;11 	right pointing filled triangle 	b 		19 		0.0 		0.08 		1.5 		1.55 	0.0
;12 	five-pointed star 		] 		19 		0.0 		0.0625 		1.3125 		1.1 	0.0
;13 	six-pointed star 		m 		19 		0.0 		0.0725 		1.3125 		1.1 	0.0
;14 	circle with center dot 		Z 		18 		0.0 		0.0 		1.3125 		0.8 	0.0
;15 	circle with cross 		[ 		37 		0.0 		0.0 		1.3125 		0.8 	0.0
;16 	filled circle 			Z 		37 		0.0 		0.0 		1.3125 		0.8 	0.0

; also set some recognizable names for predefined markers or define new ones
MARKERS at hollow_circle 		= 4
MARKERS at filled_circle 		= 16

MARKERS at hollow_square 		= 6

MARKERS at filled_square                 = NhlNewMarker(wks, "y", 35, 0.0, 0.0, 1.0, 0.9, 0.0)

MARKERS at hollow_up_pointing_triangle   = 7
MARKERS at hail_symbol 			= 7	; hollow upward-pointing triangle
MARKERS at filled_up_pointing_triangle   = NhlNewMarker(wks, "u", 34, 0.0, 0.07, 1.3125, 1.3, 0.0)

MARKERS at hollow_down_pointing_triangle = 8
MARKERS at filled_down_pointing_triangle = NhlNewMarker(wks, "u", 34, 0.0, -0.07, 1.3125, 1.3, 180.0)

MARKERS at hollow_diamond		= 9
MARKERS at filled_diamond		= NhlNewMarker(wks, "q", 35, 0.0, 0.0, 1.0, 1.0, 0.0)

MARKERS at hollow_star			= NhlNewMarker(wks, "]", 19, 0.0, 0.0, 1.0, 1.0, 0.0)
MARKERS at filled_star			= NhlNewMarker(wks, "z", 35, 0.0, 0.0, 1.0, 1.0, 0.0)

MARKERS at hollow_cloverleaf		= NhlNewMarker(wks, "{", 19, 0.0, 0.0, 1.0, 1.0, 0.0)
MARKERS at filled_cloverleaf		= NhlNewMarker(wks, "p", 35, 0.0, 0.0, 1.0, 1.0, 0.0)

MARKERS at hollow_heart			= NhlNewMarker(wks, "%", 19, 0.0, 0.0, 1.0, 1.0, 0.0)
MARKERS at filled_heart			= NhlNewMarker(wks, "r", 35, 0.0, 0.0, 1.0, 1.0, 0.0)

MARKERS at hollow_spade			= NhlNewMarker(wks, ".", 19, 0.0, 0.0, 1.0, 1.0, 0.0)
MARKERS at filled_spade			= NhlNewMarker(wks, "s", 35, 0.0, 0.0, 1.0, 1.0, 0.0)

MARKERS at tropical_storm_symbol         = NhlNewMarker(wks, "m", 35, 0.0, 0.0, 1.0, 0.7, 0.0)
MARKERS at hurricane_symbol		= NhlNewMarker(wks, "p", 37, 0.0, 0.0, 1.0, 0.8, 0.0)

MARKERS at lightning_symbol		= NhlNewMarker(wks, "l", 35, 0.0, 0.0, 1.0, 0.8, 0.0)
MARKERS at warm_ring_symbol		= NhlNewMarker(wks, "F", 34, 0.0, 0.0, 1.0, 0.7, 0.0)
MARKERS at elliptic_symbol		= NhlNewMarker(wks, "y", 19, 0.0, 0.0, 1.5, 0.6, -45.0)
MARKERS at open_eye_symbol		= NhlNewMarker(wks, "l", 19, 0.0, 0.0, 1.0, 1.3, 0.0)
MARKERS at closed_eye_symbol		= NhlNewMarker(wks, "0", 34, 0.0, 0.0, 1.0, 1.7, 0.0)
MARKERS at concentric_symbol		= NhlNewMarker(wks, "c", 35, 0.0, 0.0, 1.0, 1.3, 0.0)
MARKERS at stadium_symbol		= NhlNewMarker(wks, "O", 18, 0.0, 0.0, 1.0, 0.7, 90.0)
MARKERS at relict_eye_symbol		= 14
MARKERS at banding_eye_symbol		= MARKERS at tropical_storm_symbol

return(MARKERS)

end

;********************************************************
; Function: convert_timeoffset_to_yyyymmddhh
;
; This function takes an array of the time offset (seconds since basetime)
; and converts it to an array of nicely formatted string yyyymmddhh
;********************************************************
undef("convert_timeoffset_to_yyyymmddhh")
function convert_timeoffset_to_yyyymmddhh(time_offset[*]:numeric)

local utc_date,year,month,day,hour,yyyymmddhh_output,it

begin

opt = 0

if (.not.isatt(time_offset,"units")) then
time_offset at units = "seconds since 1970-01-01 00:00:00.0 UTC"
end if

if (isatt(time_offset,"calendar")) then
opt at calendar = time_offset at calendar
end if

utc_date = cd_calendar(time_offset,opt)

; Store return information into more meaningful variables.
year   = floattointeger(utc_date(:,0))    ; Convert to integer for
month  = floattointeger(utc_date(:,1))    ; use sprinti
day    = floattointeger(utc_date(:,2))
hour   = floattointeger(utc_date(:,3))

; Write out date string in the format yyyymmddhh.
yyyymmddhh_output = new(dimsizes(year),"string",string_FillValue)

do it = 0, dimsizes(year)-1
yyyymmddhh_output(it) = sprinti("%0.4i",year(it)) + sprinti("%0.2i",month(it)) + sprinti("%0.2i",day(it)) + sprinti("%0.2i",hour(it))
end do

return(yyyymmddhh_output)

end

;********************************************************
; Function: convert_timeoffset_to_yyyymmddhhmm
;
; This function takes an array of the time offset (seconds since basetime)
; and converts it to an array of nicely formatted string yyyymmddhhmm
;********************************************************
undef("convert_timeoffset_to_yyyymmddhhmm")
function convert_timeoffset_to_yyyymmddhhmm(time_offset[*]:integer)

local utc_date,year,month,day,hour,minute,yyyymmddhhmm_output,it

begin

opt = 0

if (.not.isatt(time_offset,"units")) then
time_offset at units = "seconds since 1970-01-01 00:00:00.0 UTC"
end if

if (isatt(time_offset,"calendar")) then
opt at calendar = time_offset at calendar
end if

utc_date = cd_calendar(time_offset,opt)

; Store return information into more meaningful variables.
year   = floattointeger(utc_date(:,0))    ; Convert to integer for
month  = floattointeger(utc_date(:,1))    ; use sprinti
day    = floattointeger(utc_date(:,2))
hour   = floattointeger(utc_date(:,3))
minute = floattointeger(utc_date(:,4))

; Write out date string in the format yyyymmddhhmm.
yyyymmddhhmm_output = new(dimsizes(year),"string",string_FillValue)

do it = 0, dimsizes(year)-1
yyyymmddhhmm_output(it) = sprinti("%0.4i",year(it)) + sprinti("%0.2i",month(it)) + sprinti("%0.2i",day(it)) + sprinti("%0.2i",hour(it)) + sprinti("%0.2i",minute(it))
end do

return(yyyymmddhhmm_output)
end
-------------- next part --------------
;********************************************************
; spline_interpolation.ncl
;
; The purpose of this program is to test different methods
; of constructing an interpolating cubic spline to a set
; of lat/lon points on a semi-regular time grid.
;
;********************************************************
;********************************************************

; Set default _FillValues that will be used throughout the program
integer_FillValue   = new(1,"integer",-9999)       ; default is -999
float_FillValue     = new(1,"float",-9999)        ; default is -999
double_FillValue    = new(1,"double",-9999)      ; default is -999
string_FillValue    = new(1,"string","")          ; default is "missing"
character_FillValue = new(1,"character",nullChar) ; default is '\0'
logical_FillValue   = new(1,"logical",_Missing)    ; default is _Missing    ; the default _FillValue used to be -1, but in NCL v5.1.0, it has changed to Missing (set by the keyword _Missing)
; logical_FillValue   = new(1,"logical",False)      ; It doesn't work to make False the global _FillValue for logicals since code breaks all over the place

set_default_fillvalue("integer", integer_FillValue)
set_default_fillvalue("float", float_FillValue)
set_default_fillvalue("double", double_FillValue)
set_default_fillvalue("string", string_FillValue)
set_default_fillvalue("character", character_FillValue)
set_default_fillvalue("logical", logical_FillValue)	; this is the new way of setting the default _FillValue throughout the program -- in 2008, this wasn't used or was commented out

;********************************************************

begin

; define input for this example using the actual Best Track values for Hurricane Katrina (2005)
input_time = (/1124820000,1124841600,1124863200,1124884800,1124906400,1124928000, \
1124949600,1124971200,1124992800,1125009000,1125014400,1125036000, \
1125057600,1125079200,1125100800,1125122400,1125144000,1125165600, \
1125187200,1125208800,1125230400,1125252000,1125273600,1125295200, \
1125313800,1125316800,1125326700,1125338400,1125360000,1125381600, \
1125403200,1125424800,1125446400,1125468000/)

input_lat  = (/23.1,23.4,23.8,24.5,25.4,26,26.1,26.2,26.2,26,25.9,25.4,25.1,24.9, \
24.6,24.4,24.4,24.5,24.8,25.2,25.7,26.3,27.2,28.2,29.3,29.5,30.2,  \
31.1,32.6,34.1,35.6,37,38.6,40.1/)

input_lon  = (/-75.1,-75.7,-76.2,-76.5,-76.9,-77.7,-78.4,-79,-79.6,-80.1,-80.3,   \
-81.3,-82,-82.6,-83.3,-84,-84.7,-85.3,-85.9,-86.7,-87.7,-88.6,     \
-89.2,-89.6,-89.6,-89.6,-89.6,-89.6,-89.1,-88.6,-88,-87,-85.3,-82.9/)

output_time = (/1124820000,1124823600,1124827200,1124830800,1124834400,1124838000, \
1124841600,1124845200,1124848800,1124852400,1124856000,1124859600, \
1124863200,1124866800,1124870400,1124874000,1124877600,1124881200, \
1124884800,1124888400,1124892000,1124895600,1124899200,1124902800, \
1124906400,1124910000,1124913600,1124917200,1124920800,1124924400, \
1124928000,1124931600,1124935200,1124938800,1124942400,1124946000, \
1124949600,1124953200,1124956800,1124960400,1124964000,1124967600, \
1124971200,1124974800,1124978400,1124982000,1124985600,1124989200, \
1124992800,1124996400,1125000000,1125003600,1125007200,1125009000, \
1125010800,1125014400,1125018000,1125021600,1125025200,1125028800, \
1125032400,1125036000,1125039600,1125043200,1125046800,1125050400, \
1125054000,1125057600,1125061200,1125064800,1125068400,1125072000, \
1125075600,1125079200,1125082800,1125086400,1125090000,1125093600, \
1125097200,1125100800,1125104400,1125108000,1125111600,1125115200, \
1125118800,1125122400,1125126000,1125129600,1125133200,1125136800, \
1125140400,1125144000,1125147600,1125151200,1125154800,1125158400, \
1125162000,1125165600,1125169200,1125172800,1125176400,1125180000, \
1125183600,1125187200,1125190800,1125194400,1125198000,1125201600, \
1125205200,1125208800,1125212400,1125216000,1125219600,1125223200, \
1125226800,1125230400,1125234000,1125237600,1125241200,1125244800, \
1125248400,1125252000,1125255600,1125259200,1125262800,1125266400, \
1125270000,1125273600,1125277200,1125280800,1125284400,1125288000, \
1125291600,1125295200,1125298800,1125302400,1125306000,1125309600, \
1125313200,1125313800,1125316800,1125320400,1125324000,1125326700, \
1125327600,1125331200,1125334800,1125338400,1125342000,1125345600, \
1125349200,1125352800,1125356400,1125360000,1125363600,1125367200, \
1125370800,1125374400,1125378000,1125381600,1125385200,1125388800, \
1125392400,1125396000,1125399600,1125403200,1125406800,1125410400, \
1125414000,1125417600,1125421200,1125424800,1125428400,1125432000, \
1125435600,1125439200,1125442800,1125446400,1125450000,1125453600, \
1125457200,1125460800,1125464400,1125468000/)

; rename to function input vars for convenience
xi = input_lon
yi = input_lat

; regularly-spaced output
;  noutput_points = 101
;  t = fspan(0,1,noutput_points)

; actual desired output grid
noutput_points = dimsizes(output_time)
t  = 1.*(output_time - output_time(0))   /  (output_time(noutput_points - 1) - output_time(0))

xo = new(noutput_points,"float",-9999)
yo = new(noutput_points,"float",-9999)

; Control parameters that apply to ftkurv are: sig, sl1, sln, sf1.
; Given a sequence of input points ( (xi(0),yi(0)), ... , (xi(npts-1),yi(npts-1))), the interpolated curve is parameterized
; by mapping points in the interval [0.,1.] onto the interpolated curve. The resulting curve has a parametric representation
; both of whose components are splines under tension and functions of the polygonal arc length. The value 0. is mapped onto
; (xi(0),yi(0)) and the value 1. is mapped onto (xi(mpts-1),yi(mpts-1)). Values outside the interval [0.,1.] are mapped to
; extrapolated values.
;
; The value for the parameter sig specifies the tension factor. Values near zero result in a cubic spline; large values
; (e.g. 50) result in nearly a polygonal line. A typical value is 1. (the default).
;
; The value for parameter sl1 is in radians and contains the slope at (xi(0),yi(0)). The angle is measured counter-clockwise
; from the X axis and the positive sense of the curve is assumed to be that moving from point 0 to point npts-1. A value for
; sl1 may be omitted as indicated by the switch sf1.
;
; The value for parameter sln is in radians and contains the slope at (xi(npts-1),yi(npts-1)). The angle is measured
; counter-clockwise from the X axis and the positive sense of the curve is assumed to be that moving from point 0 to point
; npts-1. A value for sln may be omitted as indicated by the switch sf1.

; The value of sf1 controls whether to use the values for sl1 and sln, or compute those values internally. Specifically, sf1
; = 0 if sl1 and sln are user-specified.
; = 1 if sl1 is user-specified, but sln is internally calculated.
; = 2 if sln is user-specified, but sl1 is internally calculated.
; = 3 if sl1 and sln are internally calculated.
;
; You can extrapolate values with ftkurv (that is, calculate interpolated values for abscissae outside of the domain of the input),
; but these values are, in general, unreliable.

; Modify the control factors:
ftsetp("sig",0.1)    ; tension factor: values near zero give a cubic spline, large values (e.g.., 50) result in a near-polynomial line. A typical value is 1 (the default).
;  ftsetp("sl1",1.)    ; specify the slope of the curve at the first point
;  ftsetp("sln",-1.)   ; specify the slope of the curve at the last point
ftsetp("sf1",3)     ; the value of sf1 controls whether to use the values for sl1 and sln, or compute those values
;                       internally. Specifically, sf1
;                           = 0 if sl1 and sln are user-specified.
;                           = 1 if sl1 is user-specified, but sln is internally calculated.
;                           = 2 if sln is user-specified, but sl1 is internally calculated.
;                           = 3 if sl1 and sln are internally calculated - By default the slopes at the end points are computed internally.

; interpolate the best track points to a new output grid using an interpolatory parametric cubic spline fitting function (not an approximating spline)
ftkurv(xi,yi,t,xo,yo)

; print the coordinates of the output
print("xo = " + xo + "   yo = " + yo)

; now validate the interpolated spline values at the Best Track points
print("")
print("yyyymmddhhmm      input lat   output lat    input lon    ouput lon    distance at input points")
print("============      =========   ==========    =========    =========    ========================")

do it = 0, dimsizes(input_time) - 1

indc = ind(output_time .eq. input_time(it))
diff = sqrt((yo(indc) - input_lat(it))^2. + (xo(indc) - input_lon(it))^2.)

print(convert_timeoffset_to_yyyymmddhhmm(input_time(it)) + \
"      " + sprintf("%5.3f",input_lat(it)) + \
"      " + sprintf("%5.3f",yo(indc)) + \
"        " + sprintf("%7.3f",input_lon(it)) + \
"      " + sprintf("%7.3f",xo(indc)) + \
"      " + diff)

end do

;********************************
; create plot
;********************************
plot_filename = "track_spline"

wks = gsn_open_wks("png",plot_filename)  	; Open a postscript workstation.

domain_limits = "florida_closeup"
;  domain_limits = "full"

; set the domain
if (domain_limits .eq. "florida_closeup") then
minlat = 24.
maxlat = 28.
minlon = -82.
maxlon = -78.
else
minlat = min(yi)
maxlat = max(yi)
minlon = min(xi)
maxlon = max(xi)
end if

lat_range = 6.
lon_range = 6.
lon_center = -80.

; Create arrays that only hold the points for the 00z locations
input_time_mm = str_get_cols(convert_timeoffset_to_yyyymmddhh(input_time),4,5)
input_time_dd = str_get_cols(convert_timeoffset_to_yyyymmddhh(input_time),6,7)
input_time_hh = str_get_cols(convert_timeoffset_to_yyyymmddhh(input_time),8,9)

lat0 = mask(yi,stringtointeger(input_time_hh),0)	; note that we have to use stringtointeger because the type is unknown when passing attribute arrays

; set some plotting resources
res = True

res at gsnDraw          = False   ; so we can add poly stuff
res at gsnFrame         = False   ; do not advance frame

res at gsnPaperOrientation = "portrait"

res at gsnMaximize	= True    ; Maximize plot in frame; aspect ratio will be preserved.
res at tiMainString      = "Katrina (2005)"
res at gsnCenterString   = "Comparison of Splining"

res at mpDataBaseVersion  = "Ncarg4_1"     ; Alias 'MediumRes'
res at mpFillOn              = True   		; False to turn off gray continents
res at mpOutlineOn           = True    		; turn on continental outline

;  res at mpPerimOn = True	; turn on the map perimeter

; the following set the resources for the lat/lon grid (but not the tickmarks)
res at mpGridAndLimbOn        = "True"
res at mpGridAndLimbDrawOrder = "Draw"
res at mpGridLineColor        = "grey37"

res at mpLimitMode      = "LatLon"

res at mpMinLatF        = minlat
res at mpMaxLatF        = maxlat
res at mpMinLonF        = minlon
res at mpMaxLonF        = maxlon

res at mpGridLatSpacingF  = 5	; set the spacing for the grid
res at mpGridLonSpacingF  = 5

res at gsnMajorLatSpacing = res at mpGridLatSpacingF	; set the spacing for the labels
res at gsnMajorLonSpacing = res at mpGridLonSpacingF

;  res at gsnMinorLatSpacing = 1.0	; currently this function doesn't support minor tickmarks
;  res at gsnMinorLonSpacing = 1.0

res at tmXTOn = False
res at tmYROn = False

res at tmXTLabelsOn = False
res at tmYRLabelsOn = False

;  opt = True

;  opt at SpanInternationalDateline = True
;  opt at OverrideTickmarkLabelSpacing = True
;  opt at MaxNumberLatTickmarkLabels = res at mpGridLatSpacingF
;  opt at MaxNumberLonTickmarkLabels = res at mpGridLonSpacingF

;  set_nice_map_tickmarks(res,opt)	; set up some nice tickmarks where we can ensure that the spacing for lat and lon labels are the same

; We want to draw the lat/lon cross point markers over ocean only, so we need to:
;
; - Create a map plot with filled land and transparent ocean.
; - Create a map plot with filled ocean and transparent land.
; - Attach the land plot to the ocean plot so that the ocean
;    plot will get drawn first.
; - Attach the grid markers to the ocean plot. This means that the
;   land plot will cover the grid markers.
; - Attach the other polymarkers and polylines to the land plot so
;   that they will get drawn overtop both the land and ocean plots
;
; Note: we have to use "gsn_add_annotation" to attach one map
; plot to another, because you can't use the "overlay" procedure
; to overlay two map plots.

;******************************************
; Create ocean plot by setting water fill *
; resources and making land transparent   *
;******************************************
res at mpDataSetName 	     = "Earth..4"    ; was "Earth..4"

res at mpLandFillColor        = "Transparent" ; was "GoldenRod1" or "Tan"
res at mpOceanFillColor       = "SlateGray2" ; was "LightBlue1"
res at mpInlandWaterFillColor = "SlateGray2" ; was "LightBlue1"	; was "PaleTurquoise3"

ocean_plot = gsn_csm_map(wks,res)

;**********************************************************************
; Plot little dots at the integer crossings of latitude and longitude *
;**********************************************************************
; Set up resources for polymarkers.

resp               = True

resp at gsMarkerIndex = 16          ; Filled dots
resp at gsMarkerColor = "grey37"

resp at gsMarkerSizeF = 0.001

; put in a check here to make sure the lon's go in reverse if we cross the dateline

ilatmin = floattoint(floor(res at mpMinLatF - 0.5))
ilatmax = floattoint(ceil(res at mpMaxLatF + 0.5))
ilonmin = floattoint(floor(res at mpMinLonF - 0.5))
ilonmax = floattoint(ceil(res at mpMaxLonF + 0.5))

ilat_range = ilatmax - ilatmin + 1
ilon_range = ilonmax - ilonmin + 1

dum_dots = new((ilat_range*ilon_range),"graphic")

ii = 0
do ilon = ilonmin, ilonmax
do ilat = ilatmin, ilatmax
if ( ((ilat % res at mpGridLatSpacingF) .ne. 0) .and. ((ilon % res at mpGridLatSpacingF) .ne. 0) ) then
ii = ii + 1
end if
end do
end do

; Here's another way of doing it, but this would put them everywhere and sometimes these show up on the lines

; Create a 2D array of lat/lon values, but make this a 1D array
; so we can pass them to gsn_add_polymarker.

;  lat_grid = ispan(minlat,maxlat,1)
;  lon_grid = ispan(minlon,maxlon,1)

;  nlat_grid = dimsizes(lat_grid)
;  nlon_grid = dimsizes(lon_grid)

;  lat2d_1d = ndtooned(conform_dims((/nlat_grid,nlon_grid/),lat_grid,0))
;  lon2d_1d = ndtooned(conform_dims((/nlat_grid,nlon_grid/),lon_grid,1))

; Attach polymarkers to ocean plot.
;  dum = gsn_add_polymarker(wks,ocean_plot,lon2d_1d,lat2d_1d,resp)

;********************************************
; Create land plot by setting land fill     *
; resources and making water transparent    *
; we also plot the national (and internal)  *
; boundaries so that the U.S. counties will *
; appear.                                   *
;********************************************
res at mpLandFillColor        	= "Wheat3"	    ; was "GoldenRod1" or "Tan"
res at mpInlandWaterFillColor 	= "Transparent"
res at mpOceanFillColor       	= "Transparent"

res at mpOutlineBoundarySets  	= "AllBoundaries"   ; "GeophysicalAndUSStates"  "AllBoundaries"

res at mpNationalLineColor	= "grey37"
res at mpNationalLineThicknessF  = 0.5

res at mpOutlineDrawOrder 	= "Draw"
res at mpLabelDrawOrder   	= "Draw"
res at mpPerimDrawOrder   	= "Draw"
;  res at mpFillDrawOrder 		= "PreDraw" ; Default is "PostDraw" - we're going to put some dots on the map, but we only want them to show up on the water areas
;  res at mpGridAndLimbDrawOrder 	= "Draw"
;  res at tfPolyDrawOrder    	= "PostDraw" ; Draw all primitives on top

land_and_counties_plot = gsn_csm_map(wks,res)

; Add the land plot as an annotation of the ocean plot.
; This means the ocean plot (and all attached stuff) is drawn first.

amres = True

;********************************************************************
; Create plot with state boundaries drawn in a thicker, darker gray *
;********************************************************************
res at mpFillOn                    = False

res at mpNationalLineColor         = -1	; transparent
res at mpUSStateLineColor          = "grey27"
res at mpUSStateLineThicknessF     = 1.5

state_plot = gsn_csm_map(wks,res)
annoid_state_plot = gsn_add_annotation(land_and_counties_plot,state_plot,amres)   ; Add the state plot as an annotation of the ocean plot.

;***************************************
; Create plot with national boundaries *
;***************************************
res at mpOutlineBoundarySets = "National"	; this first way will do all national boundaries
res at mpNationalLineColor    = "Black"
res at mpNationalLineThicknessF    = 2.0

;  res at mpNationalLineColor            = -1	; transparent	; this second way would do just the U.S. national boundaries
;  res at mpUSnationalLineColor          = "Black"
;  res at mpUSnationalLineThicknessF     = 1.5

res at mpGeophysicalLineColor = "grey37"
res at mpGeophysicalLineThicknessF = 1.0

national_plot = gsn_csm_map(wks,res)
annoid_national_plot = gsn_add_annotation(state_plot,national_plot,amres)   ; Add the national plot as an annotation of the ocean plot.

;***************************************************************
; Now add the Best Track elements to the land plot             *
;***************************************************************
nbest = dimsizes(input_time)

; first, set the line styles for each segment based on intensity and storm type
line_colors        		= new(nbest,"string")
line_thicknesses   		= new(nbest,"float")
line_dash_patterns 		= new(nbest,"integer")
line_dash_segment_lengths 	= new(nbest,"float")

; first add the BT track
do k = 0, nbest - 2

; set some default values for this line segment
line_colors(k)	       = "Black"
line_thicknesses(k)       = 5
line_dash_patterns(k)     = 0
line_dash_segment_lengths = 0.15

if (.not.ismissing(yi(k))) then	; don't do if the initial point is missing
add_line(wks,national_plot,xi(k:k+1),yi(k:k+1),line_colors(k),line_thicknesses(k),line_dash_patterns(k),line_dash_segment_lengths(k))	; to try the other way of doing things that avoids the warnings
end if  ; not missing data

end do ; end loop over ibest

delete(line_colors)
delete(line_thicknesses)
delete(line_dash_patterns)
delete(line_dash_segment_lengths)

;***************************************************************
; Add the interpolated Best Track elements to the land plot    *
;***************************************************************
nbest_interpolated = dimsizes(output_time)

; first, set the line styles for each segment based on intensity and storm type
line_colors        		= new(nbest_interpolated,"string")
line_thicknesses   		= new(nbest_interpolated,"float")
line_dash_patterns 		= new(nbest_interpolated,"integer")
line_dash_segment_lengths 	= new(nbest_interpolated,"float")

; set some default values for this line segment
line_colors	            = "Red"
line_thicknesses          = 3
line_dash_patterns        = 0
line_dash_segment_lengths = 0.15

; add the interpolated BT track
do k = 0, nbest_interpolated - 2
add_line(wks,national_plot,xo(k:k+1),yo(k:k+1),line_colors(k),line_thicknesses(k),line_dash_patterns(k),line_dash_segment_lengths(k))	; to try the other way of doing things that avoids the warnings
end do

size_factor = 0.9

;********************************************
; Draw black polymarkers at 00Z BT locations
;********************************************
do ii = 0, dimsizes(lon0) - 1
end do

;********************************************
; Draw white polymarkers at 12Z BT locations
;********************************************
do ii = 0, dimsizes(lon12) - 1
end do

;************************************************************************************************************************
; Draw the Chelmow-Willoughby fixes as small yellow circles
;************************************************************************************************************************
do it = 0, dimsizes(xo) - 1
end do

;***********************
; Plot some day labels *
;***********************
txres_day    = True

txres_day at gsnDraw = False  ; True
txres_day at txJust  = "TopRight"

; Specify offsets for the hour labels
xhoffset = -0.25*lon_range/20.0	; normalize the offset according to the domain size
yhoffset = -0.25*lat_range/20.0

; Now draw the day labels at each 00Z marker
txres_day at txFontHeightF = 0.013	; 0.008

dumz = new(nbest,"graphic")

month_names = (/"","Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"/)

do ib = 0, nbest - 1

if ( (ib .eq. 0) .or. (ib .eq. nbest - 1) .or. (.not.ismissing(day12(ib)) .and. (stringtointeger(day12(ib)) .eq. 1)) ) then
day_label = month_names(stringtointeger(input_time_mm(ib))) + " " + day12(ib)
else
day_label = day12(ib)
end if

dumz(ib) = gsn_add_text(wks,national_plot,day_label,lon12(ib) + xhoffset,lat12(ib) + yhoffset,txres_day)
end do

;*****************************************************************
; create a legend and attach as an annotation to the land plot   *
;*****************************************************************
lgres                    = True

lgres at lgAutoManage       = False

lgres at vpWidthF           = 0.25       ; was 0.08        ; width of legend (NDC)
lgres at vpHeightF          = 0.20       ; was 0.08        ; height of legend (NDC)
;  lgres at lgBottomMarginF    = 0.17     ; was 0.25

lgres at lgPerimFill        = 0          ; Use solid fill (0) instead of the default hollow fill
lgres at lgPerimFillColor   = "White"	; was "Background"
lgres at lgPerimThicknessF  = 1.5

lgres at lgBoxMajorExtentF  = 0.0 ; was 0.1
lgres at lgBoxMinorExtentF  = 0.2 ; was 0.2	; controls how wide the box holding the legend items (lines and markers) can be in relation to legend

;  lgres at lgBoxBackground    = "PaleTurquoise3"

;  lgres at lgLabelFont             = "helvetica"
lgres at lgLabelFontHeightF      = 0.04
lgres at lgLabelFontAspectF      = 1.2		     ; values larger than 1.0 make the font skinny and tall, values less than 1.0 make the font fat and short
;  lgres at lgLabelConstantSpacingF = 0.0	             ; amount of extra space to add between characters - values less than 0.0 not allowed
lgres at lgLabelJust		= "CenterLeft"	; default is "CenterCenter"
lgres at lgLabelOffsetF		= 0.04			; controls the amount of space (along the minor axis) between the legend lines and the legend text labels; default is 0.02 and the direction depends on the value of lgLabelJust

lgres at lgMonoItemType        = False                 ; indicates that we wish to set the item types individually
lgres at lgMonoMarkerIndex     = False
lgres at lgMonoLineThickness   = False
lgres at lgMonoMarkerThickness = False
lgres at lgMonoMarkerSize      = False

lgres at lgLineDashSegLenF  = 0.1

;                              BT (interpolated)    12Z position   00Z position    BT
lgres at lgItemTypes        = (/"MarkLines"         , "Markers"    ,"Markers"     ,"Lines"/)
lgres at lgMarkerIndexes    = (/ 1                  ,  4           , 16           , 0     /)
lgres at lgLineThicknesses  = (/ 3.0                ,  0.1         , 0.1          , 5.0   /)
lgres at lgMarkerColors     = (/"Yellow"            , "White"      ,"Black"       ,"Black"/)
lgres at lgMarkerSizes      = (/ 0.015              ,  0.015       , 0.015        , 0.015 /)*size_factor
lgres at lgLineColors       = (/"Red"               , "White"      ,"White"       ,"Black"/) ; colors for legend lines
lgres at lgDashIndexes      = (/ 0                  ,  0           , 0            , 0     /) ; dash indexes

lgres at lgItemCount        = dimsizes(lgres at lgItemTypes)

legend_labels = (/"Best Track (interpolated)",  \
"Position/date at 1200 UTC ", \
"Position at 0000 UTC ",      \
"Best Track"/)

; NOTE: the number of legend items must be correct otherwise a segmentation fault occurs

legend = gsn_create_legend(wks,lgres at lgItemCount,legend_labels,lgres)

; attach the legend to the land plot as an annotation
amres1 = True

amres1 at amZone		  = 0		; This zone starts off in the exact middle of the viewport	; comment this line out to get it inside the plot
amres1 at amJust           = "TopRight"	; was "TopRight"
amres1 at amParallelPosF   =  0.49	; was 0.49
amres1 at amOrthogonalPosF = -0.49	; was -0.49

;****************************************
; Attach a day label onto the land plot *
;****************************************
txres2 = True

txres2 at txPerimOn             = False
txres2 at txFontHeightF         = 0.008
txres2 at txJust                = "CenterCenter"

number_text = gsn_create_text(wks,"dd",txres2)

; Add the day label to the legend as an annotation
amres2 = True

amres2 at amZone		  = 0		; This zone starts off in the exact middle of the viewport	; comment this line out to get it inside the plot
amres2 at amJust           = "TopRight"
amres2 at amParallelPosF   =  0.33	; was 0.33
amres2 at amOrthogonalPosF = -0.235	; was -0.245

;********************************
; Plot a title on the land plot *
;********************************
txres_title = True

txres_title at txFont               = 22	; Helvetica bold
txres_title at txFontHeightF        = 0.025
;  txres_title at txConstantSpacingF   = 1.0
;  txres_title at txFontAspectF        = 1.2
;  txres_title at txFontThicknessF     = 1.6
txres_title at gsnDraw              = False

txres_title at txPerimOn             = True
txres_title at txPerimColor          = "Black"
txres_title at txPerimThicknessF     = 1.5
txres_title at txPerimSpaceF         = 0.4
txres_title at txBackgroundFillColor = "White"

title_text = gsn_create_text(wks,"Best Track and an interpolated spline for~C~       Katrina (2005)",txres_title)

; now add the title text as an annotation
amres_title                  = True

amres_title at amParallelPosF   = -0.13
amres_title at amOrthogonalPosF = -0.49
amres_title at amJust           = "TopCenter"
; annoid_title = gsn_add_annotation(national_plot,title_text,amres_title)

;  drawNDCGrid(wks)

draw(ocean_plot)    ; Note markers get drawn under land.
frame(wks)

end
-------------- next part --------------
A non-text attachment was scrubbed...
Name: track_spline.png
Type: image/png
Size: 66677 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20151223/9e97b0c4/attachment.png
```