如何只做出shapefile数据中对应区域的ncl图

FrankieLJY

尽头的尽头 发表于 2017-3-7 22:18
复杂一些的要写代码，做底图当然不用啦

还是写了好多代码 不过是按照官网上的例子改的 按网上的例子 导出风没有问题 但是导出我的温度这个还是有问题，不知道版主能不能帮我看看代码

尽头的尽头

FrankieLJY 发表于 2017-3-8 15:11
还是写了好多代码 不过是按照官网上的例子改的 按网上的例子 导出风没有问题 但是导出我的温度这个还是有 ...

你目的是干嘛呢，你发来我看看

FrankieLJY

尽头的尽头 发表于 2017-3-8 16:15
你目的是干嘛呢，你发来我看看

版主不好意思，这两天没有上家园
我是想两个时间做京津冀地区的平均温差值，并且通过mask只保留京津冀地区部分
下面代码是我做的两个月的平均温度的差值，虽然能画出图，但是不知道做对了没有，请批评指正
begin
    dir1="/data1/Jingjinji2000-201001/"                      ;;;;;;;;;;;a1 前一个年份
    files1= "wrfout2_2000_01"
    a1=addfiles(dir1 + files1 + ".nc", "r")

    dir2="/data1/Jingjinji2010-201001/"           ;;;;;;;;;;;a2 后一个年份
    files2= "wrfout2_2010_01"
    a2=addfiles(dir2 + files2 + ".nc", "r")

    b=addfile("./wrfout_d02_2010-01-15_00:00:00.nc", "r")；我的文件b只是为了给最后wrf_map_overlays赋值

    type="png"
    wks=gsn_open_wks(type, "ave_diff_geo2005-2010_201001")

    res=True
    res@cnFillMode="RasterFill"

    pltres=True
    mpres=True
    ;mpres@mpFillOn = True
    mpres@mpDataBaseVersion           = "Ncarg4_1"
    mpres@mpDataSetName               = "Earth..4"
    mpres@mpNationalLineColor         = "Black"
    mpres@mpUSStateLineColor          = "Black"
    mpres@mpOutlineBoundarySets       = "AllBoundaries"
   
    xhr1 =wrf_user_getvar(a1, "Q2", -1)
    tc1=dim_avg_n(xhr1, 0)
    tc1=tc1-273.16
    tc1@description="Surface Temperature"
    tc1@units="C"

    xhr2 =wrf_user_getvar(a2, "Q2", -1)
    tc2=dim_avg_n(xhr2, 0)
    tc2=tc2-273.16
    tc2@description="Surface Temperature"
    tc2@units="C"

    tc=tc2-tc1
    ;tc=tc*2
    tc@description="Surface Temperature"
    tc@units="C"

    opts=res
    opts@cnFillOn=True
    opts@ContourParameters=(/-1.,1.,0.02/)
    opts@gsnSpreadColorEnd=-10
    contour_tc=wrf_contour(b, wks, tc, opts)        ;添加的是多个文件，b该如何添加
    delete(opts)

    plot=wrf_map_overlays(b, wks, (/contour_tc/), pltres, mpres)
end

FrankieLJY

尽头的尽头 发表于 2017-3-8 16:15
你目的是干嘛呢，你发来我看看

上一段代码的图是我在一楼贴出来的图，但是京津冀以外的地区都有值，所以我想加一个mask，把兴趣区域以外的部分都给清除掉，所以我在NCL官网上找到一个案例http://www.ncl.ucar.edu/Applications/Scripts/shapefiles_5.ncl，进行修改，成如下脚本，但是一直有错误，希望版主有时间的时候帮我看看好吗，谢谢！
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/wrf/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

undef("calcPolygonAreaAndCentroid")
function calcPolygonAreaAndCentroid(x, y, firstVert, lastVert)
  local cY, cY, a, tmp
begin
  cX = 0.d
  cY = 0.d
  area = 0.d
  do i=firstVert,lastVert
    if (i .eq. lastVert) then
      j = firstVert
    else
      j = i + 1
    end if
    tmp = x(i)*y(j) - x(j)*y(i)
    area = area + tmp
    cX = cX + (x(i) + x(j))*tmp
    cY = cY + (y(i) + y(j))*tmp
  end do

  area = area / 2.0
  cX = cX / (6.0 * area)
  cY = cY / (6.0 * area)
  ; recall that the area calculation may yield a negative result,
  ; depending upon CW vs. CCW ordering of the vertices.
  return (/ abs(area), cX, cY /)
end

