Visualizing data in Cubed-Sphere grid: Difference between revisions

Line 11: Line 11:


== Matlab ==
== 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.
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, user can use the code here to transform the dat to lat-lon grid and visualize it.
 
close all;
clear all;
fname='MAX_Discharge.geosgcm_prog.20141118_2130z.nc4';
 
% user specify the file name, variables, the the lat-lon resolution, level and time.
% level is optional
 
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
 
function [grid_data,X,Y]=CS2ltlon(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 1 1],[ydim xdim nf level time]);
  else
      T=ncread(fname,'T',[1 1 1 1],[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');


== IDL ==
== IDL ==