GEOS-5

Visualizing data in Cubed-Sphere grid: Difference between revisions

← Older edit
No edit summary
 
(14 intermediate revisions by the same user not shown)
Line 1: Line 1:
{{rightTOC}}
{{rightTOC}}
<span style="color: red; font-weight: bold; font-size: 24pt">IN PROGRESS</span>
 


== Cubed-Sphere grid  background ==
== Cubed-Sphere grid  background ==
Line 14: Line 14:


#Download source codes and save them as driver.m, CSnative.m, extendFace1.m,extendFace3.m, extendFace4.m, extendFace6.m, and the example data file respectively to the same folder. The extendFace(1-6).m are used to fill the seam between cubed-sphere faces. Face2 and Face5 are not necessary to be  extended because they are redundant.
#Download source codes and save them as driver.m, CSnative.m, extendFace1.m,extendFace3.m, extendFace4.m, extendFace6.m, and the example data file respectively to the same folder. The extendFace(1-6).m are used to fill the seam between cubed-sphere faces. Face2 and Face5 are not necessary to be  extended because they are redundant.
#Users should specify the time, level , variable’s name and file name in driver.m
#Users should specify the time, level , variable’s name and file name in driver.m. The test data file can be downloaded [[Media:TEST7.geosgcm_prog.20000415_0000z.nc4]]
#Run with Matlab: % driver
#Run with Matlab: % driver


Line 21: Line 21:
==== [[Recipe: Matlab program  fills the seam| extendFace1.m]] ====
==== [[Recipe: Matlab program  fills the seam| extendFace1.m]] ====
==== [[Recipe: Matlab program  fills the seam3| extendFace3.m]] ====
==== [[Recipe: Matlab program  fills the seam3| extendFace3.m]] ====
 
==== [[Recipe: Matlab program  fills the seam4| extendFace4.m]] ====
==== [[Recipe: Matlab program  fills the seam| extendFace4.m]] ====
==== [[Recipe: Matlab program  fills the seam6| extendFace6.m]] ====
==== [[Recipe: Matlab program  fills the seam| extendFace6.m]] ====
 
 
<syntaxhighlight lang="matlab" line>
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');
</syntaxhighlight>


== Python ==  
== 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.
For the new format of the output, users can use the python codes here to transform the data to lat-lon grid and visualize the products:
 
#Install python 3 (or above) and the modules numpy, netcdf4, matplotlib and mpl_toolkits  on your computer.  
<syntaxhighlight lang="python" line>
#Download source codes CSnative.py, example1.py, example2.py and the example data file  [[Media:TEST7.geosgcm_prog.20000415_0000z.nc4]] the same folder.
import numpy as np
#To show example 1, $python example1.py
import matplotlib.pyplot as plt
#To show example 2, $python example2.py
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
==== [[Recipe: python program  reads cubed-sphere data| CSnative.py]]====
</syntaxhighlight>
==== [[Recipe: python program  example1 | example1.py]] ====
==== [[Recipe: python program  example2| example2.py]] ====


== IDL ==
== IDL ==
To use IDL codes and the pbs scripts to visualize a large number of GEOS-5 cubed-sphere products:
#Download the codes and save them as cube_to_latlon.pro, winds_cube_to_latlon.pro, convert_latlon.j (pbs script) , idlstart ( environment ) and README.
#The README file explains all the scripts and the steps visualizing the data. Although all the examples and data are on discover, it is not difficult for users to make changes to run the programs on the other systems.


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.
==== [[Recipe: README| README]] ====
==== [[Recipe: environment | idlstart]]====
==== [[Recipe: pbs script | Ops12KM_R_to_latlon.j]]====


==== [[Recipe: IDL program to read cubed-sphere data| read_and_interpolate_cube.pro]] ====
==== [[Recipe: cube_lat_laon |cube_to_latlon.pro]]====
==== [[Recipe: IDL program to convert cubed-sphere data to lat-lon data | cube_to_latlon.pro]] ====
==== [[Recipe: winds | winds_cube_to_latlon.pro]]====


== GrADS ==
== GrADS ==
Line 188: Line 58:
# sdfopen file.nc4
# sdfopen file.nc4
# run command <tt>dc</tt> to plot the cube-sphere grid NetCDF file. Users can get help by just run dc without any argument.
# run command <tt>dc</tt> to plot the cube-sphere grid NetCDF file. Users can get help by just run dc without any argument.
==Panoply==
Now the Panoply after version 4.7.0 can view the native cubed-sphere products. This software can be downloaded from https://www.giss.nasa.gov/tools/panoply/
==Converting (interpolating) cubed-sphere data to Lat-Lon data==
For some purposes it will be impractical to adapt existing analysis tools to directly work with MERRA cubed-sphere products.  In such cases, users will need to build an executable that can interpolate cubed-sphere data to a set of predefined lat-lon resolutions.  (Note that special care must be taken for vector data, e.g., (u,v).)
The building requirements and instruction can be found here  https://geos5.org/wiki/index.php?title=Building_Baselibs . The executable program cube2laton can be used to convert (or interpolate) the data
Converting data:
%  cube2latlon <in-file> <out-file> <target-resolution>
E.g.,
% cube2latlon  merra.nc  latlon.nc  4x5
Retrieved from "https://geos5.org/wiki/index.php?title=Visualizing_data_in_Cubed-Sphere_grid"