GEOS-5

Visualizing data in Cubed-Sphere grid: Difference between revisions

← Older edit
No edit summary
 
(33 intermediate revisions by 2 users not shown)
Line 1: Line 1:
#==== in progress ====
{{rightTOC}}
 
 
== Cubed-Sphere grid  background ==
== 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.
Several GEOS-5 products are now produced on the native cubed-sphere computational grid.   Although this format may be less familiar, this choice has several advantages for some user communities.  In particular, the Chemical Transport Modeling (CTM) community can use native mass-flux products to significantly improve conservation properties when running with GEOS-5. To minimize transition difficulties for GEOS-5 data consumers, this page provides instructions for visualizing cubed-sphere data products using a variety of common packages:  python, Matlab, IDL, and GRaDS, and Panoply. This page also provides instructions for downloading and building an executable that can regrid cubed-sphere products onto a traditional lat-lon grid.


== Fortran ==
== 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  
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 ==
== 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 uses the same formats as that of GCHP.
Users can find a “GridUtils” package in MATLAB ( R2015b and R2016a or above) 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 codes here to transform the data to lat-lon grid and visualize the products.


== IDL ==
#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. The test data file can be downloaded [[Media:TEST7.geosgcm_prog.20000415_0000z.nc4]]
#Run with Matlab: % driver


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 are hard-wired in the code. The users will need to provide those files and change the code accordingly.
==== [[Recipe: Matlab program  driver for visualization| driver.m]] ====
==== [[Recipe: Matlab program  converts cubed-sphere to latlon| CSnative.m]] ====
==== [[Recipe: Matlab program  fills the seam| extendFace1.m]] ====
==== [[Recipe: Matlab program fills the seam3| extendFace3.m]] ====
==== [[Recipe: Matlab program  fills the seam4| extendFace4.m]] ====
==== [[Recipe: Matlab program  fills the seam6| extendFace6.m]] ====


<syntaxhighlight lang="idl" line>
== Python ==  
; read_and_interpolate_cube.pro
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:
pro read_and_interpolate_cube, myfile, myVar, data, gLons, gLats, LonBeg, LonEnd, LatBeg, LatEnd, undef, TRIANGLES=TR, SCALE_FACTOR=SCALE_FACTOR, PIXEL=PIXEL, FLIP=FLIP, regridCubeToLatlon=regridCubeToLatlon, TIME=TIME, LEVEL=LEVEL, DATA_LEVELS=DATA_LEVELS, RANGE=RANGE, regridResolution=regridResolution, LOG=LOG, NCOLORS=NCOLORS, PROJECT=PROJECT, LISTVARS=LISTVARS, CMORPH_DATA=CMORPH_DATA, TRMM_DATA=TRMM_DATA, FILEID=INFILEID, STRETCH=STRETCH, CONSERV=CONSERV
#Install python 3 (or above) and the modules numpy, netcdf4, matplotlib and mpl_toolkits  on your computer.
#Download source codes CSnative.py, example1.py, example2.py and the example data file  [[Media:TEST7.geosgcm_prog.20000415_0000z.nc4]] the same folder.
#To show example 1, $python example1.py
#To show example 2, $python example2.py


  print, myfile
==== [[Recipe: python program reads cubed-sphere data| CSnative.py]]====
  help, myVar
==== [[Recipe: python program example1 | example1.py]] ====
;help, LonBeg, LonEnd, LatBeg, LatEnd
==== [[Recipe: python program  example2| example2.py]] ====
if (~keyword_set(TIME)) then TIME=0


exists = FILE_TEST(myfile)
== IDL ==
if (exists) then begin
To use IDL codes and the pbs scripts to visualize a large number of GEOS-5 cubed-sphere products:
  if ( keyword_set(INFILEID)) then fileID = INFILEID
#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.
  if (~keyword_set(INFILEID)) then fileID = ncdf_open(myfile, /NOWRITE)
#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.


  if (keyword_set(LISTVARS)) then begin
==== [[Recipe: README| README]] ====
    ;find the number of variables
==== [[Recipe: environment | idlstart]]====
    fileinq_struct = ncdf_inquire(fileID)
