画中国东北数据为什么会成这样

983019071

为什么我提取中国数据会成这样   哪个高手能告诉我以下这个程序哪里错了

;*************************************************
; shapefiles_9.ncl
;
; Concepts illustrated:
;   - Masking a data array based on a geographical area obtained from a shapefile
;   - Calculating an areal time series
;   - Make a time series plot
;
;*************************************************
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/csm/shea_util.ncl"
begin
;---Open shapefile and read lat/lon values.
;  dir     = ncargpath("data") + "/shp/"
  f       = addfile("bou2_4p.shp", "r")
  print(f)
  shp_lon = tofloat( f->x )
;  print(shp_lon)
  shp_lat = tofloat( f->y )
  printVarSummary(shp_lat)
  nshp    = dimsizes(shp_lon)
;---Shape file lon are -180 to 180                  
;   Make them 0-360 to match the GPCP grid
;  shp_lon = where(shp_lon.lt.0, shp_lon + 360, shp_lon)
;  shp_lat = where(shp_lat.lt.0, shp_lon + 180, shp_lat)
;---Get Max & Min lat/lon for the shape file
  min_shp_lat = min(shp_lat)
  max_shp_lat = max(shp_lat)
  min_shp_lon = min(shp_lon)
  max_shp_lon = max(shp_lon)
;---Open GPCP data file
;   Use NCL syntax to extract the smallest are that includes the shape file
  f      = addfile("cru_ts3.23.1981.1990.tmp.dat.nc","r")
  f1 = addfile("cru_ts3.23.1991.2000.tmp.dat.nc","r")
  f2 = addfile("cru_ts3.23.2001.2010.tmp.dat.nc","r")
  f3 = addfile("cru_ts3.23.2011.2014.tmp.dat.nc","r")
  
  u1     = f->tmp(:,{min_shp_lat:max_shp_lat},{min_shp_lon:max_shp_lon})
  u2 = f1->tmp(:,{min_shp_lat:max_shp_lat},{min_shp_lon:max_shp_lon})
u3 = f2->tmp(:,{min_shp_lat:max_shp_lat},{min_shp_lon:max_shp_lon})
  u4 = f3->tmp(:,{min_shp_lat:max_shp_lat},{min_shp_lon:max_shp_lon})
  
;    u1     = f->tmp(:,{38:53},{113:134})
;  u2 = f1->tmp(:,{38:53},{113:134})
; u3 = f2->tmp(:,{38:53},{113:134})
;  u4 = f3->tmp(:,{38:53},{113:134})
  printVarSummary(u1)
  dimp    = dimsizes(u1)
  printVarSummary(dimp)
  ntim    = dimp(0)
  nlat    = dimp(1)
  print(nlat)
  mlon    = dimp(2)
;---Create an array and initialize to _FillValue
  pmask1   = new(dimsizes(u1), typeof(u1), u1@_FillValue)
  pmask2   = new(dimsizes(u2), typeof(u2), u2@_FillValue)
  pmask3   = new(dimsizes(u3), typeof(u3), u3@_FillValue)
  pmask4   = new(dimsizes(u4), typeof(u4), u4@_FillValue)  
  copy_VarCoords(u1,pmask1)
  copy_VarCoords(u2,pmask2)
  copy_VarCoords(u3,pmask3)
  copy_VarCoords(u4,pmask4)
;---Keep only data within the polygon
;   Use NCL array syntax (:) to propagate to all times
  do nl=0,nlat-1
    do ml=0,mlon-1
      if(gc_inout(u1&lat(nl),u1&lon(ml),shp_lat,shp_lon)) then
         pmask1(:,nl,ml) = u1(:,nl,ml)
      end if
    end do
  end do
  delete(nl)
  delete(ml)
  
  do nl=0,nlat-1
    do ml=0,mlon-1
      if(gc_inout(u2&lat(nl),u2&lon(ml),shp_lat,shp_lon)) then
         pmask2(:,nl,ml) = u2(:,nl,ml)
      end if
    end do
  end do
  delete(nl)
  delete(ml)
  
   do nl=0,nlat-1
    do ml=0,mlon-1
      if(gc_inout(u3&lat(nl),u3&lon(ml),shp_lat,shp_lon)) then
         pmask3(:,nl,ml) = u3(:,nl,ml)
      end if
    end do
  end do
  delete(nl)
  delete(ml)
  
   do nl=0,nlat-1
    do ml=0,mlon-1
      if(gc_inout(u4&lat(nl),u4&lon(ml),shp_lat,shp_lon)) then
         pmask4(:,nl,ml) = u4(:,nl,ml)
      end if
    end do
  end do
  delete(nl)
  delete(ml)
printVarSummary(pmask1)
ave1 = new((/10,30,42/),double)
ave2 = new((/10,30,42/),double)
ave3 = new((/10,30,42/),double)
ave4 = new((/10,30,42/),double)
tmp = new((/34,30,42/),double)
do i = 0,9
   ave1(i,:,:)=dim_avg_n_Wrap(pmask1(i*12+4:i*12+8,{38:53},{113:134}),0)
   ave2(i,:,:)=dim_avg_n_Wrap(pmask2(i*12+4:i*12+8,{38:53},{113:134}),0)
   ave3(i,:,:)=dim_avg_n_Wrap(pmask2(i*12+4:i*12+8,{38:53},{113:134}),0)
end do
delete(i)
do i=0,3
   ave4(i,:,:)=dim_avg_n_Wrap(pmask4(i*12+4:i*12+8,{38:53},{113:134}),0)
