- 积分
- 6467
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2015-11-6
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
请教各位,根据http://bbs.06climate.com/forum.php?mod=viewthread&tid=11417这个网址所介绍的绘制站点资料图的方法,我用753站的降水资料也进行了绘制,脚本如下,但是如图很明显是错误的,值完全超出了中国边界,重新看了下脚本也找不出错误,希望各位大神帮忙指点一下,实在是当局者迷了。
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/WRFUserARW.ncl"
begin
fname = "/nuist/p/work/yongqing/liaoyue/rain_20.txt"
fname1 = "/nuist/p/work/yongqing/liaoyue/rain_21.txt"
fname2 = "/nuist/p/work/yongqing/liaoyue/rain_22.txt"
lines = asciiread(fname,-1,"string")
lines1 = asciiread(fname1,-1,"string")
lines2 = asciiread(fname2,-1,"string")
lat = stringtofloat(str_get_field(lines(0:),2," "))
lon = stringtofloat(str_get_field(lines(0:),3," "))
pwv = stringtofloat(str_get_field(lines(0:),4," "))
rain = stringtofloat(str_get_field(lines(0:),5," "))/10.
rain1 = stringtofloat(str_get_field(lines1(0:),5," "))/10.
rain2 = stringtofloat(str_get_field(lines2(0:),5," "))/10.
; This second file is not so tricky. The 2D lat/lon data is sorted with
; lat values first, and then lon values.
olon = new(66,"float");
olat = new(40,"float");
data = new((/40,66/),"float")
data1 = new((/40,66/),"float")
data2 = new((/40,66/),"float")
do i=0,65
olon(i) =72+i
end do
do l=0,39
olat(l) = 17+l
end do
olon!0 = "lon"
olon@long_name = "lon"
olon@units = "degrees-east"
olon&lon = olon
olat!0 = "lat"
olat@long_name = "lat"
olat@units = "degrees_north"
olat&lat = olatrain@_FillValue =3270.
rain1@_FillValue =3270.
rain2@_FillValue =3270.
rscan = (/10,5,3/) ;连续的有效半径大小,最大为10,依次递减
data = obj_anal_ic_deprecated(lon,lat,rain,olon,olat,rscan, False) ;Creanm插值
data1 = obj_anal_ic_deprecated(lon,lat,rain1,olon,olat,rscan, False)
data2 = obj_anal_ic_deprecated(lon,lat,rain2,olon,olat,rscan, False)
wks=gsn_open_wks("eps","obs_rain")
gsn_define_colormap(wks,"WhiteBlueGreenYellowRed")
res=True
;res@gsnDraw=False
;res@gsnFrame=False
;res@gsnMaximize=True
res@gsnAddCyclic = False ;由于我们的数据不是循环地球一周的,因此必须把这个置否
res@gsnLeftString="" ;去掉左上角标记
res@gsnRightString=""
res@tmYLLabelFontHeightF=0.01
res@tmXBLabelFontHeightF=0.01
;res@vpWidthF=0.34
;res@vpHeightF=0.25
res@tmXTOn=False
res@tmYROn=False
res@mpDataSetName = "Earth..4"
res@mpDataBaseVersion = "MediumRes"
res@mpOutlineOn = True
res@mpOutlineSpecifiers = (/"China:states","Taiwan"/)
res@mpMinLatF = 17 ; Asia limits
res@mpMaxLatF = 55
res@mpMinLonF = 72
res@mpMaxLonF = 136
;默认的底图投影方式是等经纬度投影,画出来的中国地图比较扁,我们常看到的中国地图,投影方式是兰伯特投影,因此需要对投影方式进行修改
res@mpProjection = "LambertConformal" ;兰伯特投影
res@mpLambertMeridianF = 110.0
res@mpLimitMode = "LatLon"
res@mpLambertParallel1F = .001 ;Default: .001 ;可以自己改一改,看看投影有什么不同,挺有趣的
res@mpLambertParallel2F = 89.999 ;Default: 89.999
;最后将填充区域设定在中国行政区域图之内,如果使用默认效果,等值线会对整个矩形区域填充颜色,因此需要去掉中国边境范围外的填充颜色
res@mpAreaMaskingOn = True ;使能填充覆盖
res@mpMaskAreaSpecifiers = (/"China:states","Taiwan"/) ;China:states
res@mpOceanFillColor = 0 ;用白色填充海洋 0是colormap的索引值
res@mpInlandWaterFillColor = 0 ;用白色填充内陆湖水
res1=res
res1@cnFillOn=True
res1@cnLinesOn=False
res1@cnLineLabelsOn=False
res1@lbLabelBarOn=True
res1@lbLabelFontHeightF=0.015
;res1@cnLevelSelectionMode="ExplicitLevels"
;res1@cnLevels=(/20,40,60,80,100,120,140,160,180,200,220,240/)
;res1@cnFillColors=(/0,20,40,60,80,100,120,140,160,180,200,220,240/)
plot=gsn_csm_contour_map(wks,data1,res1)
end
|
-
|