- 积分
- 93
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2014-5-13
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看，没有帐号？立即注册
x
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 "/space/gpfsl01/users/hanshuai/liujunjian/program/china/cnmap/cnmap.ncl"
begin
x =asciiread("x.txt",(/90/),"float")
y =asciiread("y.txt",(/90/),"float")
bias =asciiread("bias.txt",(/90/),"string")
era_bias =asciiread("era_bias.txt",(/90/),"string")
;---------------------------------------------------------------------------
nLAT = 2880 ;纬度
nLON = 5760 ;经度
file_time =(/"2015060100","2015060106","2015060112","2015060118"/)
do i=0,0 ;-------------------循环时次
;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++需要修改的地方
;pic_name=file_time(i) +"_gfs_isccp_pw_bias"
pic_name=file_time(i) +"_era_isccp_pw_bias"
;pic_name=file_time(i) +"_gfs_era_pw_bias"
print(pic_name)
;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++需要修改的地方
file_name_isccp =file_time(i)+"_isccp_pw.nc"
file_name_gfs =file_time(i)+"_gfs_pw.nc"
file_name_era =file_time(i)+"_era_pw.nc"
;---------------------------------------读取文件
f_isccp =addfile(file_name_isccp,"r")
f_gfs =addfile(file_name_gfs,"r")
f_era =addfile(file_name_era,"r")
;---------------------------------------读取变量
var_isccp =f_isccp->pw(:,:)
var_gfs =f_gfs->PW(:,:)
var_era =f_era->tcwv(:,:)
;---------------------------------------差运算
gfs_isccp =var_gfs/1.0-var_isccp*10
era_isccp =var_era/1.0-var_isccp*10
gfs_era =gfs_isccp-era_isccp
;---------------------------------------公共范围
newlon = fspan(0,359.5,nLON)
newlat = fspan(-90 ,90 ,nLAT)
newlat@units = "degrees_north"
newlon@units = "degrees_east"
;---------------------------------------GFS-ISCCP
gfs_isccp!0 ="latitude"
gfs_isccp!1 = "longitude"
gfs_isccp&latitude = newlat
gfs_isccp&longitude = newlon
;---------------------------------------ERA-ISCCP
era_isccp!0 ="latitude"
era_isccp!1 = "longitude"
era_isccp&latitude = newlat
era_isccp&longitude = newlon
;---------------------------------------GFS-ERA
gfs_era!0 ="latitude"
gfs_era!1 = "longitude"
gfs_era&latitude = newlat
gfs_era&longitude = newlon
;---------------------------------------输出NC数据
;fout_gfs=addfile(file_time(i)+"_gfs_isccp_pw.nc","c")
;fout_gfs->pw=gfs_isccp
;fout_era=addfile(file_time(i)+"_era_isccp_pw.nc","c")
;fout_era->pw=era_isccp
wks = gsn_open_wks("png",pic_name)
gsn_define_colormap(wks,"BlAqGrYeOrReVi200")
setvalues NhlGetWorkspaceObjectId()
"wsMaximumSize" : 300000000
end setvalues
res = True
res@cnLineLabelsOn = False
res@cnFillDrawOrder = "Predraw" ; fill and lines before map
;res@lbLabelBarOn = False
;res@cnFillPalette = "BlueRed" ; Blue-Red colormap
res@gsnAddCyclic = False
res@tiMainString = pic_name
res@gsnMaximize = True
res@gsnDraw = False
res@gsnFrame = False
res@gsnLeftString = "Total column water vapour"
res@gsnRightString = "kg/m**2"
;>--------------------------------------------<
; set for the map
;>--------------------------------------------<
res@mpMinLatF = 17.
res@mpMaxLatF = 55.
res@mpMinLonF = 72.
res@mpMaxLonF = 136.
res@cnLevelSelectionMode = "ExplicitLevels"
;res@cnLevels = ispan(-40,40,10);
res@cnLevels = ispan(-20,20,5);
res@mpFillOn = True
res@mpOutlineOn = False ; Use outlines from shapefile
res@cnFillDrawOrder = "PreDraw"
res@mpDataBaseVersion = "MediumRes"
res@mpDataSetName = "Earth..4"
res@mpAreaMaskingOn = True
res@mpMaskAreaSpecifiers = (/"China","Taiwan","Disputed area between India and China","India:Arunachal Pradesh"/)
res@mpLandFillColor = "white"
res@mpInlandWaterFillColor = "white"
res@mpOceanFillColor = "white"
res@mpOutlineBoundarySets = "NoBoundaries"
;>--------------------------------------------<
; set for the plot
res@cnFillOn = True
res@cnLinesOn = False
res@cnLevelSpacingF = 2.
res@gsnSpreadColors = True
res@lbLabelAutoStride = True
;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++需要修改的地方
;map=gsn_csm_contour_map(wks,gfs_isccp,res)
map=gsn_csm_contour_map(wks,era_isccp,res)
;map=gsn_csm_contour_map(wks,gfs_era,res)
;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++需要修改的地方
res@tmYROn = False ; Turn off right and top tick marks
res@tmXTOn = False ; Turn off right and top tick marks
; add China map
;>------------------------------------------------------------<
cnres = True
cnres@china = True ;draw china map or not
cnres@river = True ;draw changjiang&huanghe or not
cnres@province = True ;draw province boundary or not
cnres@nanhai = True ;draw nanhai or not
cnres@diqu = False ;draw diqujie or not
chinamap = add_china_map(wks,map,cnres)
;----------------------------------------------
txres = True
txres@txFontHeightF = 0.010
mkres = True
mkres@gsMarkerIndex = 17 ; Filled circle
mkres@gsMarkerSizeF = 1
mkres@gsMarkerColor = "red" ; polymarker color
;----------------------------------------------
do j = 0,89
txt1 =gsn_add_text(wks,map,era_bias(j),x(j),y(j),txres)
;gsn_polymarker(wks,map,x(j),y(j),mkres)
point =gsn_add_polygon(wks,map,x(j),y(j),mkres)
end do
draw(map)
frame(wks)
end do
end
主要是是这个 do j = 0,89
txt1 =gsn_add_text(wks,map,era_bias(j),x(j),y(j),txres)
;gsn_polymarker(wks,map,x(j),y(j),mkres)
point =gsn_add_polygon(wks,map,x(j),y(j),mkres)
end do
point 画不出来
warning ：
warning:TransformPostDraw: tfPolyDrawList element 3482 is invalid
warning:TransformPostDraw: tfPolyDrawList element 3483 is invalid
warning:TransformPostDraw: tfPolyDrawList element 3484 is invalid
warning:TransformPostDraw: tfPolyDrawList element 3485 is invalid
warning:TransformPostDraw: tfPolyDrawList element 3486 is invalid
warning:TransformPostDraw: tfPolyDrawList element 3487 is invalid
warning:TransformPostDraw: tfPolyDrawList element 3488 is invalid
warning:TransformDataPolygon, not enough points for a polygon
，求大神支招。
|
|
窥视卡
雷达卡
发表于 2016-1-28 22:53:29
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
抢沙发
千斤顶
显身卡
楼主