load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" begin season = "ANN" time0 = (/1951,1981,1998/) time2 = 2012 filename = systemfunc ("ls " + "./READER_SAT/*.nc") ; the txt files' names were listed in the A-Z word order npts = dimsizes(filename) lon = new((/npts/),"float") lat = new((/npts/),"float") R = new((/npts/),"float") pb = R print(npts + " stations") out_filename = "Trend_" + npts + "stations_" + season wks = gsn_open_wks("ps", out_filename) gsn_define_colormap(wks,"WhViBlGrYeOrRe") ; define a different colormap. nc1 = NhlNewColor(wks,.8,.8,.8) ; Add light gray to colormap, ; for continents. map = new(3,"graphic") dum = new(3,"graphic") ; ********************** Loop 1 ********* do j = 0,2 ; *********************************** time1= time0(j) do i = 0, dimsizes(filename)-1 fi = addfile( filename(i), "r") x = fi->sat lon(i) = x@lon lat(i) = x@lat time = x&time year1 = floattoint(time(0)) year2 = floattoint(time(dimsizes(time)-1)) year = ispan(year1,year2,1) delete(time) ;print("from " + year1 + " to " + year2) ; annual mean yann = month_to_annual(x,1) ; annual mean yann!0 = "time" yann&time = year ; calculate the linear trend if(year1.le.time1) then time = ispan(time1, time2, 1) rc = regline( time, yann({time1:time2}) ) else time = ispan(year1, time2, 1) rc = regline( time, yann({year1:time2}) ) end if df = rc@nptxy-2 prob = (1 - betainc(df/(df+rc@tval^2), df/2.0, 0.5) ) ;yReg = rc*time + rc@yintercept ; NCL array notation ; yReg is same length as x, y R(i) = rc pb(i) = prob delete(x) delete(time) delete(year) delete(yann) end do ; ************** plot ******************************** ; -------Options-------- arr = (/-0.1,-0.05,0.,0.05,0.1/) ; bin settings (bin0 = < 0., ; bin1 = 0.:4.999, etc.) colors = (/"blue","blue","blue","red","red","red"/) ; marker colors, dimsizes must ; be equal to dimsizes(arr)+1 sizes = (/3,2,1,1,2,3/)*0.01 labels = new(dimsizes(arr)+1,string) ; Labels for legend. ; ------------------------------ ; Create X and Y arrays to hold the points for each range and initialize ; them to missing values. We want to use num_distinct_markers ; different colors, so we need num_distinct_markers sets of X and ; Y points. num_distinct_markers = dimsizes(arr)+1 ; number of distinct markers lat_new = new((/num_distinct_markers,dimsizes(R)/),float,-999) lon_new = new((/num_distinct_markers,dimsizes(R)/),float,-999) ;printVarSummary(lat_new) ; Group the points according to which range they fall in. At the ; same time, create the label that we will use later in the legend. do i = 0, num_distinct_markers-1 if (i.eq.0) then indexes = ind(R.lt.arr(0)) labels(i) = "x < " + arr(0) end if if (i.eq.num_distinct_markers-1) then indexes = ind(R.ge.max(arr)) labels(i) = "x >= " + max(arr) end if if (i.gt.0.and.i.lt.num_distinct_markers-1) then indexes = ind(R.ge.arr(i-1).and.R.lt.arr(i)) labels(i) = arr(i-1) + " <= x < " + arr(i) end if ; Now that we have the set of indexes whose values fall within ; the given range, take the corresponding lat/lon values and store ; them, so later we can color this set of markers with the appropriate ; color. if (.not.any(ismissing(indexes))) then npts_range = dimsizes(indexes) ; # of points in this range. lat_new(i,0:npts_range-1) = lat(indexes) lon_new(i,0:npts_range-1) = lon(indexes) end if delete(indexes) ; Necessary b/c "indexes" may be a different ; size next time. end do ;printVarSummary(lat_new) ;=========================================================================== ; Begin plotting section. ; ;wks = gsn_open_wks("eps",out_filename ) ; Open a workstation and ;gsn_define_colormap(wks,"WhViBlGrYeOrRe") ; define a different colormap. ;nc1 = NhlNewColor(wks,.8,.8,.8) ; Add light gray to colormap, ; for continents. ; Set up some map resources ; mpres = True mpres@gsnMaximize = True ; Maximize plot in frame. mpres@gsnFrame = False ; Don't advance the frame mpres@gsnDraw = False mpres@mpFillOn = False ; ; Zoom in on United States. ; ; mpres@mpMinLatF = 25. ; mpres@mpMaxLatF = 50. ; mpres@mpMinLonF = 235. ; mpres@mpMaxLonF = 290. mpres@gsnPolar = "SH" mpres@mpMaxLatF = -50 ; mpres@mpFillColors = (/-1,-1,nc1,-1/) ;assign light gray to land masses mpres@tiMainString = "Trend: " + time1 + "-" + time2 + ", " + season map(j) = gsn_csm_map_polar(wks,mpres) ; Create logical variables to hold the marker and text resources. ; These markers are different than the XY markers, because they are not ; associated with an XY plot. You can put these markers on any plot. ; gsres = True gsres@gsMarkerIndex = 16 ; Use filled dots for markers. txres = True txres@txFontHeightF = 0.015 ; ; Loop through each grouping of markers, and draw them one set at ; a time, assigning the proper color and size with gsn_marker. ; ; At the same time, draw a legend showing the meaning of the markers. ; xleg = (/0.07,0.07,0.32,0.32,0.57,0.57,0.82,0.82/) + 0.15 ; Location of xtxt = (/0.16,0.16,0.44,0.44,0.66,0.66,0.91,0.91/) + 0.15 ; legend markers yleg = (/0.22,0.17,0.22,0.17,0.22,0.17,0.22,0.17/) ;- 0.05 ;0.03 ; and text ytxt = (/0.22,0.17,0.22,0.17,0.22,0.17,0.22,0.17/) ;- 0.05 ;0.03 ; strings. ; Loop do i = 0, num_distinct_markers-1 print(i+ " " + lat_new(i,0)) if (.not.ismissing(lat_new(i,0))) gsres@gsMarkerColor = colors(i) gsres@gsMarkerSizeF = sizes(i) print(" " + sizes(i) + " " + colors(i)) ;gsn_polymarker(wks,map(j),lon_new(i,:),lat_new(i,:),gsres) dum(j) = gsn_add_polymarker(wks,map(j),lon_new(i,:),lat_new(i,:),gsres) end if ; Add marker and text for the legend. gsres@gsMarkerColor = colors(i) gsres@gsMarkerSizeF = sizes(i) gsn_polymarker_ndc(wks, xleg(i),yleg(i),gsres) gsn_text_ndc (wks,labels(i),xtxt(i),ytxt(i),txres) end do ; Loop Over ; add significant circle gres = True ;gres@gsnDraw = False ;res@gsnFrame = False ;do i = 0,npts-1 ; print("pb = " + pb(i)) ; if(pb(i).ge.0.95) then ; gres@gsMarkerIndex = 6 ; gres@gsMarkerSizeF = 0.03 ; gres@gsMarkerColor = "black" ; gsn_polymarker(wks, map(j), lon(i),lat(i),gres) ; else ; end if ;end do end do resp = True resp@gsnMaximize = True gsn_panel(wks,map, (/1,3/), resp) ;draw(map) ;frame(wks) end