;---------------------------------------------------------------------- ; Concepte 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 ; data boundaries for desired region latmax = 25 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) 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 = 90 ; # of points; more points nicer 'circle' but more time opt = False circle = geolocation_circle(slat, slon, srad, srad_unit, N, opt) ; circle is type list printVarSummary(circle) ; list containing two variables print("-----") 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) ; [location | 15] x [radii | 1] x [circle | N] printVarSummary(circle_lon) print("-----") ;---Read in TRMM precipitation data ;---March 2000 ;; double pre(time, longitude, latitude) ; ;; pre:long_name = "precipitation" ; ;; pre:units = "mm" ; ;; pre:_FillValue = NaN ; <=== not recognized by NCL ;; pre:missing_value = NaN ; <=== ?not ieee? f2 = addfile("all_mar.nc", "r") prc2 = f2->pre(:,{lonmin:lonmax},{latmin:latmax}) print(any(isnan_ieee(prc2))) ;;; if (any(isnan_ieee(prc2))) then ; ?not" IEEE NaN? value1 = 1d20 ;;; replace_ieeenan (prc2, value, 0) ;;; prc2@missing_value = prc2@_FillValue ; missing_value not used ;;; end if delete( [/prc2@_FillValue, prc2@missing_value/] ) ; ?not? IEEE Nan prc2@_FillValue = value1 prc2@missing_value = prc2@_FillValue printVarSummary(prc2) printMinMax(prc2,0) print("-----") prc2_reorder = prc2(time|:,latitude|:,longitude|:) printVarSummary(prc2_reorder) ; [time | 31] x [latitude | 120] x [longitude | 100] printMinMax(prc2_reorder,0) print("-----") ymd = cd_calendar(prc2_reorder&time, -2) print(ymd) ; March 1-31, 2000 print("-----") ;---Explore data: Look at the data distribution of the 1st day opt_sd = True opt_sd@PrintStat = True nt= 0 ; first ?day" stat_sd = stat_dispersion(prc2_reorder(nt,:,:), opt_sd ) print("-----") ;; prc2_reorder@_FillValue = 0.0 ; not a good _FillValue prc2_reorder@_FillValue = value1 ; Extract March 27-30 ntim2 = (/26,27,28,29/) ; March 27-30 tcmar = prc2_reorder(ntim2,:,:) tcmar@long_name = "March: "+tcmar@long_name printVarSummary(tcmar) ; [time | 4] x [latitude | 120] x [longitude | 100] printMinMax(tcmar,0) print("-----") ;---Explore Data: Look at the data distribution of days March 27-30 opt_sd = True opt_sd@PrintStat = True stat_sd = stat_dispersion(tcmar, opt_sd ) print("-----") tcavgmar = dim_avg_n_Wrap(tcmar,0) tcavgmar@long_name = "March 27-30: Average precip" tcavgmar@units = prc2_reorder@units printVarSummary(tcavgmar) ; [latitude | 120] x [longitude | 100] printMinMax(tcavgmar,0) print("-----") ;--- March contribution totmar = dim_sum_n_Wrap(prc2_reorder,0) ; [latitude | 120] x [longitude | 100] totmar@long_name = "Total March precip" totmar@units = prc2_reorder@units printVarSummary(totmar) printMinMax(totmar,0) print("-----") martc = dim_sum_n_Wrap(tcmar,0) ; [latitude | 120] x [longitude | 100] martc@long_name = "March 27-30: Total precip" printVarSummary(martc) printMinMax(martc,0) print("-----") nonzerototmar = where(totmar.eq.0.0, prc2_reorder@_FillValue, totmar) copy_VarMeta(totmar,nonzerototmar) printVarSummary(nonzerototmar) printMinMax(nonzerototmar,0) print("-----") marcont = (martc / nonzerototmar)*100 copy_VarMeta(martc,marcont) printVarSummary(marcont) printMinMax(marcont,0) print("-----") ;---Select a 'nice' contour level range for all plots ;---Method 1: It will give a common contour spacing mnmxint_R1 = nice_mnmxintvl( min(tcavgmar({latmin:latmax},{lonmin:lonmax})) \ , max(tcavgmar({latmin:latmax},{lonmin:lonmax})), 18, False) mnmxint_conR1 = nice_mnmxintvl( min(marcont({latmin:latmax},{lonmin:lonmax})) \ , max(marcont({latmin:latmax},{lonmin:lonmax})), 18, False) if (.not.all(mnmxint_R1.eq.mnmxint_conR1)) then print(mnmxint_R1) print("----") print(mnmxint_conR1) print("----") print("Change cnLevels for each array!") end if ;---Graphic resources pltType = "png" ; "x11", "png", "pdf", "eps" pltDir = "./" pltName = "tc500_modi_mar" pltPath = pltDir + pltName wks = gsn_open_wks(pltType, pltPath) res = True res@gsnDraw = False res@gsnFrame = False res@gsnAddCyclic = False ; extracted region is not cyclic in longitude ;res@gsnStringFontHeightF = 0.015 ; make subtitles smaller res@mpMinLonF = lonmin res@mpMaxLonF = lonmax res@mpMinLatF = latmin res@mpMaxLatF = latmax res@mpCenterLonF = (lonmax+lonmin)/2 res@mpDataBaseVersion = "MediumRes" ; higher map resolution res@pmTickMarkDisplayMode = "Always" ; nicer tickmarks for regional plots res@cnFillOn = True res@cnLinesOn = False res@cnLineLabelsOn = False res@cnFillMode = "RasterFill" ;;res@cnFillPalette = "WhiteBlueGreenYellowRed" ;;res@cnFillPalette = "WhBlGrYeRe" res@cnFillPalette = "Example" ;---Manually: specify contour levels ;---Method 2: It will allow unequal contour spacing res@cnLevelSelectionMode = "ExplicitLevels" ; "mm/day res@cnLevels = (/1,5,10,15,20,25,30,40,50,60,70,80,90,100/) res@lbLabelBarOn = False ; turn off individual label bars resP = True ; panel resP@gsnMaximize = True ; Maximize Panel resP@gsnPaperOrientation = "landscape" resP@gsnPanelLabelBar = True ; add common colorbar ;resP@lbLabelFontHeightF = 0.007 ; make labels smaller ;resP@gsnPanelRowSpec = True ; tell panel what order to plot ;---For each center and radius: extract extent of each circle [ reduce gc_inout time ] ;---Use NCL's 'reassignment syntax [ := ] to accomodate possiblr changing array sizes cdim = dimsizes(circle_lat) nCenter = cdim(0) ; # of centers: same as: nCenter=dimsizes(slat) nRadius = cdim(1) ; # of circles: same as: nRadius=dimsizes(srad) nPoints = cdim(2) ; # of points : same as: nPoints=N ;plot = new ( (/nCenter,2/), "graphic") plot = new (2, "graphic") do i=0,nLoc-1 do j=0,nRad-1 clat = circle_lat(i,j,:) clon = circle_lon(i,j,:) ;printVarSummary(clat) ;printVarSummary(clon) ;printVarSummary(clat) min_clat = min(clat) min_clon = min(clon) max_clat = max(clat) max_clon = max(clon) ;---Subset the desired rectagle of data R1 := tcavgmar({min_clat:max_clat},{min_clon:max_clon}) conR1 := marcont({min_clat:max_clat},{min_clon:max_clon}) ;---Set points that are outside of the circle of data to missing lat2d_R1 := conform(R1,R1&latitude,0) lon2d_R1 := conform(R1,R1&longitude,1) in_circle := gc_inout(lat2d_R1,lon2d_R1,clat,clon) R1 = where(in_circle, R1, R1@_FillValue) conR1 = where(in_circle, conR1, R1@_FillValue) lat2d_conR1 := conform(conR1,conR1&latitude,0) lon2d_conR1 := conform(conR1,conR1&longitude,1) resP@gsnPanelMainString = "March 27-30, 2000: " \ + "("+sprintf("%4.1f", slat(i))+","+sprintf("%4.1f", slon(i))+")" plot(0) = gsn_csm_contour_map(wks, R1,res) plot(1) = gsn_csm_contour_map(wks,conR1,res) gsn_panel(wks,plot,(/1,2/),resP) end do ; end 'j' [Radii] end do ; end 'i' [nLoc]