;---To mark the location of the sta-------------------------------------------------------- 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" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl" begin ;**********read data************************************** ;------------- read station datas-------------------------------- filepath_08 = "/home/hy1/datado/obs/201608/201608130800.txt" lines_08 = asciiread(filepath_08,-1,"string") nsta_08 = dimsizes(lines_08) t2_08_sta = stringtofloat(str_get_field(lines_08(0:nsta_08-1),4," ")) lat_08 = stringtofloat(str_get_field(lines_08(0:nsta_08-1),3," ")) lon_08 = stringtofloat(str_get_field(lines_08(0:nsta_08-1),2," ")) ;---------------------------------------------------------------------------- filepath_14 = "/home/hy1/datado/obs/201608/201608131400.txt" lines_14 = asciiread(filepath_14,-1,"string") nsta_14 = dimsizes(lines_14) t2_14_sta = stringtofloat(str_get_field(lines_14(0:nsta_14-1),4," ")) lat_14 = stringtofloat(str_get_field(lines_14(0:nsta_14-1),3," ")) lon_14 = stringtofloat(str_get_field(lines_14(0:nsta_14-1),2," ")) ;---------------------------------------------------------------------------- filepath_20 = "/home/hy1/datado/obs/201608/201608132000.txt" lines_20 = asciiread(filepath_20,-1,"string") nsta_20 = dimsizes(lines_20) t2_20_sta = stringtofloat(str_get_field(lines_20(0:nsta_20-1),4," ")) lat_20 = stringtofloat(str_get_field(lines_20(0:nsta_20-1),3," ")) lon_20 = stringtofloat(str_get_field(lines_20(0:nsta_20-1),2," ")) ;---------------------------------------------------------------------------- ;--------------------read wrfout default data---------------------------------------------------------- wrf_file1 = addfile("/home/hy1/wrf38/WRFV3/test/em_real/wrfout20160811-default/wrfout_d03_2016-08-11_00:00:00.nc","r") t2_08_sim1 = (/wrf_user_getvar(wrf_file1,"T2",48)/) t2_08_sim1 = t2_08_sim1-273.15 t2_14_sim1 = (/wrf_user_getvar(wrf_file1,"T2",54)/) t2_14_sim1 = t2_14_sim1-273.15 t2_20_sim1 = (/wrf_user_getvar(wrf_file1,"T2",60)/) t2_20_sim1 = t2_20_sim1-273.15 ;---------------------read wrfout modify data------------------------------------------------------------------- wrf_file2 = addfile("/home/hy1/wrf38/WRFV3/test/em_real/wrfout20160811/wrfout_d03_2016-08-11_00_00_00.nc","r") t2_08_sim2 = (/wrf_user_getvar(wrf_file2,"T2",48)/) t2_08_sim2 = t2_08_sim2-273.15 t2_14_sim2 = (/wrf_user_getvar(wrf_file2,"T2",54)/) t2_14_sim2 = t2_14_sim2-273.15 t2_20_sim2 = (/wrf_user_getvar(wrf_file2,"T2",60)/) t2_20_sim2 = t2_20_sim2-273.15 printVarSummary(t2_08_sim1) ; exit ;************************************************ wks = gsn_open_wks("ps","/home/hy1/datado/stat") gsn_define_colormap(wks,"BlAqGrYeOrReVi200") ; choose colormap ;------------- Plotting sim for T2--------------------------------------- rest = True rest@gsnDraw = False rest@gsnFrame = False rest@NoHeaderFooter = True rest@Footer = False rest@cnInfoLabelOn = False rest@cnConstFLabelOn = False rest@FieldTitle = " " rest@SubFieldTitle = " " rest@UnitLabel = "" rest@InitTime=False rest@cnFillOn = True rest@cnLevelSelectionMode = "ExplicitLevels" rest@cnLevels = ispan(25,40,1) rest@cnFillColors = (/toint(fspan(2,186,17))/) rest@lbLabelBarOn = False rest@trGridType = "TriangularMesh" shp_name1="/home/qks/local/ncl/lib/ncarg/nclscripts/cnmap/huadong.shp" lnres1=True lnres1@gsLineColor = "grey" lnres1@gsLineThicknessF = 2.0 lnres1@tfPolyDrawOrder = "PostDraw" ;-- Need Change --; ;---Set up some map resources. mpres = True mpres@gsnFrame = False ; Don't advance the frame mpres@gsnDraw = False ; Don't advance the frame mpres@mpOutlineOn = True ; turn off outline mpres@mpFillOn = False mpres@tmYLLabelFontHeightF = 0.012 mpres@tmXBLabelFontHeightF = 0.03 mpres@mpMinLatF = 30.1234 mpres@mpMaxLatF = 32.1005 mpres@mpMinLonF = 120.219 mpres@mpMaxLonF = 122.553 mpres@mpOutlineBoundarySets = "NoBoundaries" mpres@mpAreaMaskingOn = True mpres@mpOceanFillColor="White" mpres@mpLandFillColor="White" mpres@mpInlandWaterFillColor="White" reslay = True plot = new((/3,3/),graphic) map1=wrf_map(wks,wrf_file1, mpres) ;画地图范围 poly1=gsn_add_shapefile_polylines(wks, map1, shp_name1,lnres1) ;画加了区县的地图 p1 = wrf_contour(wrf_file1,wks,t2_08_sim1,rest) ;等值线 plt_res = True plt_res@PanelPlot = True plt_res@FramePlot = False ;-- Need Change --; plt_res@NoTitles = True ;-- Need Change --; plt_res@CommonTitle = False ;-- Need Change --; ; mpres = True ;-- Need Chagne --; mpres@PanelPlot = True mpres@FramePlot = False ;-- Need Change --; over_id1 = wrf_map_overlays(wrf_file1,wks,(/p1/),plt_res,mpres) plot(0,0) = over_id1 ;---------------------------------------------------------------------------------------- map2=wrf_map(wks,wrf_file1, mpres) poly2=gsn_add_shapefile_polylines(wks, map2, shp_name1,lnres1) p2 = wrf_contour(wrf_file1,wks,t2_14_sim1,rest) over_id2 = wrf_map_overlays(wrf_file1,wks,(/p2/),plt_res,mpres) plot(1,0) = over_id2 ;------------------------------------------------------------------------------------------- map3=wrf_map(wks,wrf_file1, mpres) poly3=gsn_add_shapefile_polylines(wks, map3, shp_name1,lnres1) p3 = wrf_contour(wrf_file1,wks,t2_20_sim1,rest) over_id3 = wrf_map_overlays(wrf_file1,wks,(/p3/),plt_res,mpres) plot(2,0) = over_id3 ;------------------------------------------------------------------------------------------- map4=wrf_map(wks,wrf_file1, mpres) poly4=gsn_add_shapefile_polylines(wks, map4, shp_name1,lnres1) p4 = wrf_contour(wrf_file2,wks,t2_08_sim2,rest) over_id4 = wrf_map_overlays(wrf_file1,wks,(/p4/),plt_res,mpres) plot(0,1) = over_id4 ;----------------------------------------------------------------------------- map5=wrf_map(wks,wrf_file1, mpres) poly5=gsn_add_shapefile_polylines(wks, map5, shp_name1,lnres1) p5 = wrf_contour(wrf_file2,wks,t2_14_sim2,rest) over_id5 = wrf_map_overlays(wrf_file1,wks,(/p5/),plt_res,mpres) plot(1,1) = over_id5 ;--------------------------------------------------------------------------- map6=wrf_map(wks,wrf_file1, mpres) poly6=gsn_add_shapefile_polylines(wks, map6, shp_name1,lnres1) p6 = wrf_contour(wrf_file2,wks,t2_20_sim2,rest) over_id6 = wrf_map_overlays(wrf_file1,wks,(/p6/),plt_res,mpres) plot(2,1) = over_id6 ;---------------------------------------------------------------------- ; This function attaches a labelbar to the given plot. ;---------------------------------------------------------------------- function attach_labelbar(wks,map,arr[*]:numeric,colors[*]) local lbres, vph, vpw, nboxes begin getvalues map "vpHeightF" : vph "vpWidthF" : vpw end getvalues nboxes = dimsizes(colors) lbres = True ; labelbar only resources lbres@lbAutoManage = False ; Necessary to control sizes lbres@lbFillColors = colors lbres@vpWidthF = 0.7 * vpw ; labelbar width lbres@vpHeightF = 0.1 * vph ; labelbar height lbres@lbMonoFillPattern = True ; Solid fill pattern lbres@lbOrientation = "horizontal" lbres@lbPerimOn = False lbres@lbLabelAlignment = "InteriorEdges" lbres@lbTitlePosition = "Right" lbres@lbTitleDirection = "Across" lbres@lbTitleString = "~S~o~N~C" lbres@lbBoxLinesOn = False lbres@lbLabelStride = 1 lbres@pmLabelBarHeightF=0.15 lbres@pmLabelBarWidthF=0.85 lbres@lbTopMarginF = 0.2 lbres@lbLeftMarginF = -0.20 lbres@lbRightMarginF = -0.70 lbid = gsn_create_labelbar(wks,nboxes,""+arr,lbres) amres = True amres@amJust = "TopCenter" amres@amParallelPosF = 0.0 ; Center amres@amOrthogonalPosF = 0.55 ; Bottom annoid = gsn_add_annotation(map,lbid,amres) return(annoid) end ;--------------------------------------------------------------------- ;---------------------------------------------------------------------- ; Main code ;---------------------------------------------------------------------- arr = ispan(25,40,1) narr = dimsizes(arr) labels = new(narr+1,string) ; Labels for legend. num_distinct_markers = dimsizes(arr)+1 ; number of distinct markers poly_again = new(3,graphic) ;-- Need Change --; do ishici = 0,2 map=wrf_map(wks,wrf_file1, mpres) poly_again(ishici)=gsn_add_shapefile_polylines(wks, map, shp_name1,lnres1) ;-- Need Change --; ;------------------------------------------------------------------------- gsres = True gsres@gsMarkerIndex = 16 ; Use filled dots for markers. gsres@gsMarkerSizeF = 0.005 colors = tointeger(fspan(2,186,17)) ;------------------------------------------------------------------------ if(ishici.eq.0)then t2_sta = t2_08_sta t2_lat = lat_08 t2_lon = lon_08 end if if(ishici.eq.1)then t2_sta = t2_14_sta t2_lat = lat_14 t2_lon = lon_14 end if if(ishici.eq.2)then t2_sta = t2_20_sta t2_lat = lat_20 t2_lon = lon_20 end if lat_new = new((/num_distinct_markers,dimsizes(t2_sta)/),float,-999) lon_new = new((/num_distinct_markers,dimsizes(t2_sta)/),float,-999) ;---Group the points according to which range they fall in. do i = 0, num_distinct_markers-1 if (i.eq.0) then indexes = ind(t2_sta.lt.arr(0)) end if if (i.eq.num_distinct_markers-1) then indexes = ind(t2_sta.ge.max(arr)) end if if (i.gt.0.and.i.lt.num_distinct_markers-1) then indexes = ind(t2_sta.ge.arr(i-1).and.t2_sta.lt.arr(i)) end if if (.not.any(ismissing(indexes))) then npts_range = dimsizes(indexes) ; # of points in this range. lat_new(i,0:npts_range-1) = t2_lat(indexes) lon_new(i,0:npts_range-1) = t2_lon(indexes) end if delete(indexes) ; Necessary b/c "indexes" may be a different size next time. end do do i = 0, num_distinct_markers-1 if (.not.ismissing(lat_new(i,0))) gsres@gsMarkerColor = colors(i) dumstr = unique_string("marker") map@$dumstr$ = gsn_add_polymarker(wks,map,lon_new(i,:),lat_new(i,:),gsres) end if end do plot(ishici,2) = map delete(map) delete(lat_new) delete(lon_new) delete(t2_sta) delete(t2_lat) delete(t2_lon) end do print("ok") ;---Panel resources pres = True pres@gsnPanelLabelBar = True ; Add color labelbar pres@pmLabelBarHeightF = 0.08 ; Default kind of thin pres@pmLabelBarWidthF = 0.50 pres@gsnPanelRight = 0.90 ; Leave room for shaded labelbar on ; right, to be created below. pres@gsnPanelMainString = " " pres@lbTitlePosition = "Top" pres@lbTitleString = "temperature (~S~o~N~C)" pres@lbTitleFontHeightF = 0.01 pres@lbTopMarginF = 0.2 pres@gsnPanelFigureStrings = (/"a) Case1-08h","b) Case2-08h","c) Obs-08h",\ "d) Case1-14h","e) Case2-14h","f) Obs-14h",\ "g) Case1-20h","h) Case2-20h","i) Obs-20h"/) pres@gsnPanelFigureStringsPerimOn=False pres@amJust = "TopRight" plotm=ndtooned(plot) gsn_panel(wks,plotm(0:8),(/3,3/),pres) map=wrf_map(wks,wrf_file1, mpres) lbid = attach_labelbar(wks,map,arr,colors) ; Attach a labelbar ;draw(plot) ; Drawing the map will draw everything ;frame(wks) ; Advance the frame. end