NCL画到某个区域之后有空缺的格点，怎么回事呢？

阿吉木柯 · 发表于 2014-7-14 18:51:40

登录后查看更多精彩内容~

您需要登录才可以下载或查看，没有帐号？立即注册

x

我现在正画4个季节的水汽图。
下面这个是春季的没有画到某个区域的图，可以看到没有缺测，并且我导出数据，也没发现缺测值。

而我把这个数据弄到某个区域之后，竟然会有貌似两个格点的空白值，这个区域是我在arcgis中弄的.shp文件格式。如下

这个为什么啊？不知道为什么会出现这个情况
我贴出代码。请大牛们帮我看看呐。
初学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/csm/shea_util.ncl"
load "./shapefile_mask_data.ncl"
begin
f1 = addfile("d:/ERA-JRA/ERA/column.nc", "r")
;print(f1)
;printVarSummary(f1)
wpc1=f1->tcwv(0:407,:,:)
lat=f1->latitude
lon=f1->longitude
scale_factor=0.001119838395708893
add_offset=36.77590352980346
wpc=wpc1*scale_factor+add_offset
wpc!0="time"
wpc!1="lat"
wpc!2="lon"
wpcClm=clmMonTLL(wpc) ; 计算月平均 (12,nlat,mlon)
printVarSummary(wpcClm)
Wseason =wpcClm(lat|:,lon|:,month|:) ; 转换维度
springW1 =Wseason(:,:,2:4) ;春季
springW =dim_sum(springW1)
summerW1 =Wseason(:,:,5:7) ;夏季
summerW =dim_sum(summerW1)
autumnW1 =Wseason(:,:,8:10) ;春季
autumnW =dim_sum(autumnW1)
winterW1 =Wseason(:,:,0:1) ;冬季
winterW2 =Wseason(:,:,11)
winterW3 =Wseason(:,:,2:4)
do k=0,1
winterW3(:,:,k)=winterW1(:,:,k)
end do
winterW3(:,:,2)=winterW2(:,:)
winterW=dim_sum(winterW3)
;春季
springW!0="lat"
springW!1="lon"
springW&lat = fspan(90,-90,241) ; (lat_start, lat_end, lat_num) 指定纵坐标分别对应的纬度
springW&lon = fspan(0,359.25,480) ; (lon,start, lon_end, lon_num) 指定横坐标分别对应的经度
springW&lat@units = "degrees_north" ; 赋予lat方向信息，这个是必须有的
springW&lon@units = "degrees_east" ; 赋予lon方向信息，这个是必须有的
;夏季
summerW!0="lat"
summerW!1="lon"
summerW&lat = fspan(90,-90,241) ; (lat_start, lat_end, lat_num) 指定纵坐标分别对应的纬度
summerW&lon = fspan(0,359.25,480) ; (lon,start, lon_end, lon_num) 指定横坐标分别对应的经度
summerW&lat@units = "degrees_north" ; 赋予lat方向信息，这个是必须有的
summerW&lon@units = "degrees_east" ; 赋予lon方向信息，这个是必须有的
;秋季
autumnW!0="lat"
autumnW!1="lon"
autumnW&lat = fspan(90,-90,241) ; (lat_start, lat_end, lat_num) 指定纵坐标分别对应的纬度
autumnW&lon = fspan(0,359.25,480) ; (lon,start, lon_end, lon_num) 指定横坐标分别对应的经度
autumnW&lat@units = "degrees_north" ; 赋予lat方向信息，这个是必须有的
autumnW&lon@units = "degrees_east" ; 赋予lon方向信息，这个是必须有的
;冬季
winterW!0="lat"
winterW!1="lon"
winterW&lat = fspan(90,-90,241) ; (lat_start, lat_end, lat_num) 指定纵坐标分别对应的纬度
winterW&lon = fspan(0,359.25,480) ; (lon,start, lon_end, lon_num) 指定横坐标分别对应的经度
winterW&lat@units = "degrees_north" ; 赋予lat方向信息，这个是必须有的
winterW&lon@units = "degrees_east" ; 赋予lon方向信息，这个是必须有的
;************画图******************
shp_filename="D:/ERA-JRA/studyarea/zong.shp"
data_mask0=shapefile_mask_data(springW,shp_filename,True)
data_mask1=shapefile_mask_data(summerW,shp_filename,True)
data_mask2=shapefile_mask_data(autumnW,shp_filename,True)
data_mask3=shapefile_mask_data(winterW,shp_filename,True)
minlat = 30
maxlat = 60
minlon = 45
maxlon = 110
;---Start the graphics
wks = gsn_open_wks("ps","wpc-each")
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
res@mpMinLatF = minlat
res@mpMaxLatF = maxlat
res@mpMinLonF = minlon
res@mpMaxLonF = maxlon
;---Create contours over map.
gsn_define_colormap(wks,"rainbow")
res@cnFillPalette = "MPL_gist_rainbow"
res@mpMinLatF=minlat
res@mpMaxLatF=maxlat
res@mpMinLonF=minlon
res@mpMaxLonF=maxlon
res@cnFillOn=True
res@cnLinesOn=False
res@lbBoxLinesOn=False
res@pmLabelBarWidthF=0.8
res@cnFillMode="RasterFill"
res@cnRasterSmoothingOn=True
;--Resources for polyline
lnres = True
lnres@gsLineColor = "black"
lnres@gsLineThicknessF = 1.5 ; 2x thickness
map_mask0 = gsn_csm_contour_map(wks,data_mask0,res)
line_mask0 = gsn_add_shapefile_polylines(wks, map_mask0, shp_filename, lnres)
draw(map_mask0)
frame(wks)
map_mask1 = gsn_csm_contour_map(wks,data_mask1,res)
line_mask1 = gsn_add_shapefile_polylines(wks, map_mask1, shp_filename, lnres)
draw(map_mask1)
frame(wks)
map_mask2 = gsn_csm_contour_map(wks,data_mask2,res)
line_mask2 = gsn_add_shapefile_polylines(wks, map_mask2, shp_filename, lnres)
draw(map_mask2)
frame(wks)
map_mask3 = gsn_csm_contour_map(wks,data_mask3,res)
line_mask3 = gsn_add_shapefile_polylines(wks, map_mask3, shp_filename, lnres)
draw(map_mask3)
frame(wks)
end

