请问我要怎么改才能使山西的风场去掉但是地形图依旧在呢？

pkn789

请问我想去掉山西的风场，在底图上加上地形高度，我把山西填成了白色，然后把风场叠加在了上面，这样就去掉了山西的风场，但是我在叠加地形图的时候，如果先画地形，把风场叠加在上面，那么山西的地形就也会被覆盖成白色，就会如下图所示，如果我先画风场，那么地形的填色图就会覆盖风场的矢量图，请问我要怎么改才能使山西的风场去掉但是地形图依旧在呢？
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/wrf/WRFUserARW.ncl"
load "/home/wangl/public/cnmap/cnmap.ncl"
load "/home/wangl/public/cnmap/shapefile_mask_data.ncl"
begin
dir = "/home/wangl/pkn/WRFV3/test/em_real/"
filepath = systemfunc("ls " + dir + "*.grib2")
nn = dimsizes(filepath)
print(nn)
do ii = 0,nn-1
print(filepath(ii))
fnl = addfile(filepath(ii),"r")
str = str_get_cols(filepath(ii),6,16)
print(str)
lat = fnl->lat_0(::-1)
lon = fnl->lon_0
u = fnl->UGRD_P0_L100_GLL0({500*100},:,:)
v = fnl->VGRD_P0_L100_GLL0({500*100},:,:)


u = u(::-1,:)
v = v(::-1,:)
printVarSummary(u)
u = u*2.5
v = v*2.5
spd     = (u*u + v*v)^(0.5) ; m/sec
      spd@description = "Wind Speed"
      spd@units = "m/s"
printVarSummary(spd)
klon = fspan(0,359,1801)
        klat = fspan(-90,90,901)
u!0 = "lat"
u&lat = lat
        u@units = "degrees_north"
        u!1 = "lon"
        u@units = "degrees_east"
u&lon = lon
        v!0 = "lat"
        v&lat = lat
        v@units = "degrees_north"
        v!1 = "lon"
        v@units = "degrees_east"
        v&lon = lon

uu = linint2(lon,lat,u,True,klon,klat,0)
vv = linint2(lon,lat,v,True,klon,klat,0)

uu!0 = "lat"
uu&lat = klat
uu!1 = "lon"
uu&lon = klon
uu&lon@units = "degrees_east"
        uu&lat@units = "degrees_north"
vv!0 = "lat"
vv&lat = klat
vv!1 = "lon"
vv&lon = klon
vv&lon@units = "degrees_east"
        vv&lat@units = "degrees_north"
delete(lat)
delete(lon)
undef("create_map")  
function create_map(wks,title)  
local a, res2,minlat,maxlat,minlon,maxlon  
begin
  minlat = 33.8
  maxlat = 38.3
  minlon = 109.9
  maxlon = 116.5  
res               = True
res@tmXTOn = False
res@tmYROn = False
res@mpOutlineOn   = True  
res@mpFillOn      = False  
  res@mpDataBaseVersion = "MediumRes"  
  res@mpDataSetName="Earth..4"  
  res@mpOutlineSpecifiers=(/"China","Shanxi"/)   
  res@mpProvincialLineColor="white"   
  res@mpProvincialLineThicknessF =4
res@mpLimitMode           = "LatLon"  
  res@mpMinLatF             = minlat  
res@mpMaxLatF             = maxlat  
  res@mpMinLonF             = minlon  
  res@mpMaxLonF             = maxlon  
res@mpProjection = "LambertConformal"
res@mpLambertMeridianF  = 113
res@mpLimitMode = "Corners"
res@mpLambertParallel1F = 30.
res@mpLambertParallel2F = 60.
res@mpLeftCornerLatF       = 33.8
   res@mpRightCornerLatF      = 38.3
   res@mpLeftCornerLonF       = 109.9
   res@mpRightCornerLonF      = 116.5
res@pmTickMarkDisplayMode = "Always"

map = gsn_csm_map(wks,res)  
return(map)
end
wks = gsn_open_wks("pdf","3070900700")
map = create_map(wks,"Shanxi of China")  
filename = "/home/wangl/public/map/bou2_4p.shp"   
f = addfile(filename, "r")   ; Open shapefile  
NAME=(/f->NAME/)  
asciiwrite ("NAME.txt", NAME)
  shanxi=(/NAME(25)/)
; Read data off shapefile  
  
  geometry = f->geometry  
  segments = f->segments  
  geomDims = dimsizes(geometry)  
  segsDims = dimsizes(segments)  
   
  geom_segIndex = f@geom_segIndex  
  geom_numSegs  = f@geom_numSegs  
  segs_xyzIndex = f@segs_xyzIndex  
  segs_numPnts  = f@segs_numPnts  
  
  lines       = new(segsDims(0),graphic)   ; Array to hold polygons  
  numFeatures = geomDims(0)  
  
lon   = f->x  
  lat   = f->y  
plres             = True       ; resources for polylines  
  plres@gsEdgesOn   = True       ; draw border around polygons  