==== [[Recipe: pbs script | Ops12KM_R_to_latlon.j]]====
    nvars = fileinq_struct.nvars
    ; print output so far
    print, myfile
    print, nvars, ' variables in file'
    print, ' '
    print, ' '
    ;loop through variables
    for varndx = 0,nvars-1 do begin
    ; get the name, datatype, dims, number of attributes
      varstruct = ncdf_varinq(fileID,varndx)
    ; print the variable index, name, datatype, dims
      print, '--------------------------------------------------------'
      print,varndx,'    ',varstruct.name, '    ',varstruct.datatype
      print,'dims = ',varstruct.dim
    endfor ; variable loop
    stop
  endif


  finq = NCDF_INQUIRE( fileID )        ; finq /str = ndims, nvars, ngatts, recdim
==== [[Recipe: cube_lat_laon |cube_to_latlon.pro]]====
  dim_unlimited = finq.recdim          ; = -1 if undefined, otherwise index into dim array
==== [[Recipe: winds | winds_cube_to_latlon.pro]]====
  if ( finq.ndims ge 2 ) then begin
  dimstr = ' '
  dimsize = 0L
  NCDF_DIMINQ, fileID, 0, dimstr, dimsize
  ;help, dimstr
  NX=dimsize
  NCDF_DIMINQ, fileID, 1, dimstr, dimsize
  ;help, dimstr
  NY=dimsize
  endif
  NZ=0L
  NT=1L
  if ( finq.ndims ge 3 ) then begin
  NCDF_DIMINQ, fileID, 2, dimstr, dimsize
  ;help, dimstr
  if (dimstr eq 'time') then begin
    NT=dimsize
  endif else begin
    NZ=dimsize
  endelse
  endif else begin
    if ( finq.ndims ge 4 ) then begin
    NCDF_DIMINQ, fileID, 3, dimstr, dimsize
    ;help, dimstr
    NT=dimsize
    endif
  endelse
;help, NX, NY, NZ, NT


  if (keyword_set(LEVEL)) then begin
== GrADS ==
  if (NZ ge 1) then begin
    varID=ncdf_varid(fileID,'lev')
    ncdf_varget,fileID,varID,levs
  endif
  endif


  if (myVar eq 'WSPD50M') then begin
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


    varID=ncdf_varid(fileID,'U50M')
Then the users can follow the steps. First, we assume a version of GEOS-5 is available.
    ncdf_varget,fileID,varID,u
# run $ESMADIR/src/GMAO_Shared/GEOS_Util/plots/configure. This will produce a .quickplotrc in that directory
    varID=ncdf_varid(fileID,'V50M')
# source $ESMADIR/src/GMAO_Shared/GEOS_Util/plots/.quickplotrc
    ncdf_varget,fileID,varID,v
# run grads
    cdataM = SQRT(u*U + v*v)
# 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.


  endif else begin
==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/


  myVars=STRSPLIT(myVar,'*', /extract)
