NCL mask青藏高原以外的区域出错，求大神指点

笨北极熊笨

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"

begin
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;---add 4 data and calcute the correlations
f=addfile(“mean.nc","r")
pmi=f->pmi
ts=pmi(lat|:,lon|:,time|:)
;-----------------------------------------------------
f1=addfile("mon6sm1.nc","r")
sm1=f1->swvl1
lon=f1->longitude
lat=f1->latitude
ts1=sm1(latitude|:,longitude|:,time|:)
e1=esccr(ts,ts1,0)
e1!2="lat"
e1!3="lon"
e1&lat=lat
e1&lon=lon
cor1=new((/33,57/),"float")
cor1=e1(0,0,:,:,0)
cor1@_FillValue =-999
print(max(cor1))
print(min(cor1))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
MASK_INSIDE = False
filename =”/tibetan.shp"  
  f5       = addfile(filename, "r")
  mrb_lon = f5->x
  mrb_lat = f5->y
  nmrb    = dimsizes(mrb_lon)

  min_mrb_lat = min(mrb_lat)
  max_mrb_lat = max(mrb_lat)
  min_mrb_lon = min(mrb_lon)
  max_mrb_lon = max(mrb_lon)
  ; Create new data to mask. Depending on whether you want to
; mask data inside or outside the geographical outline,
; start off with your data grid all missing, or filled in.
;
  if(MASK_INSIDE) then
    data_mask1 = cor1     
  else
    data_mask1 = new(dimsizes(cor1),typeof(cor1),cor1@_FillValue)
    copy_VarCoords(cor1,data_mask1)
     end if
  lat1d       = fspan(24,48,33)
  lon1d       = fspan(63,105,57)
  lat1d@units = "degrees_north"
  lon1d@units = "degrees_east"
  ilt_mn = ind(min_mrb_lat.gt.lat1d)
  ilt_mx = ind(max_mrb_lat.lt.lat1d)
  iln_mn = ind(min_mrb_lon.gt.lon1d)
  iln_mx = ind(max_mrb_lon.lt.lon1d)
  ilt1   = ilt_mn(dimsizes(ilt_mn)-1)    ; Start of lat box
  iln1   = iln_mn(dimsizes(iln_mn)-1)    ; Start of lon box
  ilt2   = ilt_mx(0)                     ; End of lat box
  iln2   = iln_mx(0)                     ; End of lon box

  if(MASK_INSIDE) then
    do ilt=ilt1,ilt2
      do iln=iln1,iln2
        if(gc_inout(lat1d(ilt),lon1d(iln),mrb_lat,mrb_lon)) then
          data_mask1(ilt,iln) = data_mask1@_FillValue
           end if
      end do
    end do
  else
    do ilt=ilt1,ilt2
      do iln=iln1,iln2
        if(gc_inout(lat1d(ilt),lon1d(iln),mrb_lat,mrb_lon)) then
          data_mask1(ilt,iln) = cor1(ilt,iln)
              end if
      end do
    end do
  end if
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;>=========================================================<
  wks       = gsn_open_wks("eps","1")     
res                         = True            
res@tiMainString            =""
res@gsnMaximize             =False
res@gsnDraw                 = False
res@gsnFrame                = False
res@mpMinLatF               = 24
res@mpMaxLatF               = 48
res@mpMinLonF               = 63
res@mpMaxLonF               = 105
res@mpLandFillColor         = "white"
res@cnFillDrawOrder         = "Draw"
;------------------------------------------------------                 
plot(0)=gsn_csm_contour_map(wks,data_mask1,res)
  ;-----------------------------------------------------------------------
      ;---Attach the polylines  
      pres             = True  
      pres@gsLineColor = "blue"  
      pres@gsLineThicknessF =2   
      poly1 = gsn_add_shapefile_polylines(wks,plot(0),filename,pres）     
      draw(plot)  
      frame(wks)  


    end  
本来是panel图，把其他重复的都删除，这个程序是按照官网上http://www.ncl.ucar.edu/Applications/Scripts/mask_9.ncl改的，但画出来就是这个效果，求大神指点。目的是想把青藏高原以外的等值线都mask掉，感觉错位了。数据范围是24-48N,63-105E，0.75*0.75的

longlivehj

感觉是数据纬度的方向反了

嗳づ左岸の兔

我也是照着例子改的~不知道楼主遇到这个问题吗warning:ContourPlotInitialize: no valid values in scalar field; ContourPlot not possible:[errno=1101]

limoxi

想请问楼主的青藏高原.shp文件是从哪里下载的？

cpp803520

嗳づ左岸の兔 发表于 2016-3-31 19:56
我也是照着例子改的~不知道楼主遇到这个问题吗warning:ContourPlotInitialize: no valid values in scalar  ...

你好，请问你的问题解决了么，我照着15的例子修改想画出中国区域的图也出现了和你一样的问题？请问你是怎么解决的？

嗳づ左岸の兔

cpp803520 发表于 2016-5-11 09:42
你好，请问你的问题解决了么，我照着15的例子修改想画出中国区域的图也出现了和你一样的问题？请问你是怎 ...

好像是经纬度的问题，与图中的mask循环有关

freekiller

数据的南北向反了，

cpp803520

嗳づ左岸の兔 发表于 2016-5-11 11:36
好像是经纬度的问题，与图中的mask循环有关

这样啊，你是怎么解决的呢？经纬度反了么？

当时年少春衫薄

楼主，你好，我按照你的脚本画出来的是青藏高原轮廓是实的，不是轮廓，请问问题出在哪里吗？

一日的天光

mark一下吧！