ncl只画青藏高原地区出错！

蔷薇花烬

ncl新手上路，最近在做毕业设计，需要画青藏高原地区的土壤湿度分布图，于是乎，参考了家园上相关的帖子
http://bbs.06climate.com/forum.php?mod=viewthread&tid=46333
http://bbs.06climate.com/forum.php?mod=viewthread&tid=30626
http://bbs.06clim ate.com/forum.php?mod=viewthread&tid=43465
按照这些类似的帖子进行了修改，但是程序一直提示：fatal:Dimension sizes of left hand side and right hand side of assignment do not match
fatal:["Execute.c":8578]:Execute: Error occurred at or near line 46 in file
提示表明是数据的维度出现的问题，读入的nc数据是全球的资料，当程序中的经纬度定义为nc文件中原来的经纬度时是可以出图的，但是是全球的并不是高原地区的，改变经纬度就出现上述的问题了，真的是刚刚开始接触没多久，很多地方不懂，周围也没有用ncl画图的人，所以，恳请各位大大不要嫌弃，给予解答，感激不尽！！！
程序如下：
;*************************************************************
;  load  NCL draw library
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
   ; Whether to mask data inside or outside the
                         ; given geographical area.
;---Generate some dummy data over map
  
  MASK_INSIDE = False
;******************************************************
minlat      =   25
maxlat      =   40
minlon      =   73
maxlon      =   105
nlat        = 69
nlon        = 121
filename = addfile("D:/LALALA/ncl-starter/taizhan/1948/GLDAS_NOAH025_M.A194801.020.nc4","r")
data= filename->SoilMoi0_10cm_inst(0,:,:)
;*************************************************************
data@_FillValue = -9999

;---Create dummy 1D lat/lon arrays
  lat1d       = fspan(minlat,maxlat,nlat)
  lon1d       = fspan(minlon,maxlon,nlon)
  lat1d@units = "degrees_north"
  lon1d@units = "degrees_east"
  
;---Attach lat/lon coordinate array information.
  
;---------------------------------------------------------
  
  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("D:/LALALA/ncl-starter/taizhan/DBATP/DBATP_Polygon.shp","r")
  mrb_lon = f->x
  mrb_lat = f->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
;---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
;
; Get the approximate index values that contain the
; Mississippi River Basin area. We'll use this to
; narrow in on the area that we need to figure out
; the masking for. Any lat/lon values outside this
; area we don't care about.
;
; This will make our gc_inout loop below go almost 4x
; faster, because we're not checking every single
; lat/lon point to see if it's within the area of
; interest.
;
  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
;---Put data back in the areas that we don't want masked.
  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_mask(ilt,iln) = data_mask@_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_mask(ilt,iln) = data(ilt,iln)
        end if
      end do
    end do
  end if

;---Start the graphics
  wks  = gsn_open_wks("X11","mask")        ; send graphics to PNG file
  gsn_define_colormap(wks, "BlueRed")
  
  res                     = True
  res@gsnMaximize         = True; maximize plot in frame
  res@gsnDraw             = False          ; don't draw plot yet
  res@gsnFrame            = False          ; don't advance frame yet
  
   res@mpMinLatF               = minlat
   res@mpMaxLatF               = maxlat
   res@mpMinLonF               = minlon
   res@mpMaxLonF               = maxlon
   res@gsnAddCyclic        = False          ; Don't add a cyclic point.
  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 = gsn_add_polyline(wks, map_data, mrb_lon, mrb_lat, lnres)
  draw(map_data)
  frame(wks)
  map_mask      = gsn_csm_contour_map(wks,data_mask,res)
  mrb_line_mask = gsn_add_polyline(wks, map_mask, mrb_lon, mrb_lat, lnres)
  draw(map_mask)
  
  
  draw(map_data)
  frame(wks)
end
nc数据的文件信息：
Variable: data
Type: float
Total Size: 3456000 bytes
            864000 values
Number of Dimensions: 2
Dimensions and sizes:   [lat | 600] x [lon | 1440]
Coordinates:
            lat: [-59.875..89.875]
            lon: [-179.875..179.875]
Number Of Attributes: 10
  time :           0
  units :       kg m-2
  standard_name :       soil_moisture_content
  long_name :   Soil moisture
  scale_factor :         1
  add_offset :   0
  missing_value :       -9999
  _FillValue :  -9999
  vmin :         2
  vmax :        47.58773

astiny

你看下min(mrb_lat)这个函数的维度吧。。。感觉应该是有俩维？

蔷薇花烬

astiny 发表于 2017-5-7 12:59
你看下min(mrb_lat)这个函数的维度吧。。。感觉应该是有俩维？

后面的是直接按照官网上的内容改的，感觉是data这个变量的问题，但是一直改不对，我把原来提示有错误的那两行data&lat    = lat1d，data&lon    = lon1d删除了之后，出来了青藏高原的轮廓，但是没有等值线的数值，提示:ContourPlotInitialize: no valid values in scalar field; ContourPlot not possible:[errno=1101]
看以前家园的帖子有人说是mask时候的循环的问题，但是没有找到如何解决

18408251130

也不知道我说的对不对啊，我发现GLDAS_NOAH025的月平均数据是没有土壤湿度的，最后是缺测吧？

Fy-Tail

楼主，解决了问题吗，出现了一样的问题

蔷薇花烬

Fy-Tail 发表于 2017-11-13 22:58
楼主，解决了问题吗，出现了一样的问题

我的是重复设置坐标的问题

duuning

66666666666666

qingshan

楼主：我跟你出现同样的问题，请问重复坐标怎么修改啊？
fatal:Dimension sizes of left hand side and right hand side of assignment do not match
fatal:["Execute.c":8637]:Execute: Error occurred at or near line 35 in file QTP_MASK_NEW.ncl

541631787

很有用{:5_213:}{:5_213:}

yeah...

请问你这个后来修改成功了吗