==Converting (interpolating) cubed-sphere data to Lat-Lon data==
  s=SIZE(myVars)
  if (s(1) eq 1) then begin
    varID=ncdf_varid(fileID,myVar)
    if (keyword_set(LEVEL)) then begin
      print, 'Reading LEVEL: ', LEVEL
      ifind = WHERE(levs eq LEVEL)
      if (ifind(0) ne -1) then begin
        ncdf_varget,fileID,varID,cdataM,COUNT=[NX,NY,1,1], OFFSET=[0,0,ifind(0),TIME]
      endif else begin
        print, 'ERROR: bad level: ', LEVEL
        print, levs
        stop
      endelse
    endif else begin
    ;help, TIME
      ncdf_varget,fileID,varID,cdataM,COUNT=[NX,NY,1], OFFSET=[0,0,TIME]
    endelse
  endif else begin
    varID=ncdf_varid(fileID,myVars(0))
    if (keyword_set(LEVEL)) then begin
      print, 'Reading LEVEL: ', LEVEL
      ifind = WHERE(levs eq LEVEL)
      if (ifind(0) ne -1) then begin
        ncdf_varget,fileID,varID,cdataA,COUNT=[NX,NY,1,1], OFFSET=[0,0,ifind(0),0]
      endif else begin
        print, 'ERROR: bad level: ', LEVEL
        print, levs
        stop
      endelse
    endif else begin
      ncdf_varget,fileID,varID,cdataA
    endelse
    isUNDEF = WHERE(cdataA eq undef)
    if (N_ELEMENTS(isUNDEF) eq N_ELEMENTS(cdataA)) then stop
    if (isUNDEF(0) ne -1) then cdataA(isUNDEF) = 0.0 ;!VALUES.F_NAN
    varID=ncdf_varid(fileID,myVars(1))
    if (keyword_set(LEVEL)) then begin
      print, 'Reading LEVEL: ', LEVEL
      ifind = WHERE(levs eq LEVEL)
      if (ifind(0) ne -1) then begin
        ncdf_varget,fileID,varID,cdataB,COUNT=[NX,NY,1,1], OFFSET=[0,0,ifind(0),0]
      endif else begin
        print, 'ERROR: bad level: ', LEVEL
        print, levs
        stop
      endelse
    endif else begin
      ncdf_varget,fileID,varID,cdataB
    endelse
    isUNDEF = WHERE(cdataB eq undef)
    if (N_ELEMENTS(isUNDEF) eq N_ELEMENTS(cdataB)) then stop
    if (isUNDEF(0) ne -1) then cdataB(isUNDEF) = 0.0 ;!VALUES.F_NAN
    cdataM=cdataA*cdataB
    cdataA=0
    cdataB=0
  endelse


   endelse
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).)


  if (keyword_set(TRMM_DATA)) then begin
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
    nLons = 1440
    nLats = 720
    data = FLTARR(nLons,nLats)
    data(*,*) = !VALUES.F_NAN
    data(*,160:559) = cdataM
    cdataM=data
  endif


  if (keyword_set(CMORPH_DATA)) then begin
    nLons = 1440
    nLats = 720
    data = FLTARR(nLons,nLats)
    data(*,*) = !VALUES.F_NAN
    data(*,120:599) = REVERSE(cdataM,2)
    cdataM=data
    isUNDEF = WHERE(cdataM eq MIN(cdataM,/NAN))
    if (isUNDEF(0) ne -1) then cdataM(isUNDEF) = !VALUES.F_NAN
    flip_hemispheres, cdataM
  endif


  isUNDEF = WHERE(cdataM eq undef)
Converting data:
  if (N_ELEMENTS(isUNDEF) eq N_ELEMENTS(cdataM)) then stop
  if (isUNDEF(0) ne -1) then cdataM(isUNDEF) = !VALUES.F_NAN


  if (keyword_set(SCALE_FACTOR)) then begin
%  cube2latlon <in-file> <out-file> <target-resolution>
      print, 'SCALING by: ', SCALE_FACTOR
      cdataM = cdataM*SCALE_FACTOR
  endif


; Check if data is on cubed-sphere grid
E.g.,
  isCubedSphere=0
  if (NY eq 6*NX) then isCubedSphere=1
  if (ARRAY_EQUAL(cdataM,0.0)) then stop
;isNAN = WHERE(FINITE(cdataM))
;if (N_ELEMENTS(isNAN) ne N_ELEMENTS(cdataM)) then stop
;isUNDEF = WHERE(cdataM eq undef)
;if (N_ELEMENTS(isUNDEF) eq N_ELEMENTS(cdataM)) then stop
  if (~keyword_set(INFILEID)) then ncdf_close,fileID
 
  if (isCubedSphere AND keyword_set(regridCubeToLatlon)) then begin
  data=FLTARR(gLons,gLats)
  get_domain, data, LonBeg, LonEnd, LatBeg, LatEnd
  m=size(data)
  cLons=m(1)
  cLats=m(2)
  endif else begin
  cLons=gLons
  cLats=gLats
  endelse
 