复制代码

longlivehj · 发表于 2014-7-14 19:27:17

把所有多边形合成为一个后再试试。

阿吉木柯 · 发表于 2014-7-14 18:53:00

上面两张图都是画的春季的（springW），下面的代码是第二张图的。

兰溪之水 · 发表于 2014-7-14 19:14:11

问题应该出现在load "./shapefile_mask_data.ncl"
你这个函数对点判断有问题，你看缺掉在边界附近就知道了~

阿吉木柯 · 发表于 2014-7-17 09:44:58

兰溪之水发表于 2014-7-14 19:14
问题应该出现在load "./shapefile_mask_data.ncl"
你这个函数对点判断有问题，你看缺掉在边界附近就知道了 ...

确实是这个的原因，函数没问题，是我的shp文件出了点问题，现在已经解决了。谢谢啊！

阿吉木柯 · 发表于 2014-7-17 09:45:17

longlivehj 发表于 2014-7-14 19:27
把所有多边形合成为一个后再试试。

确实是我的shp文件出了问题，谢谢啊！解决了。

阿吉木柯 · 发表于 2014-7-17 10:02:02

目前我还有l两个个问题。

year

season

第一张是年平均，可以看到背景颜色是白色，第二张是季节平均，但是背景是灰色，画图res都是一样的，为什么会这样呢？

还有就是，上面图的里海，巴尔喀什湖等这些都不是我的shp。我的shp文件只包括国家，里海这些信息怎么在画图当中略去呢，就是不要这些信息，怎么设置呢？

阿吉木柯 · 发表于 2014-7-17 10:06:10

@兰溪之水 @longlivehj

longlivehj · 发表于 2014-7-17 10:27:15

阿吉木柯发表于 2014-7-17 10:02
目前我还有l两个个问题。

地图的问题是因为你设置了res@mpDataBaseVersion = "MediumRes" ；
背景的问题，最好能把你最新版本的脚本发上来看看。

阿吉木柯 · 发表于 2014-7-17 11:16:42

res@mpDataBaseVersion = "MediumRes" ；这条命令不是高分辨率？我注释掉之后图还是没有变化呢。还是有里海等标志。

后来仔细看了两张图的res，发现那个白色背景的res中多了一条

res@cnFillMode="RasterFill"
res@cnRasterSmoothingOn=True

但这个不是不平滑的命令么？
我注释掉之后，就变成灰色背景了。是不是这个不平滑的命令让灰色的不平滑之后变成了看起了是白色，实际上还不是白色的呢。

