[Re-opened] Interpolation scheme to be used in WAFS SIGWX Forecast

[blchoy]

The document Representing WAFS SIGWX Data in BUFR mentioned that "An appropriate smoothing technique such as a cubic spline should be used to smooth the jet-stream curve plotted between the core points." This also applies to the horizontal boundary of the CAT and cloud areas.

In the IWXXM implementation of the WAFS SIGWX Forecast, it is anticipated to carry forward this idea in the encoding:

<!-- A Jet Stream Object -->
<iwxxm:feature>
    <iwxxm:MeteorologicalFeature gml:id="uuid.8c3d3467-e398-4627-ab23-f9e57c2215b5">
        <gml:identifier codeSpace="http://wafs/sigwxobj">eb324455-236c-46f8-b80d-c6d42c3a634c</gml:identifier>
        <iwxxm:phenomenon xlink:href="http://codes.wmo.int/49-2/MeteorologicalFeature/JETSTREAM"/>
        <iwxxm:phenomenonGeometry>
            <gml:Curve gml:id="uuid.7b5a197a-d9bf-4f3b-ba36-701fe019b6d2" srsDimension="2" axisLabels="Lat Long" srsName="http://www.opengis.net/def/ers/EPSG/0/4326">
                <gml:segments>
                    <gml:CubicSpline>
                        <gml:posList>35.4 -104.7 36.1 -101.6 36.6 -99.6 36.9 -97.7 37.4 -94.9 38.3 -92.3 39.2 -90.3 40.7 -88.1 41.5 -86.3 42.0 -82.7 41.9 -78.9</gml:posList>
                        <gml:vectorAtStart>1.0 2.0</gml:vectorAtStart>
                        <gml:vectorAtEnd>2.0 3.0</gml:vectorAtEnd>
                    </gml:CubicSpline>
                </gml:segments>
            </gml:Curve>
        </iwxxm:phenomenonGeometry>
        ...

@borisburger in his e-mail indicated difficulties in proper implementation of cubic splines and propose the use of Bezier. @GAndrsn after studying Boris' suggestion had identify some other issues with Bezier.

Summarizing the discussion so far, I think the key issues on this topic is as follow:

1. For WAFCs a smooth interpolation scheme is needed for both the boundaries of areas and the jet core. it is essential for the interpolated curve to pass through the points indicated in the gml:posList, especially in the latter case as the direction of the jet stream is derived from the interpolated curve.
2. For software developers it is difficult, if not impossible, to properly implement cubic spline interpolation, especially if the application requires realtime adjustment to the curves like resizing of the polygon.
3. The CRS involved (i.e. EPSG:4326) further complicates this issue.

I hope I have given a correct account to the issue. Graeme and Boris feel free to express your views. Others are most welcome to join the discussion.

[borisburger]

I have read the response from Graeme and I completely agree with the view that forecasters do not ever think about the particular interpolation technique used by the underlying software tool. However an advantage of using Bezier cubic splines is that the software will directly encode what user had on screen without any room for interpretation - so the recipients of the data will draw the same shape the forecaster saw on screen.

This is because all major drawing APIs on Linux & Windows like Microsoft GDI, GDI+, Direct 2D/Direct Draw, Cairo, Skia, OpenGL etc. support cubic Bezier splines out of the box. Vendors then have a very easy job of plotting the curve with a single API call.

My main concern with the proposed gml:CubicSpline is not the difficulty in implementation - this kind of cubic spline is uncommon, but not overly difficult to implement. Why do I say that it is uncommon?

The real issue is with the underlying GM_CubicSpline in ISO 19107 Geographic Information - Spatial Schema where a very unfortunate name/terminology was chosen - using the general term "cubic spline" to actually denote a so called "C2 continuous natural cubic spline" (a standard textbook name). Quoting from ISO 19107 Geographic Information - Spatial Schema, section 7.13.6 Interface CubicSpline:

the curve must be C2, i.e. have a continuous first and second derivative at all points ... The control parameters record must contain derivativeAtStartvectorAtStart, and derivativeAtEndvectorAtEnd which are the tangent vectors at controlPoint[1] and controlPoint [n] ...

The above is a very particular and unique type of C2 continuous cubic splines (C2 = continuous 1st and 2nd derivatives). The general term "cubic spline" would normally denote a whole class of cardinal, Catmull-Rom, Kochanek-Bartels, Hermitean and other types of polynomial cubic splines. What is denoted by gml:CubicSpline is very confusingly not this whole class, but a very specific C2 continuous type. Please note in the quote from ISO 19107 that even with gml:CubicSpline you still need to specify derivatives at the curve endpoints! (but it is true you do not need to specify the derivatives at the inner control points where the segments join, because these can be fully determined from the C2 continuity requirement - but for the outer endpoints it is still necessary).

A good overview article is https://en.wikipedia.org/wiki/Spline_(mathematics) which mentions many cubic spline types and also explains the specific natural C2 cubic splines as used in GML. Also https://en.wikipedia.org/wiki/Cubic_Hermite_spline