;help, isCubedSphere, regridCubeToLatlon
  if (isCubedSphere AND keyword_set(regridCubeToLatlon)) then begin
    print, myVar
    print, min(cdataM,/NAN), max(cdataM,/NAN)
    if (keyword_set(regridResolution)) then begin
        print, 'Reducing Cubed-Sphere Resolution to: ', regridResolution
        print, NX, NY
        t1 = FLTARR(NX,NX)
        t2 = FLTARR(NX,NX)
        t3 = FLTARR(NX,NX)
        t4 = FLTARR(NX,NX)
        t5 = FLTARR(NX,NX)
        t6 = FLTARR(NX,NX)
        for n=0,5 do begin
          i1=0
          i2=NX-1
          j1=NX*n
          j2=j1+NX-1
          if (n eq 0) then t1(*,*) = cdataM(i1:i2,j1:j2)
          if (n eq 1) then t2(*,*) = cdataM(i1:i2,j1:j2)
          if (n eq 2) then t3(*,*) = cdataM(i1:i2,j1:j2)
          if (n eq 3) then t4(*,*) = cdataM(i1:i2,j1:j2)
          if (n eq 4) then t5(*,*) = cdataM(i1:i2,j1:j2)
          if (n eq 5) then t6(*,*) = cdataM(i1:i2,j1:j2)
        endfor
        NX2 = regridResolution
        XLocs = NX*FINDGEN(NX2)/FLOAT(NX2-1)
        t1 = INTERPOLATE(t1, XLocs, XLocs, /GRID);, CUBIC=-0.5)
        t2 = INTERPOLATE(t2, XLocs, XLocs, /GRID);, CUBIC=-0.5)
        t3 = INTERPOLATE(t3, XLocs, XLocs, /GRID);, CUBIC=-0.5)
        t4 = INTERPOLATE(t4, XLocs, XLocs, /GRID);, CUBIC=-0.5)
        t5 = INTERPOLATE(t5, XLocs, XLocs, /GRID);, CUBIC=-0.5)
        t6 = INTERPOLATE(t6, XLocs, XLocs, /GRID);, CUBIC=-0.5)
      ; t1 = BILINEAR(t1, XLocs, XLocs)
      ; t2 = BILINEAR(t2, XLocs, XLocs)
      ; t3 = BILINEAR(t3, XLocs, XLocs)
      ; t4 = BILINEAR(t4, XLocs, XLocs)
      ; t5 = BILINEAR(t5, XLocs, XLocs)
      ; t6 = BILINEAR(t6, XLocs, XLocs)
        cdataM=[[t1],[t2],[t3],[t4],[t5],[t6]]
      ;help, cdataM
        NX = NX2
        NY = NX2*6
        print, NX, NY
    endif
if (keyword_set(CONSERV)) then begin
    print, 'Conservative Regridding Cube to Latlon'
    NLONS=NX*4
    NLATS=NX*2 + 1
    sNX=STRING(NX,FORMAT='(i4.4)')
    sNLONS=STRING(NLONS,FORMAT='(i4.4)')
    sNLATS=STRING(NLATS,FORMAT='(i4.4)')
    TILE_FILE='/discover/nobackup/ltakacs/bcs/Ganymed-4_0/Ganymed-4_0_Ostia/Shared/DC'+sNLONS+'xPC'+sNLATS+'_CF'+sNX+'x6C.bin'
    print, TILE_FILE
    spawn, 'ls -l '+TILE_FILE
    openr, tile_lun, TILE_FILE, /GET_LUN, /F77_UNFORMATTED
    NT=0L
    NGRIDS=0
    GRIDNAME=''
    IM=0L
    JM=0L
    readu, tile_lun, NT
    readu, tile_lun, NGRIDS
    help, NT, NGRIDS
    for n=1,NGRIDS do begin
      readu, tile_lun, GRIDNAME
      readu, tile_lun, IM
      readu, tile_lun, JM
      help, IM, JM
    endfor
    IIo=LONARR(NT)
    JJo=LONARR(NT)
      Wo=FLTARR(NT)
    IIi=LONARR(NT)
    JJi=LONARR(NT)
      Wi=FLTARR(NT)
    dummy=''
    readu, tile_lun, dummy ;TYPE
    readu, tile_lun, dummy ;X
    readu, tile_lun, dummy ;Y
  ; Output Grid
    dLON  = 360.0/FLOAT(NLONS)
    dLAT  = 180.0/FLOAT(NLATS-1)
    lonsO  = (-180.0 + 0.5*dLON + FINDGEN(NLONS)*dLON)
    latsO  = ( -90.0 + FINDGEN(NLATS)*dLAT)
    readu, tile_lun,  Wo
                  IIo=Wo
    readu, tile_lun,  Wo
                  JJo=Wo
    readu, tile_lun,  Wo
  ; Input Grid
    readu, tile_lun,  Wi
                  IIi=Wi
    readu, tile_lun,  Wi
                  JJi=Wi
    readu, tile_lun,  Wi
  ; DONE
    close, tile_lun
    free_lun, tile_lun
    help, IIo, JJo, Wo
    print, MIN(IIo), MAX(IIo)
    print, MIN(JJo), MAX(JJo)
    print, MIN( Wo), MAX( Wo)
    help, IIi, JJi, Wi
    print, MIN(IIi), MAX(Iii)
    print, MIN(JJi), MAX(JJi)
    print, MIN( Wi), MAX( Wi)
    cdataO=FLTARR(NLONS,NLATS)
    cdataO(*,*) = 0.0
        FF=FLTARR(NLONS,NLATS)
        FF(*,*) = 0.0
    for N=0,NT-1 do begin
        IO=IIo(N)-1
        JO=JJo(N)-1
      ; Check if I'm in plot Range
        if ( (lonsO(IO) ge LonBeg) and (lonsO(IO) le LonEnd) AND $
            (latsO(JO) ge LatBeg) and (latsO(JO) le LatEnd) ) then begin
        II=IIi(N)-1
        JI=JJi(N)-1
        VALI = cdataM(II,JI)
        if (VALI ne undef) then begin
          W = Wo(N)
          cdataO(IO,JO) = cdataO(IO,JO) + W * VALI
              FF(IO,JO) =    FF(IO,JO) + W
        endif
        endif
    endfor
    izeros = WHERE(FF eq 0.0, COMPLEMENT=isVALID)
    cdataO(isVALID)=cdataO(isVALID)/FF(isVALID)
    cdataO(izeros) =!VALUES.F_NAN
    cdataM=cdataO
    isCubedSphere=0
    PIXEL=1
    cLons=gLons
    cLats=gLats
