ncl在mask青藏高原时出错

嗳づ左岸の兔

与http://bbs.06climate.com/forum.php?mod=viewthread&tid=30626帖子相似
我也是按照官网修改的，但是总是出现warning:ContourPlotInitialize: no valid values in scalar field; ContourPlot not possible:[errno=1101]
不知道为什么赋值不成功，求助大神~
begin

a=addfile("20080501000000.nc","r")

sm=a->sm
data1=sm(0,:,:)

  minlat      =   89.875
  maxlat      =  -89.875
  minlon      = -179.875
  maxlon      = 179.875
  nlat=720
  nlon=1440
  lat1d       = fspan(minlat,maxlat,nlat)
  lon1d       = fspan(minlon,maxlon,nlon)

  lat1d@units = "degrees_north"
  lon1d@units = "degrees_east"

  MASK_INSIDE = False    ; Whether to mask data inside or outside the
                         ; given geographical area.

;---Attach lat/lon coordinate array information.
  data=new((/720,1440/),float)
  data=data1(:,:)
  data!0      = "lat"
  data!1      = "lon"
  data&lat    = lat1d
  data&lon    = lon1d
  printVarSummary(data)

;---Open shapefile and read Mississippi River Basin lat/lon values.
  f       = addfile("Tibetan_Output.shp", "r")
  mrb_lon = f->x
  mrb_lat = f->y
  nmrb    = dimsizes(mrb_lon)      ; dimension size of 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)

;mask
if(MASK_INSIDE) then
---Start with data filled in.
    data_mask = data
else
;---Start with data all missing
    data_mask = new(dimsizes(data),typeof(data),data@_FillValue)
    copy_VarCoords(data,data_mask)
  ;end if

  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
;---Put missing values in the areas that we want masked.
    do ilt=ilt1,ilt2
      do iln=iln1,iln2
        if(gc_inout(lat1d(ilt),lon1d(iln),mrb_lat,mrb_lon)) then
         data_mask(ilt,iln) = data_mask@_FillValue
        end if
      end do
    end do
  else
;---Put data back in the areas that we don't want masked.
    do ilt=ilt1,ilt2
      do iln=iln1,iln2
        if(gc_inout(lat1d(ilt),lon1d(iln),mrb_lat,mrb_lon)) then
          data_mask(ilt,iln) = data(ilt,iln)
        end if
      end do
    end do
  end if

;---Start the graphics
  wks  = gsn_open_wks("png","mask")        ; send graphics to PNG file

res=True

res@gsnDraw             = False          ; don't draw plot yet
  res@gsnFrame            = False          ; don't advance frame yet
map_data = gsn_csm_contour(wks,data,res)
map_mask      = gsn_csm_contour(wks,data_mask,res)
lnres                  = True
  lnres@gsLineColor      = "blue"
  lnres@gsLineThicknessF = 2.0            ; 2x thickness
  mrb_line_mask = gsn_add_polyline(wks, map_mask, mrb_lon, mrb_lat, lnres)
draw(map_mask)
draw(map_data)
  frame(wks)
end

嗳づ左岸の兔

已找到错误
在 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
由于经纬度不同，必要时要手动赋值

daxiongxiaojie

请问高原地形文件在哪里找？

嗳づ左岸の兔

参考http://bbs.06climate.com/forum.p ... p;extra=&page=1系列帖子

当时年少春衫薄

楼主你好，请问你画出来的是青藏高原轮廓还是阴影呢

嗳づ左岸の兔

mask掉其他区域的

xiaocaoqiqiao

嗳づ左岸の兔 发表于 2016-4-1 08:56
已找到错误
在 ilt1   = ilt_mn(dimsizes(ilt_mn)-1)    ; Start of lat box
  iln1   = iln_mn(dimsizes ...

怎么手动赋值的？能麻烦楼主细说一下吗？

嗳づ左岸の兔

xiaocaoqiqiao 发表于 2016-7-8 15:25
怎么手动赋值的？能麻烦楼主细说一下吗？

我当时只是改了个循环，没有手动赋值
手动赋值就是将经纬度一个个输入

xiaocaoqiqiao

嗳づ左岸の兔 发表于 2016-7-12 16:48
我当时只是改了个循环，没有手动赋值
手动赋值就是将经纬度一个个输入

谢谢你的回复  我后来用别的方法解决了

当时年少春衫薄

楼主您好请问 minlat      =   89.875
  maxlat      =  -89.875是否可以解释一下