plres@gsEdgeColor = "white"
segNum = 0  
  do i=0, numFeatures-1   
  
  
     if( NAME(i).eq. shanxi) then  
         plres@gsFillColor = "white"   
  
         startSegment = geometry(i, geom_segIndex);  
         numSegments  = geometry(i, geom_numSegs)  
         do seg=startSegment, startSegment+numSegments-1  
            startPT = segments(seg, segs_xyzIndex)  
            endPT = startPT + segments(seg, segs_numPnts) - 1  
            lines(segNum) = gsn_add_polygon(wks, map, lon(startPT:endPT),  \  
                                                       lat(startPT:endPT), plres)  
            segNum = segNum + 1  
         end do  
           
      end if      
  end do  
draw(map)

undef("read_elev_data")
function read_elev_data(topo_file)
local nlat, nlon, topo_file, lat, lon
begin
;---Read data as a straight binary file
  nlat = 2160
  nlon = 4320
  setfileoption("bin","ReadByteOrder","BigEndian")
  elev = cbinread(topo_file,(/nlat,nlon/),"short")
;---Create 1D coordinate arrays
  lat       = fspan(90,-90,nlat)
  lon       = fspan(0,360,nlon)
  lat!0     = "lat"
  lon!0     = "lon"
  lat@units = "degrees_north"
  lon@units = "degrees_east"
  lat&lat   = lat
  lon&lon   = lon
;---Attach the coordinate arrays
  elev!0    = "lat"
  elev!1    = "lon"
  elev&lat  = lat
  elev&lon  = lon

  return(elev)
end
elev = read_elev_data("ETOPO5.DAT")


gsn_define_colormap(wks,"BlRe")

res = True
res@tmBorderThicknessF = 3.0
res@pmTickMarkDisplayMode = "Always"
res@tmXTOn = False
res@tmYROn = False
res@gsnPaperOrientation = "Landscape"
res@gsnLeftString = "Wind (m/s) 0900(UTC)  700Hpa "
res@mpOutlineOn = False
res@mpOutlineBoundarySets       = "GeophysicalAndUSStates"
  res@mpOutlineSpecifiers         = (/"China","Shanxi"/)
  res@mpDataSetName               = "Earth..4"   ;  
  res@mpDataBaseVersion           = "MediumRes"  
  res@gsFillColor = "white"
res@mpProjection = "LambertConformal"
res@mpLambertMeridianF  = 113
res@mpLimitMode = "Corners"
res@mpLambertParallel1F = 30.
res@mpLambertParallel2F = 60.
res@mpLeftCornerLatF       = 33.8
   res@mpRightCornerLatF      = 38.3
   res@mpLeftCornerLonF       = 109.9
   res@mpRightCornerLonF      = 116.5


  
  
res@cnLevelSelectionMode = "ManualLevels"
res@cnMinLevelValF   =0           ; set min contour level
res@cnMaxLevelValF   =2000

  res@cnFillOn           = True             ; turn on contour fill
  res@cnLevelSpacingF    = 200              ; NCL picks 2000
  res@cnFillMode         = "RasterFill"     ; much faster than AreaFill
  res@cnLinesOn          = False            ; turn off contour lines
  res@cnLineLabelsOn     = False            ; turn off line labels
  res@cnInfoLabelOn      = False            ; turn off info label
  res@lbBoxLinesOn       = False            ; turn off labelbar box lines
  res@gsnAddCyclic       = False            ; don't add longitude cyclic point
  res@mpFillOn           = False            ; turn off map fill
  res@tiMainString       = "ETOPO5.DAT"     ; main title
  res@pmLabelBarWidthF   = 0.8              ; default is too s

res_vc = True
        res_vcgsnAddCyclic = False
        res_vc@vcGlyphStyle = "WindBarb"
        res_vc@vcMinDistanceF = 0.012
        res_vc@vcWindBarbLineThicknessF = 2
        res_vc@vcLineArrowThicknessF = 2
        res_vc@vcMonoWindBarbColor =True
        res_vc@vcWindBarbColor="Yellow"
        res_vc@vcRefLengthF = 0.02
        res_vc@vcWindBarbTickLengthF = 0.4
        res_vc@vcRefAnnoOn = True
        
        res_vc@vcRefAnnoSide = "Top"
        res_vc@vcRefAnnoString2On = False
        res_vc@vcRefAnnoPerimOn = False
res_vc@mpOutlineSpecifiers         = (/"China","Shanxi"/)
res_vc@gsFillColor = "white"
plot_wind = gsn_csm_vector(wks,uu,vv,res_vc)
map2 = gsn_csm_contour_map(wks,elev,res)
overlay(map,plot_wind)
overlay(plot_wind,map)


;######   ADD MAP OF CHINA   #############

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

  plot = add_china_map(wks,map2,cnres)

draw(map2)
draw(map)
frame(wks)
      
end do
end

pkn789

只能发在下面了

pkn789

或者求山西的shapefile文件也可以，请问有吗？

pkn789

我找到山西的shapefile文件了，谢谢大家