endif else begin
    print, 'Regridding Cube to Latlon'
;  Fix missing data
    imiss = WHERE(cdataM eq undef)
    if (imiss(0) ne -1) then cdataM(imiss) = !VALUES.F_NAN
 
; Setup the mapping
    if (~keyword_set(PROJECT)) then begin
      SET_PLOT, 'Z'
      DEVICE, SET_RESOLUTION=[cLons, cLats]
      loadct, 0
      erase, 0
      maplimit=STRING(LatBeg)+STRING(LonBeg)+STRING(LatEnd)+STRING(LonEnd)
      strMap = "MAP_SET, XMARGIN=0, YMARGIN=0, /NOBORDER, LIMIT=["+STRJOIN(STRSPLIT(maplimit,/EXTRACT),',')+"]"
      print, strMap
      res=EXECUTE(strMap)
    endif
    bdataM=BYTSCL(cdataM)
    if (keyword_set(DATA_LEVELS)) then begin
      dlevs=DATA_LEVELS
      nlevs=N_ELEMENTS(dlevs)
      blevs=254L*(FINDGEN(nlevs)/(nlevs-1))+1
      ;help, dlevs, blevs
    endif else begin
      if (keyword_set(RANGE)) then begin
        cMIN=RANGE[0]
        cMAX=RANGE[1]
      endif else begin
        cMIN=MIN(cdataM, /NAN)
        cMAX=MAX(cdataM, /NAN)
      endelse
      ;help, cMIN, cMAX
      dlevs=(FINDGEN(255)/254.0)*(cMAX-cMIN) + cMIN
      if (keyword_set(LOG)) then begin
          logBins=cMIN + [0.345, 0.402, 0.468, 0.545, 0.634, 0.7381, 0.859129, 1, 1.16397, 1.35483, 1.57698, 1.83556, 2.13654, 2.48687, 2.89464, 3.36928, 3.92175, 4.5648, 5.31329, 6.18452, 7.1986, 8.37896, 9.75287, 11.3521, 13.2135, 15.3801, 17.902, 20.8374, 24.2541, 28.2311, 32.8602, 38.2483, 44.5199, 51.8198, 60.3168, 70.207, 81.7189, 95.1185, 110.715, 128.869, 150.0, 174.6, 203.2]*(cMAX-cMIN)/203.2
          XLocs = (N_ELEMENTS(logBins)-1)*FINDGEN(255)/FLOAT(255-1)
          dlevs = INTERPOLATE(logBins, XLocs, /GRID)
          print, dlevs
      endif
      blevs=(FINDGEN(255)+1)
    endelse
    NCOLORS=255
    if (keyword_set(NCOLORS)) then begin
        XLocs = (N_ELEMENTS(dlevs)-1)*FINDGEN(NCOLORS)/FLOAT(NCOLORS-1)
        dlevs = INTERPOLATE(dlevs, XLocs)
        nIncr = FLOAT(N_ELEMENTS(dlevs)-1)/3.0
        dlevs = SMOOTH(dlevs,nIncr)
        blevs=(FINDGEN(NCOLORS)+1)
    endif
    print, dlevs
 
