;calcuate trend of annual mean temperature begin ;united states mask file and average april 1st swe > 10 mm during historical period g = addfile("/stor9000/apps/users/NWSUAF/2017050898/ZXY/code/shp/usa_mask.nc", "r") ; united states USA_mask = g->USA_mask ; 1 we need this area .eq. 1 mask1 = conform_dims((/1008,224,200/), USA_mask(:,0:199), (/1,2/)) h = addfile("/stor9000/apps/users/NWSUAF/2017050898/ZXY/code/swe_maks.nc", "r") ;april 1st swe > 10 mm mask_swe = h->mask_swe ; 1 we need this area .eq. 1 mask2 = conform_dims((/1008,224,200/), mask_swe, (/1,2/)) condtnl_expr = (mask1 .eq. 1) .and. (mask2 .eq. 1) wks = gsn_open_wks("pdf", "annual_trend_tas") plot = new(6,graphic) model = (/"GFDL-ESM2G","IPSL-CM5A-MR","MIROC5","GFDL-ESM2G","IPSL-CM5A-MR","MIROC5"/) dirname = (/"GF_RCP45","IPSL_RCP45","MIROC_RCP45","GF_RCP85","IPSL_RCP85","MIROC_RCP85"/) do np = 0, 5 ;---------------------------process the data------------------------------------------ fils = systemfunc("ls /stor9000/apps/users/NWSUAF/2017050898/ZXY/data/monthly/"+dirname(np)+"/*.nc") f = addfiles(fils, "r") ListSetType(f, "cat") lat = f[0]->lat lon = f[0]->lon time=ispan(0, 1007, 1) lon!0 = "lon" lon@long_name = "lon" lon@units = "degrees_east" lon&lon = lon lat!0 = "lat" lat@long_name = "lat" lat@units = "degrees_north" lat&lat = lat time!0 = "time" time@long_name= "time" time@units = "month since 2016-01" time&time = time TAIR = where(condtnl_expr,tofloat(f[:]->TAIR),-999.9) TAIR!0 = "time" TAIR&time = time TAIR!1 = "lat" TAIR&lat = lat TAIR!2 = "lon" TAIR&lon = lon TAIR@missing_value = -999.9 TAIR@_FillValue = -999.9 tair = month_to_annual(TAIR, 1) printVarSummary(tair) var = tair(lat|:,lon|:,year|:) year = ispan(0, 83, 1) year!0 = "year" year@long_name= "year" year@units = "year since 2016" year&year = year trend = regCoef(year,var)*10.0 ;degc/decade pq = -(trend_manken(tair, False, 0)-1) ;significance level sig = pq(0,:,:) copy_VarCoords(TAIR(0,:,:), trend) printVarSummary(trend) copy_VarCoords(TAIR(0,:,:), sig) ;------------------------plot---------------------------------------------- res = True res@gsnAddCyclic = False res@gsnDraw = False res@gsnFrame = False res@gsnMaximize = True res@gsnLeftString = tostring(model(np)) res@gsnRightString = "(~F34~0~F~/decade)" res@mpProjection = "CylindricalEquidistant" res@mpDataBaseVersion = "MediumRes" res@mpLimitMode = "LatLon" res@mpMaxLatF = 50 res@mpMinLatF = 25 res@mpMaxLonF = -100 res@mpMinLonF = -125 res@mpPerimOn = True res@mpPerimLineColor = "Black" res@mpPerimDrawOrder = "PostDraw" res@mpFillOn = True ; turn off map fill res@mpDataBaseVersion = "MediumRes" res@mpDataSetName = "Earth..4" res@mpOutlineBoundarySets = "National" res@mpNationalLineThicknessF = 0.3 ;set national boundary line thickness res@mpProjection = "CylindricalEquidistant" res@mpOutlineDrawOrder = "PostDraw" ; draw continental outline last res@mpMaskAreaSpecifiers = (/"Conterminous US"/) res@mpLandFillColor = "white" res@mpInlandWaterFillColor = "white" res@mpOceanFillColor = "white" res@cnFillOn = True res@cnLinesOn = False res@cnLineLabelsOn = False res@cnLineLabelsOn = False res@cnFillPalette = "MPL_YlOrBr" res@cnFillDrawOrder = "PreDraw" ; important order res@cnFillMode = "RasterFill" res@cnMissingValFillPattern = "SolidFill" res@cnMissingValFillColor = "white" res@cnLevelSelectionMode = "ManualLevels" res@cnLevelSpacingF = 0.2 res@cnMinLevelValF = 0 res@cnMaxLevelValF = 2 ; res@lbOrientation = "Vertical" res@lbLabelBarOn = False plotA = gsn_contour_map(wks, trend, res) sres = True sres@gsnDraw = False sres@gsnFrame = False sres@cnLineLabelsOn = False sres@cnLinesOn = False sres@cnInfoLabelOn = False sres@cnFillOn = True sres@lbLabelBarOn= False sres@cnLevelSelectionMode = "ExplicitLevels" ;sres@cnFillDrawOrder="postDraw" sres@cnLevels = (/0.05/) ; set to significance level sres@cnMonoFillPattern = False ; want multiple patterns sres@cnFillPatterns = (/17,-1/) ; the patterns sres@cnMonoFillScale = False ; want different densities sres@cnFillScales = (/1.0,1.0/) ; change densities sres@cnMonoFillColor =True sres@cnFillDotSizeF = 0.00003 plotB = gsn_csm_contour(wks, sig, sres) overlay(plotA, plotB) plot(np) = plotA delete(plotA) delete(plotB) end do ;add shapefile polylines;;;;;;;;; dirin = "/stor9000/apps/users/NWSUAF/2017050898/ZXY/code/shp/shp_wus_huc2/" lnres = True lnres@gsLineColor = "Black" lnres@gsLineThicknessF = 0.01 shpid1 = new(6, graphic) shpid2 = new(6, graphic) shpid3 = new(6, graphic) shpid4 = new(6, graphic) shpid5 = new(6, graphic) shpid6 = new(6, graphic) shpid7 = new(6, graphic) shpid8 = new(6, graphic) shp1_name = dirin+"WBDHUC2-10.shp" shp2_name = dirin+"WBDHUC2-11.shp" shp3_name = dirin+"WBDHUC2-13.shp" shp4_name = dirin+"WBDHUC2-14.shp" shp5_name = dirin+"WBDHUC2-15.shp" shp6_name = dirin+"WBDHUC2-16.shp" shp7_name = dirin+"WBDHUC2-17.shp" shp8_name = dirin+"WBDHUC2-18.shp" shpid1 = gsn_add_shapefile_polylines(wks, plot, shp1_name, lnres) shpid2 = gsn_add_shapefile_polylines(wks, plot, shp2_name, lnres) shpid3 = gsn_add_shapefile_polylines(wks, plot, shp3_name, lnres) shpid4 = gsn_add_shapefile_polylines(wks, plot, shp4_name, lnres) shpid5 = gsn_add_shapefile_polylines(wks, plot, shp5_name, lnres) shpid6 = gsn_add_shapefile_polylines(wks, plot, shp6_name, lnres) shpid7 = gsn_add_shapefile_polylines(wks, plot, shp7_name, lnres) shpid8 = gsn_add_shapefile_polylines(wks, plot, shp8_name, lnres) resP = True resP@gsnMaximize = True resP@gsnPanelLabelBar = True resP@lbLabelFontHeightF = 0.01 gsn_panel(wks, plot, (/2,3/), resP) end