|
|
| Line 17: |
Line 17: |
| 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, input and output files are hard-wired in the code. The users will need to provide those files and change the code accordingly. | | 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, input and output files are hard-wired in the code. The users will need to provide those files and change the code accordingly. |
|
| |
|
| <syntaxhighlight lang="idl" line>
| | ==== [[Recipe: IDL program to read cubed-sphere data| read_and_interpolate_cube.pro]] ==== |
| ; read_and_interpolate_cube.pro
| | ==== [[Recipe: IDL program to convert cubed-sphere data to lat-lon data | cube_to_latlon.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>
| |
| | |
| <syntaxhighlight lang="idl" line>
| |
| ; cube_to_latlon.pro
| |
| pro cube_to_latlon
| |
| | |
| set_plot, 'Z'
| |
| | |
| LonBeg = -180
| |
| LonEnd = 180
| |
| LatBeg = -90
| |
| LatEnd = 90
| |
| lonCen = 0.5*(LonBeg+LonEnd)
| |
| LatCen = 0.0
| |
| undef = 1.e15
| |
| pngImgIdim=1440*4.0
| |
| pngImgJdim=pngImgIdim*0.5 + 1
| |
| boundaryImageObj = get_map_image_new(LatBeg, LatEnd, LonBeg, LonEnd, LonCen, LatCen, $
| |
| /NOBORDER, HORIZ_DIMS=[pngImgIdim, pngImgJdim])
| |
| elevfile='/discover/nobackup/projects/gmao/osse2/stage/c1440_NR/DATA/0.0625_deg/const/const_2d_asm_Nx/Y2005/M11/D01/c1440_NR.const_2d_asm_Nx.20051101.nc4'
| |
| read_and_interpolate, elevfile, 'lon', lons, pngImgIdim, pngImgJdim, LonBeg, LonEnd, LatBeg, LatEnd, undef
| |
| read_and_interpolate, elevfile, 'lat', lats, pngImgIdim, pngImgJdim, LonBeg, LonEnd, LatBeg, LatEnd, undef
| |
| myfile = "/discover/nobackup/projects/gmao/osse2/G5FCST/forecast-20121024_150000-C3072-4390102/G5FCST.inst1_2d_asm_Nx.20121026_1200z.nc4"
| |
| ;myfile = "/discover/nobackup/projects/gmao/osse2/c5760_NR/scratch-3550451.out/c5760_NR.inst10mn_2d_met1_Nx.20120616_0840z.nc4"
| |
| ;myfile = "/discover/nobackup/projects/gmao/osse2/H10UC90R72/scratch-3861154/H10UC90R72.inst30mn_2d_met1_Nx.20150215_0000z.nc4"
| |
| read_and_interpolate_cube, myfile, 'U250', u250, pngImgIdim, pngImgJdim, LonBeg, LonEnd, LatBeg, LatEnd, undef, /regridCubeToLatlon
| |
| read_and_interpolate_cube, myfile, 'V250', v250, pngImgIdim, pngImgJdim, LonBeg, LonEnd, LatBeg, LatEnd, undef, /regridCubeToLatlon
| |
| | |
| ; write netcdf file
| |
| ; Create a new NetCDF file with the filename inquire.nc:
| |
| id = NCDF_CREATE('c3072_NR_U250-V250.nc', /CLOBBER)
| |
| ; Fill the file with default values:
| |
| NCDF_CONTROL, id, /FILL
| |
| xid = NCDF_DIMDEF(id, 'x', pngImgIdim)
| |
| yid = NCDF_DIMDEF(id, 'y', pngImgJdim)
| |
| tid = NCDF_DIMDEF(id, 't', 1)
| |
| ; Define variables:
| |
| latid = NCDF_VARDEF(id, 'lat', [yid], /FLOAT)
| |
| lonid = NCDF_VARDEF(id, 'lon', [xid], /FLOAT)
| |
| timeid = NCDF_VARDEF(id, 'time', [tid], /FLOAT)
| |
| uid = NCDF_VARDEF(id, 'U250', [xid,yid], /FLOAT)
| |
| vid = NCDF_VARDEF(id, 'V250', [xid,yid], /FLOAT)
| |
| NCDF_CONTROL, id, /ENDEF
| |
| ; Input data:
| |
| NCDF_VARPUT, id, timeid, 0
| |
| NCDF_VARPUT, id, lonid, lons
| |
| NCDF_VARPUT, id, latid, lats
| |
| NCDF_VARPUT, id, uid, u250
| |
| NCDF_VARPUT, id, vid, v250
| |
| NCDF_CLOSE, id ; Close the NetCDF file.
| |
| | |
| end
| |
| | |
| </syntaxhighlight>
| |
|
| |
|
| == GrADS == | | == GrADS == |