; Concepts illustrated: ; - Calculate latitude/longitude arrays that define the location ; of a circles around multiple central locations . ; - Extracting variables from a variable of type 'list' ; - Using gc_inout to mask data inside a great circle ;---------------------------------------------------------------------- latmin = 0 latmax = 30 lonmin = 75 lonmax = 100 ;---Multiple locations for cyclone track ;---Read lat/lon data from Track file ; March filename = "track_mar.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 slat=data(:,0) ; Latitude ; lat:storm centers slon=data(:,1) ; Longitude ; lon:storm centers srad = (/ 500 /) ; storm radii (km) srad_unit = 1 ; ; 1=km, 0=degrees N = 180 ; # of points; more points nicer 'circle' opt = False circle = geolocation_circle(slat, slon, srad, srad_unit, N, opt) ; circle is type list ;printVarSummary(circle) ; list containing two variables nLoc = dimsizes(slat) ; # of locations nRad = dimsizes(srad) ; # of circles at each location circle_lat = circle[0] ; For clarity: explicitly extract list elements circle_lon = circle[1] ;printVarSummary(circle_lat) ;printVarSummary(circle_lon) ;---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) clat = new((/15,15,180/),float) ; arrays to hold great circle clon= new((/15,15,180/),float) do i =0,nLoc-1 do j = 0,nRad-1 clat(i,j,:) = circle_lat(i,j,:) clon(i,j,:) = circle_lon(i,j,:) clat1 = ndtooned(clat) clon1 = ndtooned(clon) ;printVarSummary(clat) ;printVarSummary(clon) ;printVarSummary(clat1) min_lat = min(clat1) min_lon = min(clon1) max_lat = max(clat1) max_lon = max(clon1) ;---Subset the desired rectagle of data R1 := tcavgmar({min_lat:max_lat},{min_lon:max_lon}) conR1 := marcont({min_lat:max_lat},{min_lon:max_lon}) ;---Set points that are outside of the circle of data to missing lat2d := conform(R1,R1&latitude,0) lon2d := conform(R1,R1&longitude,1) in_circle := gc_inout(lat2d,lon2d,clat1,clon1) R1 = where(in_circle,R1,R1@_FillValue) lat2d := conform(conR1,conR1&latitude,0) lon2d := conform(conR1,conR1&longitude,1) in_circle := gc_inout(lat2d,lon2d,clat1,clon1) conR1 = where(in_circle,conR1,conR1@_FillValue) end do end do ;printVarSummary(R1) ;printVarSummary(conR1) ;printMinMax(R1,0) ;printMinMax(conR1,0) ;---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@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 ;; 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,R1,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,conR1,res2) gsn_panel(wks,plot,(/1,2/),resP)