add China map 求助

1649518749

推荐：利用shp文件画底图
但是我叠加中国地图是只有国界，没有省界。求大神看看。方法来源http://bbs.06climate.com/forum.php?mod=viewthread&tid=11797
代码如下：
;   Sample script to create and plot PW for WRF-ARW model output
;   November 2009
load "$NCARG_ROOT/lib/ncarg/nclscripts/cnmap/cnmap.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
begin
;
; The WRF ARW input file.
; This needs to have a ".nc" appended, so just do it.
  a = addfile("/home/fc/datafile/wrfoutdata/wrfout_d01_2016-08-17_18:00:00","r")

; We generate plots, but what kind do we prefer?
;  type = "x11"
type = "ps"
; type = "ps"
; type = "ncgm"
  wks = gsn_open_wks(type,"pw")
;  gsn_define_colormap(wks,"gsdtol")

; Set some basic resources
  res = True
  res@MainTitle                   = "REAL-TIME WRF"
;  res@Footer = False
  pltres = True
  mpres = True
;;;;;;;;;;;;;;;;;set for map;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
res@mpMaxLatF= 55.

res@mpMinLatF= 15.

res@mpMaxLonF= 136.

res@mpMinLonF= 72.
mpres@mpFillOn                = True
mpres@mpOutlineOn             = True  ; Use outlines from shapefile
res@cnFillDrawOrder         = "PreDraw"
mpres@mpDataBaseVersion       = "MediumRes"
mpres@mpDataSetName           = "Earth..4"
mpres@mpAreaMaskingOn         = True
mpres@mpMaskAreaSpecifiers    = (/"China","Taiwan","Disputed area between India and China","India:Arunachal Pradesh"/)
mpres@mpLandFillColor         = "white"
mpres@mpInlandWaterFillColor  = "white"
mpres@mpOceanFillColor        = "white"
mpres@mpOutlineBoundarySets   = "NoBoundaries"  ;or National
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  gas_const = 287. ; J/K/kg
  Cp        = 1004. ; J/K/kg

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; What times and how many time steps are in this data set?
  times = wrf_user_getvar(a,"times",-1)  ; get 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

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; First get the variables we will need        

    pres = wrf_user_getvar(a,"pres",it)                ; pressure in pa
    ph = wrf_user_getvar(a,"PH",it)                    ; pert. geopt
    phb = wrf_user_getvar(a,"PHB",it)                  ; base geopt
    th = wrf_user_getvar(a,"theta",it)                 ; theta (T+300)
    qvapor = wrf_user_getvar(a,"QVAPOR",it)   

; Calc height, pressure and virtual temperature
    height = (ph + phb)/9.81                           ; height at full levels
    temp = th * (pres/100000) ^ (gas_const/Cp)
    vtemp = (1 + 0.61*qvapor) * temp                   ; virtual temp

    zdiff=height(0,:,:)
    zdiff=0.
    pw_sfc_ptop=height(0,:,:)
    pw_sfc_ptop=0.

    dim = dimsizes(pres)
    do k=0,dim(0)-1
      zdiff(:,:)= (height(k+1,:,:) - height(k,:,:))
      pw_sfc_ptop(:,:) =  pw_sfc_ptop(:,:) + ((pres(k,:,:)/(gas_const * vtemp(k,:,:))) * qvapor(k,:,:) * zdiff(:,:))
    end do
    pw_sfc_ptop@description = "PW"

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

    opts = res                          
    opts@cnFillOn = True

    ;opts@lbTitleOn = False
    opts@lbLabelBarOn = True
    ;opts@lbLabelsOn = False
    ;opts@pmLabelBarDisplayMode = "NoCreate"
    contour  = wrf_contour(a,wks,pw_sfc_ptop,opts)
    delete(opts)

    ; MAKE PLOTS                                       
    plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;>============================================================<
;                      add China map
;>------------------------------------------------------------<
cnres           = True
cnres@china     = True       ;draw china map or not
cnres@river     = True       ;draw changjiang&huanghe or not
cnres@province  = True       ;draw province boundary or not
cnres@nanhai    = True       ;draw nanhai or not
cnres@diqu      = False       ; draw diqujie or not

chinamap = add_china_map(wks,plot,cnres)
;>============================================================<
  end do  ;       end of the time loop

end

1649518749

这幅图，国界都没有。

xyan88

我画的和你的一样 也没有地图

1649518749

xyan88 发表于 2016-8-19 11:15
我画的和你的一样 也没有地图

如果不是用wrfout数据，画出来倒是有省界。

andrewsoong

1649518749 发表于 2016-8-19 11:16
如果不是用wrfout数据，画出来倒是有省界。

一次画多张就不行

1649518749

andrewsoong 发表于 2016-8-19 11:22
一次画多张就不行

没用啊，我还叠加了海平面气压场，还是没国界

xyan88

1649518749 发表于 2016-8-19 11:16
如果不是用wrfout数据，画出来倒是有省界。

请问那怎么解决啊

1649518749

xyan88 发表于 2016-8-19 16:14
请问那怎么解决啊

我都茫然了。你贴个脚本来看看，一起研究研究

xyan88

好的 我在测试下  好了贴上来

xyan88

1649518749 发表于 2016-8-19 16:18
我都茫然了。你贴个脚本来看看，一起研究研究

;   Sample script to create and plot PW for WRF-ARW model output
;   November 2009
load "$NCARG_ROOT/lib/ncarg/nclscripts/cnmap/cnmap.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
begin
;

  a=addfile("wrfout_d01_2015-01-01_00:00:00.nc","r")
  wks = gsn_open_wks("pdf","0819")

; Set some basic resources
res = True
mpres = True
pltres= True

res@MainTitle = "Real-TIME WRF-Chem"
res@cnFillDrawOrder = "PreDraw"
times = wrf_user_getvar(a,"times",-1)  ; get times in the file
res@cnFillOn = True
res@TimeLabel = times(0)      ; Set Valid time to use on plots

;;;;;;;;;;;;;;;;;set for map;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
mpres@mpMaxLatF= 25.
mpres@mpMinLatF= 15.
mpres@mpMaxLonF= 136.
mpres@mpMinLonF= 72.
mpres@mpFillOn                = True
mpres@mpOutlineOn             = True  ; Use outlines from shapefill
mpres@mpDataBaseVersion       = "MediumRes"
mpres@mpDataSetName           = "Earth..4"
mpres@mpAreaMaskingOn         = True
mpres@mpMaskAreaSpecifiers    = (/"China","Taiwan","Disputed area between India and China","India:Arunachal Pradesh"/)
;mpres@mpMaskAreaSpecifiers    = (/"Hunan","Anhui"/)
mpres@mpLandFillColor         = "white"
mpres@mpInlandWaterFillColor  = "white"
mpres@mpOceanFillColor        = "white"
mpres@mpOutlineBoundarySets   = "NoBoundaries"  ;or National

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
data = a->PM10(0,0,:,:)
contour  = wrf_contour(a,wks,data,res)
plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;>============================================================<
;                      add China map
;>------------------------------------------------------------<
cnres           = True
cnres@china     = True       ;draw china map or not
cnres@river     = True       ;draw changjiang&huanghe or not
cnres@province  = True       ;draw province boundary or not
cnres@nanhai    = True       ;draw nanhai or not
cnres@diqu      = False       ; draw diqujie or not

chinamap =  add_china_map(wks,plot,cnres)
;>============================================================<

end