Search code examples
pythonprojectioncoordinate-systems

Python: convert 1-D array (with equal area coordinate system) to a 2-D array (with Geographic Coordinate Reference System)


I have a 1-D array of data (e.g. Precipitation [precip]). Also, I have 1D latitude (min -90 deg., max +90 deg.) and 1D longitude (min 0, max 360 deg.) arrays representing the coordinates of this data. The coordinate system is "equal area". It is a global dataset.

My question is how can I convert this 1-D array to a 2-D array with Geographic Coordinate Reference System (i.e. equally spaced grid, both parallels, and meridians) with a spatial resolution of 1 by 1 degree, so that I would have a 180*360 array (preferably, using pyproj / xarray)?

Thanks!

The following is the information of the dataset:

xarray.Dataset

Dimensions: (eqcell: 41252)

Dimensions without coordinates: eqcell

Data variables:

lat                (eqcell) float32 dask.array chunksize=(41252,), meta=np.ndarray

lon                (eqcell) float32 dask.array chunksize=(41252,), meta=np.ndarray

precip              (eqcell) float32 dask.array chunksize=(41252,), meta=np.ndarray

Solution

  • It looks like you want scipy.interpolate.griddata. Here's the example from the documentation:


    Suppose we want to interpolate the 2-D function

    >>> def func(x, y):
    ...     return x*(1-x)*np.cos(4*np.pi*x) * np.sin(4*np.pi*y**2)**2
    

    on a grid in [0, 1]x[0, 1]

    >>> grid_x, grid_y = np.mgrid[0:1:100j, 0:1:200j]
    

    but we only know its values at 1000 data points:

    >>> points = np.random.rand(1000, 2)
    >>> values = func(points[:,0], points[:,1])
    

    This can be done with griddata – below we try out all of the interpolation methods:

    >>> from scipy.interpolate import griddata
    >>> grid_z0 = griddata(points, values, (grid_x, grid_y), method='nearest')
    >>> grid_z1 = griddata(points, values, (grid_x, grid_y), method='linear')
    >>> grid_z2 = griddata(points, values, (grid_x, grid_y), method='cubic')
    

    One can see that the exact result is reproduced by all of the methods to some degree, but for this smooth function the piecewise cubic interpolant gives the best results:

    >>> import matplotlib.pyplot as plt
    >>> plt.subplot(221)
    >>> plt.imshow(func(grid_x, grid_y).T, extent=(0,1,0,1), origin='lower')
    >>> plt.plot(points[:,0], points[:,1], 'k.', ms=1)
    >>> plt.title('Original')
    >>> plt.subplot(222)
    >>> plt.imshow(grid_z0.T, extent=(0,1,0,1), origin='lower')
    >>> plt.title('Nearest')
    >>> plt.subplot(223)
    >>> plt.imshow(grid_z1.T, extent=(0,1,0,1), origin='lower')
    >>> plt.title('Linear')
    >>> plt.subplot(224)
    >>> plt.imshow(grid_z2.T, extent=(0,1,0,1), origin='lower')
    >>> plt.title('Cubic')
    >>> plt.gcf().set_size_inches(6, 6)
    >>> plt.show()
    

    plot
    (source: scipy.org)