|
|
| (32 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, input and output 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')
| |
| 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)
| | Then the users can follow the steps. First, we assume a version of GEOS-5 is available. |
| s=SIZE(myVars)
| | # run $ESMADIR/src/GMAO_Shared/GEOS_Util/plots/configure. This will produce a .quickplotrc in that directory |
| if (s(1) eq 1) then begin
| | # source $ESMADIR/src/GMAO_Shared/GEOS_Util/plots/.quickplotrc |
| varID=ncdf_varid(fileID,myVar)
| | # run grads |
| if (keyword_set(LEVEL)) then begin
| | # sdfopen file.nc4 |
| print, 'Reading LEVEL: ', LEVEL
| | # run command <tt>dc</tt> to plot the cube-sphere grid NetCDF file. Users can get help by just run dc without any argument. |
| 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
| | ==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/ |
|
| |
|
| if (keyword_set(TRMM_DATA)) then begin
| | ==Converting (interpolating) cubed-sphere data to Lat-Lon 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 | | 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).) |
| 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)
| | 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 |
| 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
| | Converting data: |
| 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
| | % cube2latlon <in-file> <out-file> <target-resolution> |
| 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
| | E.g., |
| 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 ==
| |
| | |
| 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.
| |