; map the data
if (1) then begin
    iplot=0
  ; CUBE FACES
    for tile=1,6 do begin
 
    inDomain=0
    print, 'Checking for data in Cube Tile '+STRTRIM(STRING(tile),1)
    gridfile='/discover/nobackup/projects/gmao/g6dev/FMS_MOSAICS/mosaic_20090716/CUBED_GRIDS/c48_grid.tile'+STRTRIM(STRING(tile),1)+'.nc'
    if (keyword_set(STRETCH)) then gridfile='/discover/nobackup/projects/gmao/g6dev/FMS_MOSAICS/mosaic_20160301/CUBED_GRIDS/c48_grid_stretch.tile'+STRTRIM(STRING(tile),1)+'.nc'
    fileID = ncdf_open(gridfile, /NOWRITE)
    varID=ncdf_varid(fileID,'x')
    ncdf_varget,fileID,varID,Xverts
    iNeg=WHERE(Xverts gt 180)
    if (iNeg(0) ne -1) then Xverts(iNeg)=Xverts(iNeg)-360
    varID=ncdf_varid(fileID,'y')
    ncdf_varget,fileID,varID,Yverts
    ncdf_close, fileID
    ifind = WHERE( (Xverts ge LonBeg) AND (Xverts le LonEnd) )
    if (ifind(0) ne -1) then inDomain=1
    ifind = WHERE( (Yverts ge LatBeg) AND (Yverts le LatEnd) )
    if (ifind(0) ne -1) then inDomain=inDomain+1
    print, MIN(Xverts), MAX(Xverts)
 
    if (inDomain eq 2) then begin
      gridfile='/discover/nobackup/projects/gmao/g6dev/FMS_MOSAICS/mosaic_20090716/CUBED_GRIDS/c'+STRTRIM(STRING(NX),1)+'_grid.tile'+STRTRIM(STRING(tile),1)+'.nc'
      if (keyword_set(STRETCH)) then gridfile='/discover/nobackup/projects/gmao/g6dev/FMS_MOSAICS/mosaic_20160301/CUBED_GRIDS/c'+STRTRIM(STRING(NX),1)+'_grid_stretch.tile'+STRTRIM(STRING(tile),1)+'.nc'
      fileID = ncdf_open(gridfile, /NOWRITE)
      varID=ncdf_varid(fileID,'x')
      ncdf_varget,fileID,varID,Xverts
      iNeg=WHERE(Xverts gt 180)
      if (iNeg(0) ne -1) then Xverts(iNeg)=Xverts(iNeg)-360
      varID=ncdf_varid(fileID,'y')
      ncdf_varget,fileID,varID,Yverts
    ; help, Xverts, Yverts
      print, MIN(Xverts), MAX(Xverts)
    ; print, MIN(Yverts), MAX(Yverts)
      print, 'Plotting Cube Tiles '+STRTRIM(STRING(tile),1)
        for i = 0, NX-1 do begin
            if ((i mod 288) eq 0) then print, 'Plotting Cube Tiles '+STRTRIM(STRING(tile),1)+' i='+STRTRIM(STRING(i),1)
            for j = 0, NX-1 do begin
              iS=i*2
              jS=j*2
              iC=i
              jC=j + (tile-1)*NX
            ; Check if I'm in plot Range
              if ( (MAX(Xverts(iS:iS+2,jS:jS+2)) ge LonBeg) and (MIN(Xverts(iS:iS+2,jS:jS+2)) le LonEnd) ) then begin
              if ( (MAX(Yverts(iS:iS+2,jS:jS+2)) ge LatBeg) and (MIN(Yverts(iS:iS+2,jS:jS+2)) le LatEnd) ) then begin
              xsym = [Xverts(iS,jS),Xverts(iS,jS+2),Xverts(iS+2,jS+2),Xverts(iS+2,jS)]
              ysym = [Yverts(iS,jS),Yverts(iS,jS+2),Yverts(iS+2,jS+2),Yverts(iS+2,jS)]
              ifind=WHERE(dlevs ge cdataM(iC,jC))
              myColor = ifind(0)+1
              if (FINITE(cdataM(iC,jC)) eq 0) then myColor = 0
 
              if ( (keyword_set(STRETCH)) AND $
                    ((MIN(Yverts(iS:iS+2,jS:jS+2)) gt 85) OR MAX(Yverts(iS:iS+2,jS:jS+2)) lt -85) ) then begin
            ;    SORTING = SORT(xsym)
            ;    ysym = ysym(SORTING)
            ;    xsym = xsym(SORTING)
                if (iplot eq 0) then begin
                ;help, iplot
                ;print, xsym, ysym
                ;polyfill, xsym, ysym, /data, color=myColor
                endif
                iplot=iplot+1
              endif else begin
                polyfill, xsym, ysym, /data, color=myColor
              endelse
            ;  plots,  xsym, ysym, /data, color=myColor, thick=3
 
              endif
              endif
            endfor
        endfor
      ncdf_close, fileID
    endif
    endfor ; tiles