下面贴我的代码帮我看下啊。
这个是年平均的

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"
load "./shapefile_mask_data.ncl"
begin
f1 = addfile("d:/ERA-JRA/ERA/column.nc", "r")
wpc1=f1->tcwv(0:407,:,:)
lat=f1->latitude
lon=f1->longitude
scale_factor=0.001119838395708893
add_offset=36.77590352980346
wpc=wpc1*scale_factor+add_offset
wpc!0="time"
wpc!1="lat"
wpc!2="lon"
wpcClm=clmMonTLL(wpc) ; 计算月平均 (12,nlat,mlon)
printVarSummary(wpcClm)
Wyear =wpcClm(lat|:,lon|:,month|:)
yearW1 =Wyear(:,:,0:11)
yearW =dim_sum(yearW1)
;赋坐标
yearW!0="lat"
yearW!1="lon"
yearW&lat = fspan(90,-90,241) ; (lat_start, lat_end, lat_num) 指定纵坐标分别对应的纬度
yearW&lon = fspan(0,359.25,480) ; (lon,start, lon_end, lon_num) 指定横坐标分别对应的经度
yearW&lat@units = "degrees_north" ; 赋予lat方向信息，这个是必须有的
yearW&lon@units = "degrees_east" ; 赋予lon方向信息，这个是必须有的
;************画图******************
shp_filename="D:/ERA-JRA/studyarea/zong.shp"
data_mask0=shapefile_mask_data(yearW,shp_filename,True)
minlat = 30
maxlat = 60
minlon = 45
maxlon = 110
;---Start the graphics
wks = gsn_open_wks("ps","column-year")
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
;---Create contours over map.
gsn_define_colormap(wks,"rainbow")
res@cnFillPalette = "temp_19lev"
res@mpMinLatF=minlat
res@mpMaxLatF=maxlat
res@mpMinLonF=minlon
res@mpMaxLonF=maxlon
res@cnFillOn=True
res@cnLinesOn=False
res@lbBoxLinesOn=False
; res@pmLabelBarWidthF=0.8
res@cnFillMode="RasterFill"
res@cnRasterSmoothingOn=True
;--Resources for polyline
lnres = True
lnres@gsLineColor = "black"
lnres@gsLineThicknessF = 1.5 ; 2x thickness
map_mask0 = gsn_csm_contour_map(wks,data_mask0,res)
line_mask0 = gsn_add_shapefile_polylines(wks, map_mask0, shp_filename, lnres)
draw(map_mask0)
asciiwrite("data_mask0.txt",data_mask0)
end

复制代码

