登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
其中红色区域就是添加shapefile的内容参见官网 ; Example script to produce plots for a WRFreal-data run, ; with the ARW coordinate dynamics option. 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/cnmap/cnmap.ncl" load"$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl" begin ; TheWRF ARW input file. ; Thisneeds to have a ".nc" appended, so just do it. a =addfile("/home/fc/datafile/wrfoutdata/wrfout_d01_2016-08-17_18:00:00","r") ; Wegenerate plots, but what kind do we prefer? ; type ="x11" ; type ="pdf" type = "ps" ; type ="ncgm" wks = gsn_open_wks(type,"surface") ; Set somebasic resources res = True ;res@MainTitle = "REAL-TIME WRF" pltres = True pltres@PanelPlot = True ;这句话不写,图片经纬度会出问题,这句话找了好久 mpres = True ;;;;;;;;;;;;;;;;;setfor map;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; mpres@pmTickMarkDisplayMode= "Always" ;显示经纬度 mpres@mpLimitMode = "LatLon" mpres@mpMinLatF= 15 mpres@mpMaxLatF= 55 mpres@mpMinLonF= 72 mpres@mpMaxLonF= 136 mpres@mpFillOn = True mpres@mpOutlineOn = False ; Use outlines from shapefile mpres@mpDataBaseVersion = "MediumRes" mpres@mpDataSetName = "Earth..4" mpres@mpAreaMaskingOn = True mpres@mpMaskAreaSpecifiers =(/"China","Taiwan","Disputed area between India andChina","India:Arunachal Pradesh"/) mpres@mpLandFillColor = "white" mpres@mpInlandWaterFillColor = "white" mpres@mpOceanFillColor = "white" mpres@mpOutlineBoundarySets = "NoBoundaries" ;or National ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Whattimes and how many time steps are in the data set? times =wrf_user_getvar(a,"times",-1) ; get all times in the file ntimes = dimsizes(times) ; number of times in the file ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; do it = 0,ntimes-1,2 ; TIME LOOP print("Working on time: " +times(it) ) res@TimeLabel = times(it) ; Set Valid time to use on plots ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Firstget the variables we will need slp =wrf_user_getvar(a,"slp",it) ; slp wrf_smooth_2d( slp, 3 ) ; smooth slp tc =wrf_user_getvar(a,"tc",it) ; 3D tc td =wrf_user_getvar(a,"td",it) ; 3D td u =wrf_user_getvar(a,"ua",it) ; 3D U at mass points v =wrf_user_getvar(a,"va",it) ; 3D V at mass points td2 = wrf_user_getvar(a,"td2",it) ; Td2 in C tc2 =wrf_user_getvar(a,"T2",it) ; T2 in Kelvin tc2 = tc2-273.16 ; T2 in C u10 =wrf_user_getvar(a,"U10",it) ; u at 10 m, mass point v10 =wrf_user_getvar(a,"V10",it) ; v at 10 m, mass point ; tf2 = 1.8*tc2+32. ; Turn temperature intoFahrenheit ; tc2@description = "SurfaceTemperature" ; tc2@units = "~S~o~N~C" ; td_f = 1.8*td2+32. ; Turn temperature intoFahrenheit ; td2@description = "Surface Dew Point Temp" ; td2@units = "~S~o~N~C" ; u10 = u10*1.94386 ; Turn wind into knots ; v10 =v10*1.94386 ; u10@units = "kts" ; v10@units = "kts" ; dt=tc2-td2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Plotting options for T opts = res opts@cnFillOn = True opts@ContourParameters = (/ 5., 40., 5./) opts@cnFillColors =(/"cadetblue1","steelblue2","green","green4",\ "gold2","orange","red","red3","red4"/) opts@gsnSpreadColorEnd = -3 ; End third from the last color in color map contour_tc = wrf_contour(a,wks,tc2,opts) delete(opts) ; MAKE PLOTS plot =wrf_map_overlays(a,wks,(/contour_tc/),pltres,mpres) ;================add china map ==================================; ;**************************************************************************** ; add shapefiles ;**************************************************************************** ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;国界 shp1="/home/fc/ncl/lib/ncarg/nclscripts/chinamap/bou1_4l.shp" lnres1 = True lnres1@gsLineColor ="black" lnres1@gsLineThicknessF = 2.0 ; 2x thickness shp_plot1 = gsn_add_shapefile_polylines(wks,plot,shp1,lnres1) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;省级界线 shp2="/home/fc/ncl/lib/ncarg/nclscripts/chinamap/bou2_4p.shp" lnres2 = True lnres2@gsLineColor ="black" lnres2@gsLineThicknessF = 1.5 shp_plot2 = gsn_add_shapefile_polylines(wks,plot,shp2,lnres2) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;市级界线 shp3="/home/fc/ncl/lib/ncarg/nclscripts/chinamap/diquJie_polyline.shp" lnres3 = True lnres3@gsLineColor ="black" lnres3@gsLineThicknessF = 1.0 shp_plot3 =gsn_add_shapefile_polylines(wks,plot,shp3,lnres3) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;河流 ;shp4="/home/fc/ncl/lib/ncarg/nclscripts/chinamap/hyd1_4l.shp" ; lnres4 = True ; lnres4@gsLineColor = "blue" ;lnres4@gsLineThicknessF = 1.0 ; shp_plot4 =gsn_add_shapefile_polylines(wks,plot,shp4,lnres4) maximize_output(wks,False) ;这句话不写,add shapefiles 就没意义! ;================add china map ==================================; end do ; END OF TIME LOOP end
|