endif else begin
    varID=ncdf_varid(fileID,'CLONS')
    ncdf_varget,fileID,varID,LONS
    varID=ncdf_varid(fileID,'CLATS')
    ncdf_varget,fileID,varID,LATS
    ncdf_close,fileID
    LONS=LONS*(180.0/!DPI)
    LATS=LATS*(180.0/!DPI)
    mySym=SYM(1)
    for i=0,NX-1 do begin
    if ((i mod 288) eq 0) then print, 'Plotting Cube Points i='+STRTRIM(STRING(i),1)
    for j=0,NY-1 do begin
        ifind=WHERE(dlevs ge cdataM(i,j)) + 1
      ;plots,LONS[i,j],LATS[i,j],COLOR=ifind(0),PSYM=3
        plots,LONS[i,j],LATS[i,j],COLOR=ifind(0),PSYM=mySym,SYMSIZE=0.5
    endfor
    endfor
endelse
    cdataM=FLOAT(TVRD())
  ;help, cdataM
    if (0) then begin
    for i=0,cLons-1 do begin
    for j=0,cLats-1 do begin
        ifind=WHERE(blevs ge cdataM(i,j))
        cdataM(i,j) = dlevs(ifind(0))
    endfor
    endfor
    endif else begin
      cdataB=cdataM
      for n=1,N_ELEMENTS(dlevs) do begin
        ifind = WHERE(cdataB ge n)
        if (ifind(0) ne -1) then begin
          cdataM(ifind) = dlevs(n-1)
        endif
      endfor
    endelse
    imiss = WHERE(cdataB eq 0)
    if (imiss(0) ne -1) then cdataM(imiss)=!VALUES.F_NAN
    help, imiss
    ifind = WHERE(cdataM lt dlevs(0))
    cdataM(ifind) = dlevs(0)
    print, min(cdataM,/NAN), max(cdataM,/NAN)
  endelse
  endif
 
 
endif else begin
  print, myVar
  print, myfile
  stop
endelse
imiss = WHERE(cdataM eq undef)
if (imiss(0) ne -1) then cdataM(imiss) = !VALUES.F_NAN
s= SIZE(cdataM)
is3D=0
if (s(0) eq 3) then is3D=1
is2D=0
if (s(0) eq 2) then is2D=1
is1D=0
if (s(0) eq 1) then is1D=1
nLevs=0
if (is3D) then nLevs=s(3)
 
