Visualizing data in Cubed-Sphere grid: Difference between revisions

← Older edit Newer edit →

@@ Line 61: / Line 61: @@
 </syntaxhighlight>
+== Python ==
+For the new format of the output, user can use the python code here to transform the data to lat-lon grid and visualize it.
+<syntaxhighlight lang="python" line>
+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 = 'MAX_Discharge.geosgcm_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
+</syntaxhighlight>
 == IDL ==