CODEWITHSUNDEEP
pythonscipyinterpolation
ryanjdillon
ryanjdillon

Reputation: 18978

How to back-calculate latitude and longitude from scipy.interpolate.RectSphereBivariateSpline?

Using THIS example from the RectSphereBivariateSpline function in the scipy.interpolate sub-package, I would like to back calculate arrays with latitude and longitude in degrees, and an array with the interpolated data value for each coordinate on the grid.

The interpolated data object RectSphereBivariateSpline creates appears to be the u and v components of the data values, with one value for each degree change of the grid (Latitude = 180 and Longitude = 360 for this example).

Is that right?

How might I back-calculate latitude, longitude, and their respective data values for plotting?

import numpy as np
from scipy.interpolate import RectSphereBivariateSpline

def geo_interp(lats,lons,data,grid_size_deg):
        '''We want to interpolate it to a global one-degree grid'''
        deg2rad = np.pi/180.
        new_lats = np.linspace(grid_size_deg, 180, 180) * deg2rad
        new_lons = np.linspace(grid_size_deg, 360, 360) * deg2rad
        new_lats, new_lons = np.meshgrid(new_lats, new_lons)

        '''We need to set up the interpolator object'''
        lut = RectSphereBivariateSpline(lats*deg2rad, lons*deg2rad, data)

        '''Finally we interpolate the data. The RectSphereBivariateSpline 
        object only takes 1-D arrays as input, therefore we need to do some reshaping.'''
        data_interp = lut.ev(new_lats.ravel(),
                        new_lons.ravel()).reshape((360, 180)).T
        return data_interp

if __name__ == '__main__':
        import matplotlib.pyplot as plt

        '''Suppose we have global data on a coarse grid'''
        lats = np.linspace(10, 170, 9) # in degrees
        lons = np.linspace(0, 350, 18) # in degrees
        data = np.dot(np.atleast_2d(90. - np.linspace(-80., 80., 18)).T,
                        np.atleast_2d(180. - np.abs(np.linspace(0., 350., 9)))).T

        '''Interpolate data to 1 degree grid'''
        data_interp = geo_interp(lats,lons,data,1)

        '''Looking at the original and the interpolated data, 
        one can see that the interpolant reproduces the original data very well'''
        fig = plt.figure()
        ax1 = fig.add_subplot(211)
        ax1.imshow(data, interpolation='nearest')
        ax2 = fig.add_subplot(212)
        ax2.imshow(data_interp, interpolation='nearest')
        plt.show()

I thought I might use vector addition (i.e. Pythagorean theorem), but this doesn't work as I only have one value for each degree change, not point

pow((pow(data_interp[0,:],2.0)+pow(data_interp[:,0],2.0)),1/2.0)

Upvotes: 0

Views: 788

Answers (1)

niels
niels

Reputation: 145

It seems that the latitude and longitude vectors that can be used for plotting are generated in our "geo_interp" function ("net_lats" and "new_lons"). If you require these in your main program, you should either declare these vectors outside of your function, or you should have the function return these generated vectors to be used in the main program.

Hope this helps.

Upvotes: 1

Related Questions