下面这个是季节平均的

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"
load "./shapefile_mask_data.ncl"
begin
f1 = addfile("d:/ERA-JRA/ERA/column.nc", "r")
;print(f1)
;printVarSummary(f1)
wpc1=f1->tcwv(0:407,:,:)
lat=f1->latitude
lon=f1->longitude
scale_factor=0.001119838395708893
add_offset=36.77590352980346
wpc=wpc1*scale_factor+add_offset
wpc!0="time"
wpc!1="lat"
wpc!2="lon"
wpcClm=clmMonTLL(wpc) ; 计算月平均 (12,nlat,mlon)
printVarSummary(wpcClm)
Wseason =wpcClm(lat|:,lon|:,month|:) ; 转换维度
springW1 =Wseason(:,:,2:4) ;春季
springW =dim_sum(springW1)
summerW1 =Wseason(:,:,5:7) ;夏季
summerW =dim_sum(summerW1)
autumnW1 =Wseason(:,:,8:10) ;春季
autumnW =dim_sum(autumnW1)
winterW1 =Wseason(:,:,0:1) ;冬季
winterW2 =Wseason(:,:,11)
winterW3 =Wseason(:,:,2:4)
do k=0,1
winterW3(:,:,k)=winterW1(:,:,k)
end do
winterW3(:,:,2)=winterW2(:,:)
winterW=dim_sum(winterW3)
;春季
springW!0="lat"
springW!1="lon"
springW&lat = fspan(90,-90,241) ; (lat_start, lat_end, lat_num) 指定纵坐标分别对应的纬度
springW&lon = fspan(0,359.25,480) ; (lon,start, lon_end, lon_num) 指定横坐标分别对应的经度
springW&lat@units = "degrees_north" ; 赋予lat方向信息，这个是必须有的
springW&lon@units = "degrees_east" ; 赋予lon方向信息，这个是必须有的
;夏季
summerW!0="lat"
summerW!1="lon"
summerW&lat = fspan(90,-90,241) ; (lat_start, lat_end, lat_num) 指定纵坐标分别对应的纬度
summerW&lon = fspan(0,359.25,480) ; (lon,start, lon_end, lon_num) 指定横坐标分别对应的经度
summerW&lat@units = "degrees_north" ; 赋予lat方向信息，这个是必须有的
summerW&lon@units = "degrees_east" ; 赋予lon方向信息，这个是必须有的
;秋季
autumnW!0="lat"
autumnW!1="lon"
autumnW&lat = fspan(90,-90,241) ; (lat_start, lat_end, lat_num) 指定纵坐标分别对应的纬度
autumnW&lon = fspan(0,359.25,480) ; (lon,start, lon_end, lon_num) 指定横坐标分别对应的经度
autumnW&lat@units = "degrees_north" ; 赋予lat方向信息，这个是必须有的
autumnW&lon@units = "degrees_east" ; 赋予lon方向信息，这个是必须有的
;冬季
winterW!0="lat"
winterW!1="lon"
winterW&lat = fspan(90,-90,241) ; (lat_start, lat_end, lat_num) 指定纵坐标分别对应的纬度
winterW&lon = fspan(0,359.25,480) ; (lon,start, lon_end, lon_num) 指定横坐标分别对应的经度
winterW&lat@units = "degrees_north" ; 赋予lat方向信息，这个是必须有的
winterW&lon@units = "degrees_east" ; 赋予lon方向信息，这个是必须有的
;;;;;*****画图*****
shp_filename="D:/ERA-JRA/studyarea/zong.shp"
data_mask0=shapefile_mask_data(springW,shp_filename,True)
data_mask1=shapefile_mask_data(summerW,shp_filename,True)
data_mask2=shapefile_mask_data(autumnW,shp_filename,True)
data_mask3=shapefile_mask_data(winterW,shp_filename,True)
minlat = 30
maxlat = 60
minlon = 45
maxlon = 110
;---Start the graphics
wks = gsn_open_wks("ps","column-season-zong")
res = True
res@gsnDraw=False
res@gsnFrame=False
res@gsnMaximize = True ; maximize plot in frame
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
;---Create contours over map.
gsn_define_colormap(wks,"rainbow")
res@cnFillPalette = "temp_19lev"
res@mpMinLatF=minlat
res@mpMaxLatF=maxlat
res@mpMinLonF=minlon
res@mpMaxLonF=maxlon
res@cnFillOn=True
res@cnLinesOn=False
res@cnLineLabelsOn=False
; res@cnLevelSelectionMode="ManualLevels"
; res@cnMinLevelValF=20
; res@cnMaxLevelValF=280
; res@cnLevelSpacingF=10
res@lbBoxLinesOn=False
res@lbLabelBarOn=False ;关闭各自的labelbar
;--Resources for polyline
lnres = True
lnres@gsLineColor = "black"
lnres@gsLineThicknessF = 1.5 ; 2x thickness
map_masks=new(4,graphic)
map_masks(0) = gsn_csm_contour_map(wks,data_mask0,res)
line_mask0 = gsn_add_shapefile_polylines(wks, map_masks(0), shp_filename, lnres)
map_masks(1) = gsn_csm_contour_map(wks,data_mask1,res)
line_mask1 = gsn_add_shapefile_polylines(wks, map_masks(1), shp_filename, lnres)
map_masks(2) = gsn_csm_contour_map(wks,data_mask2,res)
line_mask2 = gsn_add_shapefile_polylines(wks, map_masks(2), shp_filename, lnres)
map_masks(3) = gsn_csm_contour_map(wks,data_mask3,res)
line_mask3 = gsn_add_shapefile_polylines(wks, map_masks(3), shp_filename, lnres)
pres =True ;设置panel的resources
pres@gsnMaximize= True
pres@gsnPanelLabelBar=True
pres@pmLabelBarWidthF=0.8
pres@lbBoxLinesOn=False
; pres@txString="Precipitation(mm)"
pres@gsnPanelFigureStrings=ispan(1,4,1)
gsn_panel(wks,map_masks,(/2,2/),pres)
panelres=True
panelres@gsnPanelLableBar=True
end

复制代码

谢谢@longlivehj啊！现在就是希望背景是白色以及没有里海等标志，只有我的shp文件。

		自动登录	找回密码
密码			立即注册