Visualizing data in Cubed-Sphere grid

From GEOS-5
Jump to navigation Jump to search
  1. ==== 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

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.

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 are hard-wired in the code. The users will need to provide those files and change the code accordingly.

<syntaxhighlight lang="idl" line>

read_and_interpolate_cube.pro

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 

print, myfile
help, myVar
help, LonBeg, LonEnd, LatBeg, LatEnd

if (~keyword_set(TIME)) then TIME=0

exists = FILE_TEST(myfile) if (exists) then begin 

 if ( keyword_set(INFILEID)) then fileID = INFILEID
 if (~keyword_set(INFILEID)) then fileID = ncdf_open(myfile, /NOWRITE)

 if (keyword_set(LISTVARS)) then begin
   ;find the number of variables
   fileinq_struct = ncdf_inquire(fileID)
   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
 dim_unlimited = finq.recdim           ; = -1 if undefined, otherwise index into dim array
 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
  if (NZ ge 1) then begin
   varID=ncdf_varid(fileID,'lev')
   ncdf_varget,fileID,varID,levs
  endif
 endif

 if (myVar eq 'WSPD50M') then begin

    varID=ncdf_varid(fileID,'U50M')
    ncdf_varget,fileID,varID,u
    varID=ncdf_varid(fileID,'V50M')
    ncdf_varget,fileID,varID,v
    cdataM = SQRT(u*U + v*v)

 endif else begin

 myVars=STRSPLIT(myVar,'*', /extract)
 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

 if (keyword_set(TRMM_DATA)) then begin
   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)
 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
     print, 'SCALING by: ', SCALE_FACTOR
     cdataM = cdataM*SCALE_FACTOR
 endif
Check if data is on cubed-sphere grid
 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: 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&oldid=3579"