if ( (not isCubedSphere) OR (isCubedSphere and keyword_set(regridCubeToLatlon)) ) then begin
if (~keyword_set(PROJECT)) then begin
; Interpolate
s= SIZE(cdataM)
nLons=DOUBLE(s(1))
nLats=DOUBLE(s(2))
if ( (nLons ne cLons) OR (nLats ne cLats) ) then begin
  if (is3D) then begin
    XLocs = nLons*FINDGEN(cLons)/FLOAT(cLons-1)
    YLocs = nLats*FINDGEN(cLats)/FLOAT(cLats-1)
    ZLocs = FINDGEN(nLevs)
    cdataM = INTERPOLATE(cdataM, XLocs, YLocs, ZLocs, /GRID)
  endif else begin
    if (is2D) then begin
      if (keyword_set(PIXEL)) then begin
      print, 'PIXEL Interpolation'
      XLocs = LONG(nLons*FINDGEN(cLons)/FLOAT(cLons-1))
      YLocs = LONG(nLats*FINDGEN(cLats)/FLOAT(cLats-1))
      cdataM = INTERPOLATE(cdataM, XLocs, YLocs, /GRID)
      endif else begin
      print, 'SMOOTH Interpolation'
      ;if (nLons GT cLons) then begin
      ;  oLons = nLons
      ;  oLats = nLats
      ;  iLons = oLons/2.0
      ;  iLats = (oLats-1)/2.0 + 1.0
      ;  while (iLons ge cLons) do begin 
      ;    help, cdataM
      ;    print, iLons, iLats, cLons, cLats
      ;    XLocs = (oLons*FINDGEN(iLons)/FLOAT(iLons-1))
      ;    YLocs = (oLats*FINDGEN(iLats)/FLOAT(iLats-1))
      ;    cdataM = INTERPOLATE(cdataM, XLocs, YLocs, /GRID)
      ;    oLons = iLons
      ;    oLats = iLats
      ;    iLons = oLons/2.0
      ;    iLats = (oLats-1)/2.0 + 1.0
      ;  endwhile
      ;endif else begin
        XLocs = (nLons*FINDGEN(cLons)/FLOAT(cLons-1))
        YLocs = (nLats*FINDGEN(cLats)/FLOAT(cLats-1))
        cdataM = INTERPOLATE(cdataM, XLocs, YLocs, /GRID)
      ;  help, XLocs, YLocs
      ;  help, cdataM
      ;endelse
      endelse
    endif
  endelse
endif
;help, cdataM
 
; Subdomain
if ( (not isCubedSphere) OR (isCubedSphere and ~keyword_set(regridCubeToLatlon)) ) then begin
if (is3D) then begin
  if (keyword_set(FLIP)) then flip_hemispheres, cdataM
  for k=0,nLevs-1 do begin
    data2D=cdataM(*,*,k)
    get_domain, data2D, LonBeg, LonEnd, LatBeg, LatEnd
    if (k eq 0) then begin
        s=size(data2D)
        sLons=s(1)
        sLats=s(2)
        data3D=FLTARR(sLons,sLats,nLevs)
    endif
    ; Flip the vertical so 0:nLevs-1 goes from surface:top
    data3D(*,*,nLevs-k-1)=data2D
  endfor
  data=cdataM
  data2D=0
  data3D=0
endif else begin
  if (is2D) then begin
    if (keyword_set(FLIP)) then flip_hemispheres, cdataM
    get_domain, cdataM, LonBeg, LonEnd, LatBeg, LatEnd
    data=cdataM
  endif else begin
  if (myVar eq "lon") then begin
      dlon = FLOAT(LonEnd-LonBeg)/FLOAT(cLons)
      data = LonBeg + dlon*FINDGEN(cLons)
  endif
  if (myVar eq "lat") then begin
      dlat = FLOAT(LatEnd-LatBeg)/FLOAT(cLats)
      data = LatBeg + dlat*FINDGEN(cLats)
  endif
  endelse
endelse
cdataM=0
endif else begin
  data=cDataM
endelse
 
endif else begin
  data=cDataM
endelse
 
endif else begin ; PROJECT
  data=cDataM
endelse
 
print, myVar
help, data
print, min(data, /NAN), max(data, /NAN)
 
end
 
<syntaxhighlight lang="idl" line>
; cube_to_latlon.pro
 
== 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:
% cube2latlon merra.nc  latlon.nc  4x5
1) run src/GMAO_Shared/GEOS_Util/plots/configure. This will produce a .quickplotrc in that directory
2) source .quickplotrc
3) run grads.
4) run command dc to plot the cube-sphere grid NetCDF file. Users can get help by just run dc without any argument.
Retrieved from "https://geos5.org/wiki/index.php?title=Visualizing_data_in_Cubed-Sphere_grid"