ssrfpack unif function
kwagstyl opened this issue · 12 comments
Hi,
Thank you so much for creating this tool. It's a great resource to have in python.
In trying to replicate spherical interpolation with splines under tension from Renka,
I have so far been able to run the nearest, linear and cubic interpolations from my smaller number of spherical sample points to the larger, more densely sampled sphere. The issue I'm running into is that cubic interpolations introduces some massive and implausible over and undershoots.
As I understand it, the way of dealing with this is in ssrfpack is either smoothing or interpolating with tension.
In the original documentation the steps for interpolating with tension are as follows:
(1) Reserve storage as described above and store input parameter values: N, X, Y, Z, W, PLAT, and PLON. (2) Create a Delaunay triangulation of the set of N nodes (subroutine TRMESH). (3) Compute estimated gradients at the nodes (subroutine GRADG or GRADL). (4) Optionally, compute tension factors SIGMA chosen to preserve local shape properties (subroutine GETSIG) or to satisfy arc constraints (subroutine SIG0, SIG1, or SIG2). (5) Compute interpolated values in FF (subroutine UNIF).
I think step (5) is currently not available in stripy._ssrfpack
Was the UNIF subroutine excluded? It's there in ssrfpack.pyf but I think it's commented out.
Do you know how I could add it to stripy?
Thanks,
Konrad
The UNIF routine is in SSRFPACK and so it could be wrapped. Would we need to expose the subsidiary routines as well ? And which routines would we use to obtain the equivalent result at the original (or a different) triangulation. How does this differ from the wrapper for SMSURF ?
@brmather in the first instance ... and if you don't have a clear idea, should we ask Renka what he would recommend ?
SUBROUTINE UNIF (N,X,Y,Z,F,LIST,LPTR,LEND,IFLGS,SIGMA,
. NROW,NI,NJ,PLAT,PLON,IFLGG, GRAD, FF,
. IER)
INTEGER N, LIST(*), LPTR(*), LEND(N), IFLGS, NROW, NI,
. NJ, IFLGG, IER
REAL X(N), Y(N), Z(N), F(N), SIGMA(*), PLAT(NI),
. PLON(NJ), GRAD(3,N), FF(NROW,NJ)
C
C***********************************************************
C
C From SSRFPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 07/25/96
C
C Given a Delaunay triangulation of a set of nodes on the
C unit sphere, along with data values and tension factors
C associated with the triangulation arcs, this routine
C interpolates the data values to a uniform grid for such
C applications as contouring. The interpolant is once con-
C tinuously differentiable. Extrapolation is performed at
C grid points exterior to the triangulation when the nodes
C do not cover the entire sphere.
C
C On input:
C
C N = Number of nodes. N .GE. 3 and N .GE. 7 if
C IFLAG .NE. 1.
C
C X,Y,Z = Arrays containing Cartesian coordinates of
C the nodes. X(I)**2 + Y(I)**2 + Z(I)**2 = 1
C for I = 1 to N.
C
C F = Array containing data values. F(I) is associ-
C ated with (X(I),Y(I),Z(I)).
C
C LIST,LPTR,LEND = Data structure defining the trian-
C gulation. Refer to STRIPACK
C Subroutine TRMESH.
C
C IFLGS = Tension factor option:
C IFLGS .LE. 0 if a single uniform tension
C factor is to be used.
C IFLGS .GE. 1 if variable tension is desired.
C
C SIGMA = Uniform tension factor (IFLGS .LE. 0), or
C array containing tension factors associated
C with arcs in one-to-one correspondence with
C LIST entries (IFLGS .GE. 1). Refer to Sub-
C programs GETSIG, SIG0, SIG1, and SIG2.
C
C NROW = Number of rows in the dimension statement of
C FF.
C
C NI,NJ = Number of rows and columns in the uniform
C grid. 1 .LE. NI .LE. NROW and 1 .LE. NJ.
C
C PLAT,PLON = Arrays of length NI and NJ, respective-
C ly, containing the latitudes and
C longitudes of the grid lines.
C
C IFLGG = Option indicator:
C IFLGG = 0 if gradient estimates at the ver-
C tices of a triangle are to be
C recomputed for each grid point in
C the triangle and not saved.
C IFLGG = 1 if gradient estimates are input in
C GRAD.
C IFLGG = 2 if gradient estimates are to be
C computed once for each node (by
C GRADL) and saved in GRAD.
C
C The above parameters are not altered by this routine.
C
C GRAD = 3 by N array whose columns contain the X, Y,
C and Z components (in that order) of the grad-
C ients at the nodes if IFLGG = 1, array of
C sufficient size if IFLGG = 2, or dummy para-
C meter if IFLGG = 0.
C
C Gradient estimates may be computed by Subroutines GRADL or
C GRADG if IFLGG = 1.
C
C FF = NROW by NCOL array with NROW .GE. NI and NCOL
C .GE. NJ.
C
C On output:
C
C GRAD = Array containing estimated gradients as de-
C fined above if IFLGG = 2 and IER .GE. 0.
C GRAD is not altered if IFLGG .NE. 2.
C
C FF = Interpolated values at the grid points if IER
C .GE. 0. FF(I,J) = F(PLAT(I),PLON(J)) for I =
C 1,...,NI and J = 1,...,NJ.
C
C IER = Error indicator:
C IER = K if no errors were encountered and K
C grid points required extrapolation for
C K .GE. 0.
C IER = -1 if N, NI, NJ, or IFLGG is outside its
C valid range.
C IER = -2 if the nodes are collinear.
C IER = -3 if extrapolation failed due to the
C uniform grid extending too far beyond
C the triangulation boundary.
C
C STRIPACK modules required by UNIF: GETNP, JRAND, LSTPTR,
C STORE, TRFIND
C
C SSRFPACK modules required by UNIF: APLYR, APLYRT, ARCINT,
C ARCLEN, CONSTR,
C FVAL, GIVENS, GRADL,
C HVAL, INTRC1,
C ROTATE, SETUP,
C SNHCSHI've created a wrapper for UNIF in PR #55.
plon = np.linspace(0, np.pi, 10)
plat = np.linspace(0, np.pi/2, 5)
interp = mesh.interpolate_to_grid(plon, plat, field)gives the same result as
lonq, latq = np.meshgrid(plon, plat)
interp, ier = mesh.interpolate(lonq.ravel(), latq.ravel(), field, order=1)
interp = interp.reshape(plat.size, plon.size)But what I think @kwagstyl wants is to input values of sigma, which is an optional array for storing tension factors. The best thing to do is call a function that stores the optimal tension sigma on the object that would be fed into all smoothing and interpolation routines. However, this would require wrapping SIG0, SIG1, SIG2 and GETSIG.
@kwagstyl - how are you at wrapping Fortran? 😉
@brmather
Thanks so much for doing that.
You are right, the tension sigma would necessary to figure this out.
I'm currently starting at zero for wrapping Fortran, but from your pull request integrating unif it looks like removing the ! comments from ssrfpack.pyf.
Apart from doing the same for the other functions and checkingit runs locally, is there more to it?
Happy to give that a try
Konrad
The same options should be made available to the smoothing routine, right, as we should provide the option to obtain these values point-by-point or for a whole triangulation - the regular grid is somewhat limiting. I would suggest doing the wrapping and then consider whether a more python-flavoured interface would be more helpful.
@kwagstyl it's not too difficult to wrap the functions, just beware of the intent(in) / intent(out) statements which tells python whether a variable is input or output. The other wrapped routines can be used as a guide alongside the original paper by Renka. You can submit a PR on the same UNIF branch and tag it in this issue. We can squash any bugs in the PR.
@lmoresi I'm not sure there is a fortran function for local smoothing, but FVAL handles local interpolation and GRADL local gradients. Ideally, I would like self.sigma=0 at __init__ then if the user wants to apply tension they can call a wrapper for SIG0 etc, which overwrites self.sigma with another array of length 6n-12.
@lmoresi I'm not sure there is a fortran function for local smoothing, but
FVALhandles local interpolation andGRADLlocal gradients. Ideally, I would likeself.sigma=0at__init__then if the user wants to apply tension they can call a wrapper forSIG0etc, which overwritesself.sigmawith another array of length6n-12.
I wasn't suggesting a point-wise computation, only that we should offer sampling at any points not just at points on a regular grid. Do the tensioned splines guarantee node values are honoured (so infinite tension -> bilinear interpolation) ?
I wasn't suggesting a point-wise computation, only that we should offer sampling at any points not just at points on a regular grid.
We already have unstructured interpolation wrapped in the interpolate method, just sigma needs to be input to the Fortran subroutines for tensioned splines.
Do the tensioned splines guarantee node values are honoured (so infinite tension -> bilinear interpolation) ?
Good question - another notebook in the examples documentation should probably explore this.
After beating #56 into #55 the tensioned splines stuff is ready to go. @kwagstyl could you test if it works correctly? The spline tensions can be updated with update_tension_factors and the smoothing and gradient operations should work. sig0, sig1, and sig2 are all accessible under ssrfpack. i.e.
from stripy import _ssrfpack
_ssrfpack.sig0(...)Unfortunately, interpolation is not working with this yet 😞 sigma is hardcoded to equal zero in stripack.f90 so the fortran needs some gentle modification 🔨
@brmather @lmoresi
The code runs and the output looks like something in between linear and cubic, so that seems to be working :)
The Unif function seems to automatically create a grid mesh, and I'm interpolating to a set of points.
So I'm using it like this:
#set up triangulation with original samples lons_sampled (180 samples)
mesh=stripy.sTriangulation(lons_sampled,lats_sampled)
#set up regular dense grid
lats_grid=np.linspace(-np.pi/2,np.pi/2,500)
lons_grid=np.linspace(-np.pi,np.pi,500)
#Do interpolation with splines under tension
ff = mesh.interpolate_to_grid(lons_grid, lats_grid,z_sampled)
#create triangulation for grid
mesh_g = stripy.sTriangulation(lo,la,permute=True)
#local linear interpolation to irregular spherical coordinates of interest
fs=mesh_g.interpolate(lons_sph, lats_sph, ff.T.ravel(),order=1)
Also, to get the orientation right, I had to transpose the grid output (ff.T.ravel()).
Hi @kwagstyl
Thanks for testing the code. I've overhauled the interpolate_cubic method to accept tensioned splines, so you can interpolate unstructured coordinates without creating an intermediary mesh object from the grid coordinates.
Good spot for finding the transposed array 👏 - switching between C and Fortran ordering can be confusing so I've enforced the C ordering in the interpolate_to_grid method.
We will merge all of this in to the master branch and create a new stripy release shortly.