Visualizing data in Cubed-Sphere grid

IN PROGRESS

Cubed-Sphere grid background

The GEOS-5 now has products that are stored natively in cubed-sphere grid. The utilities listed here may help users visualize the data in cubed-sphere grid or map the data to lat-lon grid.

Fortran

To map the cubed-sphere grid data to lat-lon gird data, the program c2l_CFIO_offline.x is built with AGCM at src/GEOSgcs_GridComp/GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSsuperdyn_GridComp/FVdycoreCubed_GridComp. The example configure file c2l_CFIO_offline.rc is also generated. It should be run as mpirun –np 6 ./c2l_CFIO_offline.x. This will produce a lat-lon grid data file whose name is configured in c2l_CFIO_offline.rc. Users then can use their familiar tools to visualize the lat-lon grid data.

Matlab

Users can find a “GridUtils” package in MATLAB ( R2015b and R2016a) at http://wiki.seas.harvard.edu/geos-chem/index.php/GEOS-Chem_HP_Output_Data#MATLAB_.28Cubed-sphere_and_regular_data.29 . The instructions are straight forward. Although it is intended for GCHP, it can be used for all the other outputs that use the same formats as that of GCHP. For the new format of the output, users can use the matlab code here to transform the data to lat-lon grid and visualize it.

close all;
clear all;

% users specify the file name, variables, the the lat-lon resolution, level and time. 
% level is optional 

fname='example_prog.20141118_2130z.nc4';
varname='T';
Nlat=100;
Nlon=600;
level=1;
time=1;

[grid_data,X,Y]=CS2latlon(fname,varname,Nlat,Nlon,time,level);
pcolor(Y',X',squeeze(grid_data)'),shading('flat'),colorbar

% save below as CS2latlon.m

function [grid_data,X,Y]=CS2latlon(fname,varname,Nlat,Nlon,time,level)
   ncid=netcdf.open(fname,'NOWRITE');
   vinfo=ncinfo(fname,'T')
   S=size(vinfo.Size);
   ydim=vinfo.Size(1);
   xdim=vinfo.Size(2);
   nf=vinfo.Size(3);
   netcdf.close(ncid);

   lats=ncread(fname,'lats');
   lons=ncread(fname,'lons');
   if( S(2)==5)
      T=ncread(fname,'T',[1 1 1 level time],[ydim xdim nf level time]);
   else
      % surface, no level
      T=ncread(fname,'T',[1 1 1 time],[ydim xdim nf time]);
   end

   length = xdim*ydim*nf;

   lats=reshape(lats,1,length);
   lons=reshape(lons,1,length);
   T=reshape(T,1,length);
   lats=lats(~isnan(T));
   lons=lons(~isnan(T));
   T=T(~isnan(T));

   latlim=linspace(min(lats),max(lats),Nlat);
   lonlim=linspace(min(lons),max(lons),Nlon);
   [X,Y]=meshgrid(latlim,lonlim);

   grid_data=griddata(lats,lons,T,X,Y,'v4');

Python

For the new format of the output, users can use the python code here to transform the data to lat-lon grid and visualize it.

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
import CSgrid

def main() :

   Time=0
   Level=0
   Nlat=500
   Nlon=1000
   filename = 'example_prog.20141118_2130z.nc4'
   #filename = 'MAX_Discharge.geosgcm_surf.20141118_2130z.nc4'
   var='T'
   #var = 'PHIS'
   # convert data trough interpolation
   grid_T=CSgrid.Convert2LatLon(filename,var,Nlat,Nlon,Time,Level)
   #grid_T=CSgrid.Convert2LatLon(filename,var,Nlat,Nlon,Time)
   # color plot
   plt.figure()
   ax=plt.gca()
   plt.xlabel('lon')
   plt.ylabel('lat')
   im=plt.imshow(grid_T.T,extent=(0,360,-90,90),origin='lower')
   divider = make_axes_locatable(ax)
   cax = divider.append_axes("right", size="5%", pad=0.05)
   plt.colorbar(im, cax=cax)
   plt.title(var)
   plt.show()
if __name__ == '__main__' :
   main()

# save below as CSgrid.py

import sys
from netCDF4 import Dataset
import numpy as np
from scipy.interpolate import griddata

def Convert2LatLon(filename,var,Nlat,Nlon,Time,Level=None):
   fh=Dataset(filename,mode='r')

   Times=fh.dimensions['time'].size
   #Levs=fh.dimensions['lev'].size
   nf=fh.dimensions['nf'].size
   Ydim=fh.dimensions['y'].size
   Xdim=fh.dimensions['x'].size
   if Level is None:
      T=fh.variables[var][Time,:,:,:]
   else :
      T = fh.variables[var][Time,Level, :, :, :]
   lon3d=fh.variables['lons'][:,:,:]
   lat3d=fh.variables['lats'][:,:,:]
   fh.close()

   lat3d=np.reshape(lat3d,(nf*Ydim*Xdim,))
   lon3d=np.reshape(lon3d,(nf*Ydim*Xdim,))
   T=np.reshape(T,(nf*Ydim*Xdim,))
   Tm=T
   # get rid of non values
   if(hasattr(T,'mask')):
      Tm=T[~T.mask]
      lat3d=lat3d[~T.mask]
      lon3d=lon3d[~T.mask]

   lat=np.linspace(min(lat3d),max(lat3d),Nlat)
   lon=np.linspace(min(lon3d),max(lon3d),Nlon)

   grid_x,grid_y=np.meshgrid(lat,lon)
   points=(lat3d,lon3d)
   grid_T=griddata(points,Tm,(grid_x,grid_y),'linear')

   # filled in extrapolated value with nearest

   lat3d =np.reshape(grid_x,(Nlon*Nlat,))
   lon3d =np.reshape(grid_y,(Nlon*Nlat,))
   Tm=np.reshape(grid_T,(Nlon*Nlat,))
   lat3d=lat3d[~np.isnan(Tm)]
   lon3d=lon3d[~np.isnan(Tm)]
   Tm=Tm[~np.isnan(Tm)]
   points=(lat3d,lon3d)

   grid_T=griddata(points,Tm,(grid_x,grid_y),'nearest')

   return grid_T

IDL

The script read_and_interpolate_cube.pro can read data in the cubed-sphere grid in NetCDF files, interpolate the data , map the data in lat-lon format and plot. The script cube_to_latlon.pro uses read_and_interpolate_cube.pro to convert the variables in cubed-sphere grid to lat-lon grid and write the variables to NetCDF files. It should be noted that the tile files, input and output files are hard-wired in the code. The users will need to provide those files and change the code accordingly.

GrADS

To use GrADS to display cube-sphere grid data in a NetCDF file, the path should be in the PATH environment: /discover/nobackup/projects/gmao/share/dasilva/opengrads/Contents/grads

Then the users can follow the steps. First, we assume a version of GEOS-5 is available.

1. run $ESMADIR/src/GMAO_Shared/GEOS_Util/plots/configure. This will produce a .quickplotrc in that directory
2. source $ESMADIR/src/GMAO_Shared/GEOS_Util/plots/.quickplotrc
3. run grads
4. sdfopen file.nc4
5. run command dc to plot the cube-sphere grid NetCDF file. Users can get help by just run dc without any argument.