Spherical Interpolation of data

I have a set of data points which are defined at various latitudes and longitudes. I want to interpolate the data points along the surface of a sphere, so I basically have a data point at every latitude and longitude coordinate (to .1 place)

At the moment, I have this code - which works:

        PTABLE = readtable('coordinates.dat');
        POINTS = table2array(PTABLE);
        VMRTABLE = readtable('values.dat');
        VMR = table2array(VMRTABLE(:,2));

        LON = POINTS(:,1);
        LAT = POINTS(:,2);
        VAL = VMR(:,1);

        % Convert coordinates to radians
        LATr = LAT*pi/180;
        LONr = LON*pi/180;

        % Create a longitude/latitude grid with 0.1 degree spacing
        nrows = 1800;
        ncols = 3600;
        X = linspace(0,360,ncols);
        Y = linspace(-90,90,nrows);
        [lon, lat] = meshgrid(X,Y);

        % Convert to radians and make longitude measured "westwards" (-ve sign)
        lon = -lon*pi/180;
        lat = lat*pi/180;

        % Convert the coordinate data points to cartesian
        [x, y, z] = sph2cart(LONr,LATr,1);

        % Convert the gridded points to cartesian
        [xg, yg, zg] = sph2cart(lon,lat,0.998);

        % Perform the interpolation, with no extrapolation
        F=scatteredInterpolant(x,y,z,VAL,'linear','none');
        vq = F(xg,yg,zg);

        save('data_interp.dat', 'vq', '-ascii')

However, in order to perform the interpolation, I am essentially creating a 3D domain defined by my coordinate data points, which is then interpolated directly between points (not along the surface of a sphere). This is why I need to make the radius of the gridded circle slightly smaller, otherwise all points on it are "outside" the domain of the data points, and so purely extrapolation is happening, which I do not want.

My question is this:

Is there a method which can directly interpolate my values across the surface of the sphere? I have looked into slerp, and using quaternions, but there doesn't seem to be a way to actually interpolate values, rather than rotations via that method.

I hope this is clear, and any help is much appreciated!

[Michael Nolan]

I've been trying this same thing. If you try to convert to cartesian and use scatteredInterpolant, it works OK if you pick a random point: F(.1,.1,.1) in the notation of the scatteredInterpolant documentation, but if you try to actually choose a point on the sphere: F(cos(lon) * cos (lat), sin(lon)*cos(lat), sin(lat)), matlab hangs with 100% CPU. I'm guessing that this has to do with the degenracy of the triangles, but haven't been able to solve it. Adding some noise to the radius of the interpolated values makes it not fail, but the results aren't good. I tried adding a point in the center of the sphere to try to break the degeneracy, but that didn't help: Matlab still hangs when asked for even one point on the sphere.

For my current purpose, I'm giving up and just doing in in cylindrical coordinates, but that's clearly wrong in general.

Answer 1

[John D'Errico]

So, you have some function, known at points only on the surface of a sphere, where the locations are given by Latitude and Longitude. Now you want to do interpolation. This is no problem if your points don't cross the 0/360 degree boundary, or lie around the poles. Were that not a problem, you would just use standard interpolation methods in the 2-d space of latitude & longitude.

But what do you do when points cross those boundaries? Um, yes, there are methods. I don't have the time to write in any depth now though. Fermat once said something very similar, I recall. ;-)

The wrong approach it to convert to cartesian coordinates, and then try to do interpolation, because that assumes the points lie in 3-dimensions, when in reality, the points lie on the surface of a sphere.

It would seem more logical to compute inter-point great-circle distances between the points. (Easy to compute. See the link.) Now you can use any of various distance based interpolations. This works because great-circle distance crosses the 0/360 boundaries in any direction. Problems at the poles go away. Sorry if I can't go any further right now. But that is the gist of it. I'll try to look in later.

[Jason Sharkey]

There are typically several hundred data points, and so the spacing between them has been relatively small. So the as a general rule, the distance between them in cartesian space and the distance along a great arc between them have been relatively similar. So, I have been satisfied with my solution thus far, but I was wondering if there was a method, or a function in MATLAB to do it more accurately - it would be a much worse problem if my points were not close together.

I'm surprised that I can't find anyone who has wanted to do something similar online.

[Pratik Dash]

Has there been any update to this ? I am also in need of a function that interpolates on the surface of an ellipsoid. What I have are polar and azimuthal angles and my function is time-dependent ? So I will have to perform 3D interpolation.