The issue I see with this is that not everyone studied computational geometry or computer graphics. I am afraid some vendors due to the misleading (overly broad) term gml:CubicSpline will not realise, that in GML (based on ISO 19107 GM_CubicSpline) this is a very special type of "C2 continuous natural cubic splines" and instead they will just use the same spline algorithm they have been using for the past 20 years - not displaying what was originally intended.

A brief summary about gml:CubicSpline (GM_CubicSpline):

1. The curve passes through all its control points and is smooth in mathematical sense = both 1st and 2nd derivative are continous along the entirety of the curve
2. Derivatives at outer endpoints of open curves (e.g. jetstreams) still need to be provided for every gml:CubicSpline
3. Implementations are unfortunately not readily available in major drawing APIs.
4. Implementation requires solving a 3-diagonal system of linear equations to figure out the 1st and 2nd derivatives at control points, and then convert to a Bezier to satisfy the graphics APIs that expect Beziers. The number of equations in the system is equal to the number of spline control points (100 control points leads to a system of 100 equations with 100 unknowns). It is however a textbook algorithm which can be found online and is taught in geometry/CG classes at universities, however it is not clear if people will realise this, because:
5. I have never seen C2 continuous natural cubic splines offered in graphics software like Autodesk AutoCAD, Fusion 360, or vector drawing software like Inkscape. I do not own a license for Adobe Illustrator, so I do not know if this type of C2 spline is offered there. The most commonly quoted reason for graphics software not offering C2 natural cubics is their "lack of local control" = if you move a single control point, you need to recalculate the entire 3-diagonal system of equations from scratch and the entirety of the curve changes striving to satisfy the C2 continuity. This property makes the C2 natural cubic splines a bit too unpredictable for graphics artists.

I am not fully opposed to using gml:CubicSpline. It is an intriguing concept and in GML besides gml:Bezier and gml:LinearRing there is no other reasonable choice. I fully recognize the advantage of not having to specify the additional control points like in Bezier to determine derivatives (however at outer endpoints they still need to be provided!). And the compact representation offered by gml:CubicSpline where the spline passes through all the specified control points is also a useful property.

I am just saying that the choice of gml:CubicSpline has a potential downside, that implementors will get confused by the overly generic "cubic spline" term and fall back on their existing algorithms, instead of going the more difficult route of solving the 3-diagonal system of equations for each curve. For those that will take the extra step this will make the software less efficient because of the additional computational complexity as opposed to Bezier or LinearRing. But I guess the overall impact of the extra math will not be so large considering all the other computations that need to be performed to plot the data on screen.

Until now I was trying to represent the broader community of software vendors, but I would appreciate if members of other companies wrote their point of view. Specifically for IBL, for manual drawing of SIGWX-like charts by the forecaster (in our drawing tools at least) we will still have to fall back on just C1 continuous cubic splines. Our company is currently using Cardinal cubic splines for SIGWX charts (https://en.wikipedia.org/wiki/Cubic_Hermite_spline) with c=0 tension parameter and a specific knot sequence. It is a simpler type of a cubic spline which makes approximation easier (approximation of jets or clouds drawn freehand on a tablet, or approximation of isocontours based on e.g. NWP MSLP pressure or geopotential). For the C2 natural cubics we do not know how to implement such numerical approximation algorithms that fit a 2D freehand curve. So if SIGWX switches over to C2 splines we will likely need to support also C2 continuous gml:CubicSpline and then approximate it with the Cardinal C1 cubic spline if users will want to continue tweaking the SIGWX charts for their area of interest (we have many users who base their mid-level SIGWX charts on the WAFC SIGWX data by importing the jets/clouds/CATs into graphical editor and moving control points around).

Kind regards,
Boris

[blchoy]

I see. That means we need to skip using gml:CubicSpline and define our own representation of a curve with points along it and interpolation scheme like C1 cubic spline?

[GAndrsn]

I've had a read through all of the info and attachments in the discussion so far, and I think I have a broad understanding of the issues. I never studied geometry in depth, however, so I'm only gradually wrapping my head around Bezier curves and all the different types of splines.

I can see the advantages of curves defined as Bezier curves being totally unambiguous, whereas points described simply as defining a "spline" could mean a variety of curves. The only reason I have used cubic splines in the way I have so far is that this is what I inherited when I took on the SigWx code, and I didn't have much of an understanding of the complexities of curves at the time (and still don't fully, to be honest!).

I don't object to trying to implement using Bezier curves. I googled how to convert splines to Bezier curves earlier, and funnily enough came across the same link in Jan's email which I didn't look at until later. I can have a go and seeing what effect implementing this has, although it may take me a little while because, as I say, my knowledge of curve geometry is currently relatively limited.

One of my concerns would be making sure that all of the eventual users of the data are absolutely clear that going from manually-drawn objects in BUFR data to algorithmically-drawn objects in IWXXM that they can't just use their old spline functions, as the Bezier control points will presumably be outside of the hazard area perimeter, rather than on the perimeter itself. I'll discuss this with my colleagues who deal with the users of this sort of data.