begin
    dir1="/data1/Jingjinji2000-201001/"                      ;;;;;;;;;;;a1 前一个年份
    ;files1=systemfunc("ls wrfout_d02_2010-01-*")
    files1= "wrfout2_2000_01"
    a1=addfiles(dir1 + files1 + ".nc", "r")

    dir2="/data1/Jingjinji2010-201001/"           ;;;;;;;;;;;a2 后一个年份
    ;files2=systemfunc("ls wrfout_d02_2005-01-*")
    files2= "wrfout2_2010_01"
    a2=addfiles(dir2 + files2 + ".nc", "r")

    b=addfile("./wrfout_d02_2010-01-15_00:00:00.nc", "r")

    wks=gsn_open_wks("png", "a_test")
;======================================
;对各个参数进行计算
;======================================
    xhr1 =wrf_user_getvar(a1, "T2", -1)
    tc1=dim_avg_n(xhr1, 0)
    tc1=tc1-273.16
    tc1@description="Surface Temperature"
    tc1@units="C"

    xhr2 =wrf_user_getvar(a2, "T2", -1)
    tc2=dim_avg_n(xhr2, 0)
    tc2=tc2-273.16
    tc2@description="Surface Temperature"
    tc2@units="C"

    tc=tc2-tc1
    ;tc=tc*2
    tc@description="Surface Temperature"
    tc@units="C"

    tc@lat2d = b->XLAT(0,:,:)         ; direct assignment
    tc@lon2d = b->XLONG(0,:,:)

    sCritical    = 6.0
    tc@_FillValue = 255

  ; ----------------------------------------------
  ; Get the national boundary from a shapefile.
  ; We'll use this as a mask for the data.
  ; ----------------------------------------------
  natBdry = addfile("./JJJ.shp", "r")
  
  maskedS = new(dimsizes(b),typeof(b),b@_FillValue)
  maskedS@lat2d = b@lat2d
  maskedS@lon2d = b@lon2d

  ; Get the size of the lat/lon grid of variable "s"
  sDims = dimsizes(b)
  iNumLat = sDims(1)
  iNumLon = sDims(2)

  ; Put data in the areas that we don't want masked.
  do j=0,iNumLat-1
    print("masking row " + j)
    do i=0,iNumLon-1
      if(gc_inout(b->XLAT(0,j,i), b->XLONG(0,j,i), natBdry->y, natBdry->x)) then
        maskedS(:,j,i) = tc(:,j,i)
      end if
    end do
  end do

; We want a very specific colormap...
  cmap = (/ (/1., 1., 1./), (/ 0., 0., 0./),          \  ; background/foreground
           (/ 231./255.,  250./255.,  254./255. /),  \  ; color for ocean
           (/ 230./255.,  236./255.,  236./255. /),  \  ; color for land outside AOI
            (/ .3,         .3,         .3        /),  \  ; color for city labels
            (/ 1.,         1.,         1.        /),  \  ; colors for data mapping...
            (/ 221./255.,  184./255.,  129./255. /),  \
            (/ 254./255.,  157./255.,  014./255. /),  \
            (/ 253./255.,  075./255.,  251./255. /),  \
            (/ 116./255.,  042./255.,  253./255. /),  \
            (/ 255./255.,  001./255.,  001./255. /),  \
            (/ 034./255.,  062./255.,  255./255. /)   \
         /)
  gsn_define_colormap(wks, cmap)

  res                       = True             ; plot mods desired
  res@gsnFrame              = False            
  res@gsnDraw               = False
  res@gsnMaximize           = True
  res@cnFillColors          = cmap(5:,:)
  res@cnLevelSelectionMode  = "ExplicitLevels"
  res@cnLevels              = (/ 5.4, 6.2, 6.9, 7.4, 7.8, 8.6 /)
  res@cnMinLevelValF        =  5.0
  res@cnMaxLevelValF        =  10.0
  res@cnFillOn              = True             ; color plot desired
  res@cnLinesOn             = False            ; turn off contour lines
  res@cnLineLabelsOn        = False            ; turn off contour labels

  ; -----------------------------------------------
  ; Set map extent to zoom in on Pakistan
  ; -----------------------------------------------
  res@mpMinLonF = 113.
  res@mpMaxLonF = 120.
  res@mpMinLatF = 35.
  res@mpMaxLatF = 43.
  res@mpGeophysicalLineColor = 1
  res@mpOceanFillColor       = 2                ; color for sea
  res@mpLandFillColor        = 3                ; color for land outside of borders
  res@mpOutlineBoundarySets  = "Geophysical"    ; draw coastlines for entire region
  res@mpGridAndLimbOn        = True
  res@mpGridLineColor        = (/ .25, .25, .25 /)
  res@mpGridLatSpacingF      = 4.0
  res@mpGridLonSpacingF      = 4.0

  ; ----------------------------------------------
  ; Turn on lat / lon labeling
  ; ----------------------------------------------
  res@pmTickMarkDisplayMode = "Always"         ; turn on tickmarks
  res@tmXTOn = False            ; turn off top   labels
  res@tmYROn = False            ; turn off right labels
   
  ; -----------------------------------------------
  ; Loop over all times and levels ( uncomment )
  ; Demo: one arbitrarily closen time and level  
  ; -----------------------------------------------
  dims  = dimsizes(tc)                          ; dimensions of x
  ntim  = dims(0)                              ; number of time steps
  klev  = dims(1)                              ; number of "bottom_top" levels

  nt    =  0
