;---------------------------------------------------------------------- ; Concepts illustrated: ; - Using nggcog to create a great circle ; - Using gc_inout to mask data inside a great circle ;---------------------------------------------------------------------- ; This script was contributed by Jake Huff, a Masters student in the ; Climate Extremes Modeling Group at Stony Brook University. ; ; It shows how to use nggcog to generate a great circle, and gc_inout ; to mask the data inside the circle. ;---------------------------------------------------------------------- begin latmin = 0 latmax = 30 lonmin = 75 lonmax = 100 ;---Read lat/lon data from Track file ; March filename = "martrack.txt" ncols = numAsciiCol(filename) ; Calculate the number of columns data = readAsciiTable(filename,ncols,"float",0) ; Given the # of columns, we can use readAsciiTable to read this file lat=data(:,0) ; Latitude lon=data(:,1) ; Longitude ;---Set all values outside the lat/lon box of interest to missing. ;value = -9999.9 ;lat@_FillValue = value ;lon@_FillValue = value ;lat = where(lat .ge. latmin .and. lat .le. latmax,lat,lat@_FillValue) ;lon = where(lon .ge. lonmin .and. lon .le. lonmax,lon,lon@_FillValue) ;---Read in TRMM precipitation data ;March f2 = addfile("all_mar.nc", "r") prc2 =f2->pre(:,{75:100},{0:30}) prc2_reorder = prc2(time|:,latitude|:,longitude|:) value1 = 1d20 prc2_reorder@missing_value = value1 ; Force change to CF recognized value prc2_reorder@_FillValue = value1 ; March ntim2 = (/26,27,28,29/) tcmar = prc2_reorder(ntim2,:,:) tcavgmar = dim_avg_n_Wrap(tcmar,0) ; March contribution totmar = dim_sum_n_Wrap(prc2_reorder,0) martc = dim_sum_n_Wrap(tcmar,0) nonzerototmar = where(totmar.eq.0.0, prc2_reorder@_FillValue, totmar) copy_VarMeta(totmar,nonzerototmar) marcont = (martc / nonzerototmar)*100 copy_VarMeta(martc,marcont) ;---Start the graphics plot = new(2,graphic) wks_type = "png" wks_type@wkWidth = 2500 wks_type@wkHeight = 2500 wks = gsn_open_wks(wks_type,"tc500km_mar1") res = True res@gsnMaximize = True res@gsnAddCyclic = False res@mpLimitMode = "LatLon" ; select subregion res@mpMinLonF = lonmin res@mpMaxLonF = lonmax res@mpMinLatF = latmin res@mpMaxLatF = latmax ;res@mpCenterLonF = (lonmin+lonmax)/2 res@mpProjection = "Mercator" res@mpDataBaseVersion = "HighRes" ; use GMT coastline res@pmTickMarkDisplayMode = "Always" ; nicer tickmarks for regional plots res@cnFillOn = True res@cnLinesOn = False res@cnInfoLabelOn = False ; turn off contour info label res@cnLinesOn = False ; turn off contour lines res@cnLineLabelsOn = False ; turn off line labels res@mpFillOn = False ; turn off gray continents res@lbLabelBarOn = True ; add single label bar res@lbOrientation = "horizontal" ;res@cnFillMode = "RasterFill" res@cnFillPalette = "precip3_16lev" res@trGridType = "TriangularMesh" res@gsnMaximize = True res@gsnPaperOrientation = "auto" res@gsnLeftString = " " res@gsnRightString = " " res@mpShapeMode = "FreeAspect" res@vpWidthF = 0.9 res@vpHeightF = 0.7 res@tmYROn = False ; turn off right and top tickmarks res@tmXTOn = False ;---Get valid lat/lon values for cyclone ii = ind(.not.ismissing(lat(0:1)).and..not.ismissing(lon(0:1))) nind = dimsizes(ii) ; Loop through non-missing lat/lon tracks for the cyclone and create a ; filled contour plot within the calculated lat/lon circle. clat = new(2,float) ; arrays to hold great circle clon = new(2,float) radius = 5 do ni=0,nind-1 j = ii(ni) lat_location = lat(j) lon_location = lon(j) nggcog(lat_location,lon_location,radius,clat,clon) ; Calculate great circle min_lat = min(clat) min_lon = min(clon) max_lat = max(clat) max_lon = max(clon) ;---Subset the desired rectagle of data martrmm := tcavgmar({min_lat:max_lat},{min_lon:max_lon}) contmar := marcont({min_lat:max_lat},{min_lon:max_lon}) ;---Set points that are outside of the circle of data to missing lat2d := conform(martrmm,martrmm&latitude,0) lon2d := conform(martrmm,martrmm&longitude,1) in_circle := gc_inout(lat2d,lon2d,clat,clon) martrmm = where(in_circle,martrmm,martrmm@_FillValue) lat2d := conform(contmar,contmar&latitude,0) lon2d := conform(contmar,contmar&longitude,1) in_circle := gc_inout(lat2d,lon2d,clat,clon) contmar = where(in_circle,contmar,contmar@_FillValue) end do ;; Cyclone Tracking ;March filename4 = "00_bob01.txt" ncols = numAsciiCol(filename4) data4= readAsciiTable(filename4,ncols,"float",0) y4=data4(:,0) x4=data4(:,1) ; polyline resources resp = True ; mods yes resp@gsLineThicknessF = 12.0 ; line thickness resp@gsLineColor = "black" ;Add markers to the trajectories. amres = True ; marker resources for best track amres@gsMarkerIndex = 16 ; marker style (filled circle) amres@gsMarkerSizeF = 6.0 ; marker size amres@gsMarkerColor = "red" ; maker color ;creat panel resources resP = True resP@gsnMaximize = True resP@gsnPanelLabelBar =False ; No common label bar ;resP@gsnPanelMainString = "Yearly TCs mean rainfall" ;---Create a plot res@gsnLeftString = "e) Mean TCs rainfall (mm): March" res@cnLevelSelectionMode = "ExplicitLevels" ; set explicit contour levels res@cnLevels = (/1,5,10,20,30,40,50,60,70,80/) plot(0) = gsn_csm_contour_map(wks,martrmm,res) d1 = gsn_add_polyline(wks,plot(0),x4,y4,resp) d2 = gsn_add_polymarker(wks,plot(0),x4(0),y4(0),amres) res2 = True res2@gsnMaximize = True res2@gsnAddCyclic = False res2@mpLimitMode = "LatLon" ; select subregion res2@mpMinLonF = lonmin res2@mpMaxLonF = lonmax res2@mpMinLatF = latmin res2@mpMaxLatF = latmax res2@mpProjection = "Mercator" res2@mpDataBaseVersion = "Highres" ; use GMT coastline res2@pmTickMarkDisplayMode = "Always" ; nicer tickmarks for regional plots res2@cnFillOn = True res2@cnLinesOn = False res2@cnInfoLabelOn = False ; turn off contour info label res2@cnLinesOn = False ; turn off contour lines res2@cnLineLabelsOn = False ; turn off line labels res2@mpFillOn = False ; turn off gray continents res2@lbLabelBarOn = True ; add single label bar res2@lbOrientation = "horizontal" ;res2@cnFillMode = "RasterFill" res2@cnFillPalette = "precip3_16lev" res2@trGridType = "TriangularMesh" res2@gsnMaximize = True res2@gsnPaperOrientation = "auto" res2@gsnLeftString = " " res2@gsnRightString = " " res2@mpShapeMode = "FreeAspect" res2@vpWidthF = 0.9 res2@vpHeightF = 0.7 res2@tmYROn = False ; turn off right and top tickmarks res2@tmXTOn = False res2@gsnLeftString = "f) TCs contribution (%): March" res2@cnLevelSelectionMode = "ExplicitLevels" ; set explicit contour levels res2@cnLevels = (/1,5,10,15,20,25,30,40,50,60,70,80,90/) plot(1) = gsn_csm_contour_map(wks,contmar,res2) gsn_panel(wks,plot,(/1,2/),resP) end