[Python] How to compute the smoothed curve defined by gml:CubicSpline in a WAFS SIGWX Forecast?

I'm studying the gml:CubicSpline smoothing algorithm used in the IWXXM WAFS SIGWX Forecast, based on the WAFS-Example.xml sample.

For example, consider the following GML curve:

enter code here

<gml:Curve gml:id="uuid.012c21ce-6d29-4769-bdca-d042b306fca0" srsDimension="2" axisLabels="Lat Long" srsName="http://www.opengis.net/def/crs/EPSG/0/4326">
<gml:segments>
<gml:CubicSpline>
<gmlosList> -41.623153 -158.734129 -42.921036 -158.310418 -43.614695 -156.565840 -44.448049 -154.972043 -45.933174 -154.748367 -46.733087 -153.171172 -46.985416 -151.164300 -46.884731 -149.037755 -47.995659 -147.690252 -48.199330 -145.711761 -46.904387 -145.018761 -45.627836 -145.977674 -44.393574 -146.977162 -42.959553 -147.449231 -41.878655 -148.831074 -41.959083 -150.789464 -41.053172 -152.365167 -40.289571 -154.037198 -40.159991 -155.942752 -40.378993 -157.783406 -41.623153 -158.734129 </gmlosList>
<gml:vectorAtStart> -35.75874 -5.02593 </gml:vectorAtStart>
<gml:vectorAtEnd> -35.75874 -5.02593 </gml:vectorAtEnd>
</gml:CubicSpline>
</gml:segments>
</gml:Curve>


I used Python's SciPy library to compute the smoothed coordinates, and obtained the following smoothed curve:


python code:

import numpy as np
from scipy.interpolate import CubicSpline
import matplotlib.pyplot as plt

def parse_data_to_list(s):
s_arr = s.split(' ')
s_arr_tuple = [tuple(s_arr[i:i + 2]) for i in range(0, len(s_arr), 2)]
re_list = list()
for i in s_arr_tuple:
re_list.append(list(map(float, i)))
return re_list

list_data = parse_data_to_list('-41.623153 -158.734129 -42.921036 -158.310418 -43.614695 -156.565840 -44.448049 -154.972043 -45.933174 -154.748367 -46.733087 -153.171172 -46.985416 -151.164300 -46.884731 -149.037755 -47.995659 -147.690252 -48.199330 -145.711761 -46.904387 -145.018761 -45.627836 -145.977674 -44.393574 -146.977162 -42.959553 -147.449231 -41.878655 -148.831074 -41.959083 -150.789464 -41.053172 -152.365167 -40.289571 -154.037198 -40.159991 -155.942752 -40.378993 -157.783406 -41.623153 -158.734129')

points = np.array(list_data)

# vectorAtStart 和 vectorAtEnd
vector_at_start = np.array([-35.75874 ,-5.02593]) # dx/dt = 1, dy/dt = 0
vector_at_end = np.array([-35.75874 ,-5.02593]) # dx/dt = 0, dy/dt = -1

# build parameter t
diffs = np.diff(points, axis=0)
segment_lengths = np.linalg.norm(diffs, axis=1)
t = np.concatenate([[0], np.cumsum(segment_lengths)])

# x(t), y(t)
x = points[:, 0]
y = points[:, 1]

# clamped cubic spline
cs_x = CubicSpline(t, x, bc_type=((1, vector_at_start[0]), (1, vector_at_end[0])))
cs_y = CubicSpline(t, y, bc_type=((1, vector_at_start[1]), (1, vector_at_end[1])))
# cs_x = CubicSpline(t, x)
# cs_y = CubicSpline(t, y)

# sample
t_sampled = np.linspace(t[0], t[-1], 200)
x_smooth = cs_x(t_sampled)
y_smooth = cs_y(t_sampled)

# render
plt.figure(figsize=(8, 5))
plt.plot(x, y, 'o', label='Original Points', color='red')
plt.plot(x_smooth, y_smooth, '-', label='CubicSpline with vectorAtStart/End', color='blue')
plt.quiver(x[0], y[0], vector_at_start[0], vector_at_start[1], angles='xy', scale_units='xy', scale=1, color='green', label='vectorAtStart')
plt.quiver(x[-1], y[-1], vector_at_end[0], vector_at_end[1], angles='xy', scale_units='xy', scale=1, color='orange', label='vectorAtEnd')
plt.axis('equal')
plt.grid(True)
plt.legend()
plt.title("GML CubicSpline with vectorAtStart and vectorAtEnd")
plt.xlabel("X")
plt.ylabel("Y")
plt.show()


did I do something wrong?

Continue reading...