;;do nt=0,ntim-1                               ; uncomment for loop
;;  do ll=0,klev-1
       res@tiMainString  = Times(nt)
       res@lbLabelBarOn  = False
       plot    = gsn_csm_contour_map(wks,maskedS(nt,:,:),res)

draw(plot)
frame(wks)
end

FrankieLJY

尽头的尽头 发表于 2017-3-8 16:15
你目的是干嘛呢，你发来我看看

上一段代码的图是我在一楼贴出来的图，但是京津冀以外的地区都有值，所以我想加一个mask，把兴趣区域以外的部分都给清除掉，所以我在NCL官网上找到一个案例http://www.ncl.ucar.edu/Applications/Scripts/shapefiles_5.ncl，进行修改，成如下脚本，但是一直有错误，希望版主有时间的时候帮我看看好吗，谢谢！
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/wrf/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

undef("calcPolygonAreaAndCentroid")
function calcPolygonAreaAndCentroid(x, y, firstVert, lastVert)
  local cY, cY, a, tmp
begin
  cX = 0.d
  cY = 0.d
  area = 0.d
  do i=firstVert,lastVert
    if (i .eq. lastVert) then
      j = firstVert
    else
      j = i + 1
    end if
    tmp = x(i)*y(j) - x(j)*y(i)
    area = area + tmp
    cX = cX + (x(i) + x(j))*tmp
    cY = cY + (y(i) + y(j))*tmp
  end do

  area = area / 2.0
  cX = cX / (6.0 * area)
  cY = cY / (6.0 * area)
  ; recall that the area calculation may yield a negative result,
  ; depending upon CW vs. CCW ordering of the vertices.
  return (/ abs(area), cX, cY /)
end

begin
    dir1="/data1/Jingjinji2000-201001/"                      ;;;;;;;;;;;a1 前一个年份
    ;files1=systemfunc("ls wrfout_d02_2010-01-*")
    files1= "wrfout2_2000_01"
    a1=addfiles(dir1 + files1 + ".nc", "r")

    dir2="/data1/Jingjinji2010-201001/"           ;;;;;;;;;;;a2 后一个年份
    ;files2=systemfunc("ls wrfout_d02_2005-01-*")
    files2= "wrfout2_2010_01"
    a2=addfiles(dir2 + files2 + ".nc", "r")

    b=addfile("./wrfout_d02_2010-01-15_00:00:00.nc", "r")

    wks=gsn_open_wks("png", "a_test")
;======================================
;对各个参数进行计算
;======================================
    xhr1 =wrf_user_getvar(a1, "T2", -1)
    tc1=dim_avg_n(xhr1, 0)
    tc1=tc1-273.16
    tc1@description="Surface Temperature"
    tc1@units="C"

    xhr2 =wrf_user_getvar(a2, "T2", -1)
    tc2=dim_avg_n(xhr2, 0)
    tc2=tc2-273.16
    tc2@description="Surface Temperature"
    tc2@units="C"

    tc=tc2-tc1
    ;tc=tc*2
    tc@description="Surface Temperature"
    tc@units="C"

    tc@lat2d = b->XLAT(0,:,:)         ; direct assignment
    tc@lon2d = b->XLONG(0,:,:)

    sCritical    = 6.0
    tc@_FillValue = 255

  ; ----------------------------------------------
  ; Get the national boundary from a shapefile.
  ; We'll use this as a mask for the data.
  ; ----------------------------------------------
  natBdry = addfile("./JJJ.shp", "r")
  
  maskedS = new(dimsizes(b),typeof(b),b@_FillValue)
  maskedS@lat2d = b@lat2d
  maskedS@lon2d = b@lon2d

  ; Get the size of the lat/lon grid of variable "s"
  sDims = dimsizes(b)
  iNumLat = sDims(1)
  iNumLon = sDims(2)

  ; Put data in the areas that we don't want masked.
  do j=0,iNumLat-1
    print("masking row " + j)
    do i=0,iNumLon-1
      if(gc_inout(b->XLAT(0,j,i), b->XLONG(0,j,i), natBdry->y, natBdry->x)) then
        maskedS(:,j,i) = tc(:,j,i)
      end if
    end do
  end do