[blchoy]

I think we do not need to limit ourselves to gml:CubicSpline and gml:BSpline but other interpolation schemes that can be described by GML, as long as it is the kind of smooth curve that is both fit-for-purpose (meeting WAFCs' expectation) and achievable (something that developers can do).

In case there is an agreed discussion time I can arrange a Webex meeting. I normally available 12-16UTC during weekdays (but not for today 28 May, as I have another discussion with Anna).

[GAndrsn]

@borisburger I don't have an in-depth understanding of the differences between types of spline. My assumption has been that the scipy.interpolate functions produce a B-Spline, which is a type of spline, and that a third-order B-Spline could be described as a cubic spline. As such when we've represented the data in IWXXM we've used the cubic spline type, but that has been based on advice from Choy. As a meteorologist the details here really aren't my area of expertise!

By default, the B-Splines produced by the scipy functions don't interpolate the input points, but the argument s (a "smoothing condition") can be set to 0, in which case it does (see details here). There are references at the bottom of the page by P. Dierckx that might clarify what's going on, but they require payment to access, and I doubt I would understand them anyway.

I am happy to join a Webex next week, Tues-Thurs (I am not in next Mon or Fri). My only meetings so far are Tues and Weds, 13-14UTC, but any other time is fine.

[borisburger]

I am really sorry for being slow to respond, this week is a bit crazy for me due to a milestone presentation we are preparing for a customer. I am available:

• Wednesday 07-13 UTC and 14-16 UTC.
• Thursday 12-16 UTC (at least I hope the presentation we will be doing for the customer will be over by 12UTC)

[blchoy]

So the commonly available slots are Wed 14-16UTC and Thu 13-16UTC.

May I suggest we try Wed 1430-1600UTC first? I will send a Webex invitation now.

[borisburger]

I confirm my availability Wed 14:30-1600 UTC.

[blchoy]

Thank you for Graeme and Boris joining the discussion yesterday. I would like to highlight the outcomes and follow up actions below. Feel free to comment or modify.

1. The interpolation function in SciPy seems to be able to comply with the C2 requirement of a natural cubic spline which gml:CubicSpline is referring. A number of documents were referenced during the discussion and if Boris could drop a few lines what they are Choy may be able to invite some of his colleagues to confirm.
2. It may be worth trying to reproduce the curves Graeme produced with a different software to confirm this could be done. Choy will see if there is any volunteer to conduct a test.
3. The existing WAFC's guidance on rendering may need to beef up to include necessary details to facilitate reproduction of the curves. Graeme will work with Karen on this afterwards.
4. There was a discussion on other possible issues, including "dateline crossing features" and "features close to the north or south pole". It appeared that as individual object could specify its own CRS it may help to better describe their geometry. As this would not affect the IWXXM schemas, Choy and Graeme would re-visit this closer to the time when WAFCs is going releasing data for operation trial.

Finally there is hope for us to safely close this topic. :)

I would also like to forward Boris' comment that what Graeme has been doing with regard to the generation of weather objects is very impressive. I can't agree more with him.

[GAndrsn]

Thank you to Choy for re-opening this discussion, so Boris can potentially contribute as well. I emailed Choy with the following question, which I hope one of you can answer:

Can you clarify for me exactly what the end vectors should represent? At the moment I think I’ve been calculating them as unit vectors in the direction of the start/end of the spline, but I suspect this is incorrect. I have a vague recollection of mention of these vectors representing derivatives of the spline. Could you clarify this for me? I would think that whether the end vectors are correct or incorrect will rapidly become apparent if we try reading the IWXXM data into software that can interpret glm:CubicSpline data.

[borisburger]

I am sorry for not continuing this discussion. The two vectors should really represent the derivatives at start and end. I am attaching an excerpt from ISO 19107:2019 Geographic Information - Spatial Schema:

[GAndrsn]

Great, thanks for the quick response, I will look at how I implement this in my IWXXM writing code.

[GAndrsn]

Using the scipy.interpolate.splev() function with the 'der' keyword set to 1 appears to create a first derivative vector that gives the desired tangent vector at any point along the spline. See output below for interpolation of curve y = -(0.2*(x-34.5))**2 between 0 < x < 70:

Having done a bit more testing, I believe that the magnitude of the derivative vectors is caused by the derivative vector representing (dx/ds, dy/ds), where s is the parameter that varies from 0 to 1 along the length of the spline. As such, for a quadratic curve from 0 to 70, ds will be less than 1/100th of the length of the spline per unit length of the spline, and so the magnitude of the derivative vector will be greater than 100. How this is scaled will need to be considered for the implementation of this type of derivative for lat/lon data.

I can apply this new approach to calculating end vectors to my SigWx IWXXM code, but I won't know I've got it right until the data is read into a proper visualisation software package and I can see how the data is interpreted.