NCLmasking青藏高原时经纬度反了

当时年少春衫薄

画青藏高原感热通量等值线图的时候出现了同经度热通量相反的情况，请问应该怎么办～
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
  MASK_INSIDE = False    ; Whether to mask data inside or outside the
                         ; given geographical area.

;---Generate some dummy data over map

  minlat      =   40
  maxlat      =  25
  minlon      = 73
  maxlon      =  105

  nlat            = 28
nlon=64
diri1= "/home/fanweiwei/NCEP/heat/ECMWF/nianretl/feirun/"
fils1 = systemfunc("ls " + diri1 + "sl*")
f1= addfiles (fils1  , "r")
ListSetType (f1, "join")   

x12= f1[:]->SSHF_GDS0_SFC(:,59:151,:,:)
y12=dim_avg_n(x12,(/0,1/))
l1=-y12/86400
diri2 = "/home/fanweiwei/NCEP/heat/ECMWF/nianretl/run/"
fils2= systemfunc("ls " + diri2 + "sl*")
f2= addfiles (fils2  , "r")
ListSetType (f2, "join")   
x22= f2[:]->SSHF_GDS0_SFC(:,60:152,:,:)
y22=dim_avg_n(x22,(/0,1/))
l2=-y22/86400
data=(l1*24+l2*8)/32

;---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

;---Open shapefile and read Mississippi River Basin lat/lon values.
  f       = addfile("/home/fanweiwei/NCEP/heat/TP/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
bvbb
      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","ss")        ; send graphics to PNG file

  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@gsnAddCyclic        = False          ; Don't add a cyclic point.

  res@mpDataBaseVersion   = "MediumRes"    ; slightly better resolution

;---Zoom in on North America.
  res@mpMinLatF           = minlat
  res@mpMaxLatF           = maxlat
  res@mpMinLonF           = minlon
  res@mpMaxLonF           = maxlon

res@cnFillOn=False
res@cnLinesOn=True
res@cnMonoLineColor=True


res@cnLineLabelsOn=True
res@cnLevels=(/10,15,20,25,30,35,40,45,50,55/)
res@cnLineLabelFontAspectF=2.0
res@cnLineLabelFontHeightF=0.008
res@cnMonoLevelFlag=True
res@cnLevelFlag="LineAndLabel"
res@cnLineLabelPlacementMode="Constant"
res@cnLineDashSegLenF=0.05

  res@tiMainString        = "F"
res@cnLevelSpacingF =10
;---Create contours over map.
  map_data = gsn_csm_contour_map(wks,data,res)

;---Resources for polyline
  lnres                  = True
  lnres@gsLineColor      = "black"
  lnres@gsLineThicknessF = 4.0            ; 2x thickness

;---Attach MRB outline to map.
  mrb_line = gsn_add_polyline(wks, map_data, mrb_lon, mrb_lat, lnres)

;---Draw plot and advance frame
  draw(map_data)
  frame(wks)

;---Draw contours of masked data
  res@tiMainString = "TP "

  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)
  frame(wks)

;---Draw the part of the lat/lon grid that was masked out.
;   Only draw the map this time (no contours).

  delete(res@gsnAddCyclic)

  res@tiMainString = "Area where lat/lon values were masked"
  map = gsn_csm_map(wks,res)


  draw(map)
  frame(wks)
end

当时年少春衫薄

talkd

hi您好 请问哪里能下到高原边界的shape file呢？谢谢

当时年少春衫薄

talkd 发表于 2016-11-21 10:43
hi您好 请问哪里能下到高原边界的shape file呢？谢谢

您好这是我下载的文件

talkd

当时年少春衫薄 发表于 2016-11-21 21:47
您好这是我下载的文件

谢谢!没看懂你的问题？哪里反了？

当时年少春衫薄

talkd 发表于 2016-11-22 08:33
谢谢!没看懂你的问题？哪里反了？

   数据反了，好像lat1d       = fspan(minlat,maxlat,nlat)改成   lat1d       = fspan(maxlat，minlat,nlat)就好了

123出去玩

楼主你脚本里面的bvbb是什么意思

当时年少春衫薄

123出去玩 发表于 2016-12-7 13:03
楼主你脚本里面的bvbb是什么意思

我也不知道什么意思，模仿别人的脚本的，那部分除了把变量修改了其余的都没改动

123出去玩

当时年少春衫薄 发表于 2016-12-8 11:09
我也不知道什么意思，模仿别人的脚本的，那部分除了把变量修改了其余的都没改动

我参考了楼主的脚本改了改，可以出图，感谢。

当时年少春衫薄

123出去玩 发表于 2016-12-8 18:22
我参考了楼主的脚本改了改，可以出图，感谢。

嗯嗯不谢～