; We want a very specific colormap...
  cmap = (/ (/1., 1., 1./), (/ 0., 0., 0./),          \  ; background/foreground
           (/ 231./255.,  250./255.,  254./255. /),  \  ; color for ocean
           (/ 230./255.,  236./255.,  236./255. /),  \  ; color for land outside AOI
            (/ .3,         .3,         .3        /),  \  ; color for city labels
            (/ 1.,         1.,         1.        /),  \  ; colors for data mapping...
            (/ 221./255.,  184./255.,  129./255. /),  \
            (/ 254./255.,  157./255.,  014./255. /),  \
            (/ 253./255.,  075./255.,  251./255. /),  \
            (/ 116./255.,  042./255.,  253./255. /),  \
            (/ 255./255.,  001./255.,  001./255. /),  \
            (/ 034./255.,  062./255.,  255./255. /)   \
         /)
  gsn_define_colormap(wks, cmap)

  res                       = True             ; plot mods desired
  res@gsnFrame              = False            
  res@gsnDraw               = False
  res@gsnMaximize           = True
  res@cnFillColors          = cmap(5:,:)
  res@cnLevelSelectionMode  = "ExplicitLevels"
  res@cnLevels              = (/ 5.4, 6.2, 6.9, 7.4, 7.8, 8.6 /)
  res@cnMinLevelValF        =  5.0
  res@cnMaxLevelValF        =  10.0
  res@cnFillOn              = True             ; color plot desired
  res@cnLinesOn             = False            ; turn off contour lines
  res@cnLineLabelsOn        = False            ; turn off contour labels

  ; -----------------------------------------------
  ; Set map extent to zoom in on Pakistan
  ; -----------------------------------------------
  res@mpMinLonF = 113.
  res@mpMaxLonF = 120.
  res@mpMinLatF = 35.
  res@mpMaxLatF = 43.
  res@mpGeophysicalLineColor = 1
  res@mpOceanFillColor       = 2                ; color for sea
  res@mpLandFillColor        = 3                ; color for land outside of borders
  res@mpOutlineBoundarySets  = "Geophysical"    ; draw coastlines for entire region
  res@mpGridAndLimbOn        = True
  res@mpGridLineColor        = (/ .25, .25, .25 /)
  res@mpGridLatSpacingF      = 4.0
  res@mpGridLonSpacingF      = 4.0

  ; ----------------------------------------------
  ; Turn on lat / lon labeling
  ; ----------------------------------------------
  res@pmTickMarkDisplayMode = "Always"         ; turn on tickmarks
  res@tmXTOn = False            ; turn off top   labels
  res@tmYROn = False            ; turn off right labels
   
  ; -----------------------------------------------
  ; Loop over all times and levels ( uncomment )
  ; Demo: one arbitrarily closen time and level  
  ; -----------------------------------------------
  dims  = dimsizes(tc)                          ; dimensions of x
  ntim  = dims(0)                              ; number of time steps
  klev  = dims(1)                              ; number of "bottom_top" levels

  nt    =  0
;;do nt=0,ntim-1                               ; uncomment for loop
;;  do ll=0,klev-1
       res@tiMainString  = Times(nt)
       res@lbLabelBarOn  = False
       plot    = gsn_csm_contour_map(wks,maskedS(nt,:,:),res)

draw(plot)
frame(wks)
end

尽头的尽头

参考这里http://bbs.06climate.com/forum.p ... &extra=page%3D1

FrankieLJY

尽头的尽头 发表于 2017-3-10 21:05
参考这里http://bbs.06climate.com/forum.php?mod=viewthread&tid=45252&extra=page%3D1

谢谢版主 我去学习一下

南宫蓝天

pigzero527 发表于 2017-3-3 09:52
可以用shapefile_mask_data函数对数据做mask，会把shape外区域都变成缺测值，画图时就会只显示区域内的了， ...

在NCL官网www.ncl.ucar.edu/Document/Functions/list_alpha.shtml里面，并没有shapefile_mask_data这个函数存在啊？在哪里找的？

pigzero527

南宫蓝天 发表于 2017-7-28 17:08
在NCL官网www.ncl.ucar.edu/Document/Functions/list_alpha.shtml里面，并没有shapefile_mask_data这个函 ...

这个函数是另外有人写好脚本加进去的，有个shapefile_mask_data.ncl加到lib里就可以用了，不过具体在哪儿找到的我已经忘了

irides

请问楼主京津冀地区的shp文件在哪里能下载啊？