Spherical Grid Conversion

• ACTM Performer: AIBEDO Team (PARC, UVic and UW);

• Author: Haruki Hirasawa (hhirasawa@uvic.ca)

We use pygsp package to interpolate 2D-gridded data into a spherical grid format. This is done in two steps:

Step 1: Generate a text file (“skeleton”) with the desired Icosahedral level using our python script link here.

This code block generates a text file that will be used to generate the spherical sample for level 6. To generate a text file for another grid level, please change the sphlevel in the code.

Step 2: Once the text file is generated in step 1, we use the cdo (Climate Data Operator) command line tool to generate the interpolated netCDF file. Please see here for instructions to download cdo.

The following script is given in command line to generate the interpolated file for model training:

cdo remapbil,icosphere_6.txt in.nc out.nc

Here, in.nc is the 2D-gridded file from ESM model ensembles or Reanalysis datasets, and out.nc is the name of the interpolated file that will be used for model training.

import pygsp as pg
import numpy as np
### Generates a text file with centre points of a icosahedral sphere grid
### Formatted to pass as an argument for 'cdo remap,GRIDFILE in.nc out.nc'
def xyz2lonlat(x:float,y:float,zLfloat, radius=6371.0e6):
    """From cartesian geocentric coordinates to 2D geographic coordinates."""
    """ lon: [-180, 180], lat:[-90, 90] """
    latitude = np.arcsin(z / radius)/np.pi*180
    longitude = np.arctan2(y, x)/np.pi*180
    return longitude, latitude
def gen_icosphere(level:int,radius:float = 6371.0e6):
    '''
    Generates points on a icosphere grid (from Soo Kim)
    level (int): iteration level
    radius (float): radius of the earth
    '''
    #(1) generate graph of SphereIcosahedral
    graph = pg.graphs.SphereIcosahedral(2**level)
    #(2) extract lon, lat coordinate
    vertices = graph.coords # xyz coordinates
    #print("number of vertices for level "+str(level)+" (2**"+str(level)+")is "+str(graph.n_vertices))
    lon = []
    lat = []
    radius = 1
    # convert cartesian points to spherical lon/lat
    for tmp_xyz in vertices:
        tmp_lon, tmp_lat = xyz2lonlat(tmp_xyz[0],tmp_xyz[1],tmp_xyz[2], radius=radius)
        lon.append(tmp_lon)
        lat.append(tmp_lat)
    return lon, lat
sphlevel = 6
lon,lat = gen_icosphere(sphlevel,radius = 6371.0e6)
f = open("isosphere_{0}.txt".format(sphlevel), "a")
# write grid to file
f.write("gridtype = unstructured\n")
f.write("gridsize = {0}\n".format(len(lon)))
f.write("# Longitudes\n")
f.write("xvals = " + ' '.join(map(str, lon))+"\n")
f.write("# Latitudes\n")
f.write("yvals = " + ' '.join(map(str, lat))+"\n")