GEOS-5

Visualizing data in Cubed-Sphere grid: Difference between revisions

Newer edit →
Created page with "#==== 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..."
 
Line 13: Line 13:
== IDL ==
== IDL ==


Dr. Bill Putman has created many IDL programs to visualize GEOS-5’s outputs. Users may find some useful tools at  /discover/nobackup/projects/gmao/osse2/IDL3 . The program read_and_interpolate_cube.pro can read data in the cube-sphere grid in NetCDF file, interpolate the data , map the data in Lat-Lon format and plot. The program cube_to_latlon.pro can use read_and_interpolate_cube.pro to convert the variables in cube-sphere grid to lat-lon grid and write the variable to NetCDF files.
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 ==
== GrAds ==
Retrieved from "https://geos5.org/wiki/index.php?title=Visualizing_data_in_Cubed-Sphere_grid"