; Concepts illustrated: ; - Reading TRMM netCDF file ; - Change NaN _FillValue to NCL a CF _FillValue ; - Use NCL's distributed 1x1 land-sea mask top create a mask file ; CAVEAT: landsea_mask uses a 1x1 basemap. TRMM 3B42 is 0.25 ; - Writing mask array to a NetCDF for use in other scripts ; - Use the mask to plot over land and lakes only ;*************************************************************** ; double pre ( time, longitude, latitude ) ; long_name : precipitation ; units : mm ; _FillValue : nan <=== There are no Nan (nan) in the file ; missing_value : nan ;*************************************************************** netCDF = False ; True => create netCDF file for subsequent use ; False=> create mask for loacl use PLOT_MASK = False PLOT_DATA = True ;--- landsea.nc is distributed with NCL ;--- CAVEAT: landsea.nc has 1.0 degree resolution ;--- TRMM: 0.25 degree resolution ;--- Read TRMM grid ;--- Must reorder ;--- NaN (nan) is used for _FillValue and missing_value f = addfile("yearsum2000.nc","r") TRMM = f->pre ; TRMM(time,longitude,latitude) TRMM_reorder = TRMM(time|:,latitude|:,longitude|:) ; (time,latitude,longitude) ; type double value = 1d20 TRMM_reorder@missing_value = value ; Force change to CF recognized value TRMM_reorder@_FillValue = value printVarSummary(TRMM_reorder) printMinMax(TRMM_reorder,0) print("---") ;--- Look at data distribution: use to specify plot contour intervals opt = True opt@PrintStat = True statTRMM = stat_dispersion(TRMM_reorder, opt ) print("---") ;--- landsea.nc is distributed with NCL ;--- read 1x1 land sea mask basemap file ;--- create mask array; NOTE basemap is 1 deg while TRMM is 0.25 deg ;--- make mask contain byte values for smaller file f_lsm = addfile("$NCARG_ROOT/lib/ncarg/data/cdf/landsea.nc","r") lsmTRMM = tobyte( landsea_mask(f_lsm->LSMASK,TRMM_reorder&latitude,TRMM_reorder&longitude) ) lsmTRMM@long_name = "Land-Sea Mask: TRMM 3B42" ;printVarSummary(lsmTRMM) copy_VarCoords(TRMM_reorder(0,:,:), lsmTRMM) ;printMinMax(lsmTRMM,0) ;print("---") ;====================== netCDF ============================================== if (netCDF) then dirNC = "./" filNC = "TRMM_3B42.LandSeaMask.nc" pthNC = dirNC + filNC system("/bin/rm -f "+pthNC) ; remove any pre-existing file ncdf = addfile(pthNC, "c") ; open output netCDF file fAtt = True ; assign file attributes fAtt@title = "NCL 1-degree basemap onto 0.25 deg TRMM 3B42 grid" fAtt@Conventions = "None" fAtt@creation_date = systemfunc ("date") fileattdef( ncdf, fAtt ) ; copy file attributes, Defines global (file) attributes associated with a supported file. ncdf->LSMASK_3B42 = lsmTRMM end if ;======================= PLOT =============================================== if (PLOT_MASK) then wks = gsn_open_wks("png","3B42_mask") ; send graphics to PNG file res = True res@mpFillOn = False ; do not color-fill the map res@mpPerimDrawOrder = "PostDraw" ; draw the map perimeter last res@mpMaxLatF = max(TRMM_reorder&latitude) res@mpMinLatF = min(TRMM_reorder&latitude) res@cnFillOn = True ; turn on color fill res@cnLinesOn = False ; turn off the contour lines res@cnLineLabelsOn = False ; turn off the contour line labels res@cnFillMode = "RasterFill" ; use raster fill res@cnRasterSmoothingOn = True res@cnFillPalette = "wh-bl-gr-ye-re" res@cnLevelSelectionMode = "ExplicitLevels" ; explicitly set the levels via cnLevels res@cnLevels = (/1.,2.,3.,4./) ; set the levels res@cnFillColors = (/60,100,20,140,5/) ; set the colors that will be used to color fill res@lbLabelStrings = ispan(0,4,1) ; labels for the labelbar boxes res@lbLabelAlignment = "BoxCenters" ; put the labels in the center of the label bar boxes res@lbTitleString = "0=ocean, 1=land, 2=lake, 3=small island, 4=ice shelf" ; labelbar title res@lbTitleFontHeightF = 0.0125 ; labelbar title font height res@pmLabelBarOrthogonalPosF = .125 ; move whole thing down res@gsnRightString = "CAVEAT: 1 deg basemap used: TRMM 0.25 deg" plt = gsn_csm_contour_map(wks,lsmTRMM,res) end if ; PLOT_MASK if (PLOT_DATA) then wks = gsn_open_wks("png","3B42_year2000") ; send graphics to PNG file plot = new(3,graphic) res2 = True res2@gsnDraw = False ; do not draw the plot res2@gsnFrame = False ; do not advance the frame res2@mpFillOn = False ; do not color-fill the map res2@mpPerimOn = True ; turn the map perimeter on res2@mpPerimDrawOrder = "PostDraw" ; draw the map perimeter last res2@mpMaxLatF = max(TRMM_reorder&latitude) res2@mpMinLatF = min(TRMM_reorder&latitude) res2@cnLinesOn = False ; turn off the contour lines res2@cnLineLabelsOn = False ; turn off the contour line labels res2@cnFillOn = True ; turn on color fill res2@cnLevelSelectionMode = "ExplicitLevels" ;;(0) [2] Min=0 ; from stat_dispersion ;;(0) [3] LowDec=116.459997182712 ; 100 ;;(0) [4] LowOct=159.1799952182919 ; 150 ;;(0) [5] LowSex=235.1699936967343 ; 225 ;;(0) [6] LowQuartile=395.992788657546 ; 400 ;;(0) [7] LowTri=545.1599856931716 ; 550 ;;(0) [8] Median=834.2699758410454 ; 800 ;;(0) [9] HighTri=1181.339965449646 ; 1200 ;;(0) [10] HighQuartile=1418.129985593259 ; 1400 ;;(0) [11] HighSex=1758.269965138286 ; 1750 ;;(0) [12] HighOct=2022.599935388193 ; 2000 ;;(0) [13] HighDec=2251.303572073579 ; 2250 ==> 2500 ;; ; 3000 res2@cnLevels = (/100,150,225,400,550,800,1200,1400,1750,2000,2500,3000/) ; "mm total" res2@cnFillPalette = (/"Azure2","PaleTurquoise","PaleGreen" \ ,"SeaGreen3" ,"Yellow", "Orange","HotPink" \ ; contour colors ,"Red","Violet", "Purple", "Magenta", "Brown","White" /); one more color than contour levels ;; ,"Red","Violet", "Purple", "Magenta", "Brown","Black" /); one more color than contour levels res2@cnFillMode = "RasterFill" ; use raster fill res2@cnMissingValFillPattern = 0 res2@cnMissingValFillColor = "Foreground" ; "White" ; "Background" res2@lbLabelBarOn = False ; turn off individual cb's nt = 0 ; 1st timestep plot(0) = gsn_csm_contour_map(wks,TRMM_reorder(nt,:,:),res2) ; apply the Land-Sea-Mask prcMASK = where(lsmTRMM.eq.1 .or. lsmTRMM.eq.2, TRMM_reorder(nt,:,:), TRMM_reorder@_FillValue) ; Land only copy_VarMeta(TRMM_reorder(nt,:,:), prcMASK) res2@gsnCenterString = "After Land Only Mask" plot(1) = gsn_csm_contour_map(wks,prcMASK,res2) prcMASK = where(lsmTRMM.eq.0, TRMM_reorder(nt,:,:), TRMM_reorder@_FillValue) ; Ocean only res2@gsnCenterString = "After Ocean Only Mask" plot(2) = gsn_csm_contour_map(wks,prcMASK,res2) panres = True ; create a panel resource list panres@gsnMaximize = True ; maximize the size of the paneled plots panres@gsnPanelYWhiteSpacePercent = 1.0 ; leave space between the 2 plots panres@gsnPanelLabelBar = True ; add common colorbar ;panres@lbLabelFontHeightF = 0.007 ; make labels smaller gsn_panel(wks,plot,(/3,1/),panres) end if ; PLOT_DATA