end do
tmp(0:9,:,:)= ave1(0:9,:,:)
tmp(10:19,:,:)= ave2(0:9,:,:)
tmp(20:29,:,:)= ave3(0:9,:,:)
tmp(30:33,:,:)= ave4(0:3,:,:)
tmp!0="time"
tmp&time = ispan(1981,2014,1)
printVarSummary(tmp)
;      EOF
  neof   = 3        ; number of EOFs
  optEOF = True      
  optEOF@jopt = 0   ; This is the default; most commonly used; no need to specify.
;;optEOF@jopt = 1   ; **only** if the correlation EOF is desired
  optETS = False
  latS   =  38.
  latN   =  54.
  lonL   = 113.
  lonR   =  134.
  yrStrt = 1981
  yrLast = 2014
  
  wSLP = tmp
  ;wSLP   = tmp*conform(tmp, clat(256:289), 1)
  wSLP@long_name = "tmp: "+wSLP@long_name
  x      = wSLP({lat|latS:latN},{lon|lonL:lonR},time|:)
  eof    = eofunc_Wrap(x, neof, optEOF)      
  eof_ts = eofunc_ts_Wrap (x, eof, optETS)
  printVarSummary( eof )                         ; examine EOF variables
  printVarSummary( eof_ts )
  
  ; Normalize time series: Sum spatial weights over the area of used
; =================================================================
;dimx   = dimsizes( x )
;mln    = dimx(1)
;sumWgt = mln*sum( clat({lat|latS:latN}) )
;eof_ts = eof_ts/sumWgt
;  yyyymm = cd_calendar(eof_ts&time,-1)/100
  
  time = fspan(1981,2014,34)
  wks = gsn_open_wks("pdf","EOF5-9mon")
  gsn_define_colormap(wks,"BlWhRe")       ; choose colormap
  plot = new(neof,graphic)                ; create graphic array
                                          ; only needed if paneling
; EOF patterns
  res                      = True         
  res@gsnDraw              = False        ; don't draw yet
  res@gsnFrame             = False        ; don't advance frame yet
;---This resource not needed in V6.1.0
  res@gsnSpreadColors      = True         ; spread out color table
  res@gsnAddCyclic         = False        ; plotted dataa are not cyclic

  res@mpFillOn             = False        ; turn off map fill
  res@mpMinLatF            = latS         ; zoom in on map
  res@mpMaxLatF            = latN
  res@mpMinLonF            = lonL
  res@mpMaxLonF            = lonR
  res@mpDataBaseVersion = "MediumRes"  
res@mpLimitMode           = "LatLon"
res@mpLandFillColor  = "white"
res@mpDataSetName="Earth..4" ;huizhizhongguoditu
res@mpOutlineBoundarySets  = "Geophysical"
res@mpOutlineSpecifiers=(/"China","jilin","liaoning","heilongjiang","Nei Mongol"/)
res@mpProvincialLineColor="black"   
res@mpProvincialLineThicknessF =1
  res@cnFillOn             = True         ; turn on color fill
  res@cnLinesOn            = False        ; True is default
;res@cnLineLabelsOn       = False        ; True is default
  res@lbLabelBarOn         = False        ; turn off individual lb's
                                          ; set symmetric plot min/max
  symMinMaxPlt(eof, 16, False, res)       ; contributed.ncl
; panel plot only resources
  resP                     = True         ; modify the panel plot
  resP@gsnMaximize         = True         ; large format
  resP@gsnPanelLabelBar    = True         ; add common colorbar
  resP@lbLabelAutoStride   = True         ; auto stride on labels
  resP@txString             = "tmp: "+yrStrt+"-"+yrLast+".may-ste"
   
   do n=0,neof-1
     res@gsnLeftString  = "EOF "+(n+1)
     res@gsnRightString = sprintf("%5.1f", eof@pcvar(n)) +"%"
     plot(n)=gsn_csm_contour_map_ce(wks,eof(n,:,:),res)
  end do
  gsn_panel(wks,plot,(/neof,1/),resP)     ; now draw as one plot
; second plot
;*******************************************
; EOF time series  [bar form]
  rts           = True
  rts@gsnDraw   = False       ; don't draw yet
  rts@gsnFrame  = False       ; don't advance frame yet
  rts@gsnScale  = True        ; force text scaling               
; these four rtsources allow the user to stretch the plot size, and
; decide exactly where on the page to draw it.
  rts@vpHeightF = 0.40        ; Changes the aspect ratio
  rts@vpWidthF  = 0.85
  rts@vpXF      = 0.10        ; change start locations
  rts@vpYF      = 0.75        ; the plot

  rts@tiYAxisString = "Centigrade"                    ; y-axis label      
  rts@gsnYRefLine           = 0.              ; reference line   
  rts@gsnXYBarChart         = True            ; create bar chart
  rts@gsnAboveYRefLineColor = "red"           ; above ref line fill red
  rts@gsnBelowYRefLineColor = "blue"          ; below ref line fill blue
; panel plot only resources
  rtsP                      = True            ; modify the panel plot
  rtsP@gsnMaximize          = True            ; large format
  rtsP@txString             = "tmp: "+yrStrt+"-"+yrLast+".may-ste"
; create individual plots
  do n=0,neof-1
     rts@gsnLeftString  = "EOF "+(n+1)
     rts@gsnRightString = sprintf("%5.1f", eof@pcvar(n)) +"%"
     plot(n) = gsn_csm_xy (wks,time,eof_ts(n,:),rts)
  end do
  gsn_panel(wks,plot,(/neof,1/),rtsP)     ; now draw as one plot
end

freekiller

提问题最好是给出错误，而不是贴这么长的脚步。这样没谁回去帮你查

983019071

我不知道哪里错了，所以无法给出错误在哪里？