爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

搜索
查看: 10519|回复: 6

[作图] ncl 如何画跨越0度经线的图(大西洋65W-15E)

[复制链接]
发表于 2017-10-5 06:20:08 | 显示全部楼层 |阅读模式

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

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

x
用ncep的sst画热带大西洋的eof(25N-25S,65W-15E)的时候图片出来没有等值线,因为是跨0度经线,可能需要设置吧,不知道怎么设置参数我自己的脚本

begin
f1=addfile("NCEP1950-2014sstDJF.nc","r")

sst=f1->sst
;sst@missing_value=-1e+30
b=sst
copy_VarCoords(sst(0,:,:),b(0,:,:))
b!0="time"
time=ispan(1950,2014,1)
b&time=time
b&time@long_name="year"
b&time@units="year"
b@missing_value=-1e+30
x=b
;x=dim_standardize_n(b, 1, 0)
copy_VarCoords(b,x)
printVarSummary(x)
optEOF = False
optETS = False      
neval  =3
eof= eofunc_Wrap(x({lat|-25:25},{lon|295:15},{time|:}),neval,optEOF)
eof_ts=eofunc_ts_Wrap (x({lat|-25:25},{lon|295:15},{time|:}),eof, optETS)
eof_ts=eof_ts/100
eof=eof*100

wks = gsn_open_wks("png","TAeof1")
res2=True
res2@gsnDraw = False
res2@gsnFrame = False
res2@gsnAddCyclic        = False
;res2@gsnAddCyclic        = True
;res2@mpCenterLonF         = -30
res2@mpLimitMode       = "LatLon"
;res2@mpCenterLonF         = 180  
res2@mpMinLatF         = -25        
res2@mpMaxLatF         = 25
res2@mpMinLonF         = -65
res2@mpMaxLonF         = 15
res2@cnFillOn       = True
;res2@cnFillPattern ="Explicit"
;res2@cnLevelSelectionMode="ManualLevels"
;res2@cnMinLevelValF=2.4
;res2@cnMaxLevelValF=0.8
;res2@cnLevelSpacingF=0.4
res2@cnLinesOn =False
res2@tmXTOn =False
res2@tmYROn=False
;res2@cnFillColors=(/146,161,177,193,209,226,242/)
res2@cnSmoothingOn      =True
res2@cnRasterSmoothingOn= True
res2@tmXBLabelFontHeightF  = 0.015
res2@tmYLLabelFontHeightF  = 0.015
res2@lbOrientation = "Vertical"

;res2@cnMissingValFillColor="white"


rts2           = True
rts2@gsnDraw = False
rts2@gsnFrame = False
rts2@vpHeightF = 0.30        ; Changes the aspect ratio
rts2@vpWidthF  = 0.75
rts2@vpXF      = 0.80        ; change start locations
rts2@vpYF      = 0.75        ; the plot
rts2@tiYAxisString = ""                  
rts2@tiXAxisString = "year"
rts2P                      = True            ; modify the panel plot
rts2P@gsnMaximize          = True            ; large format
;rts2@trYMaxF=2
;rts2@trYMinF=-1.2

rts2@gsnYRefLine           = 0.              ; reference line   
rts2@gsnXYBarChart         = True            ; create bar chart
rts2@gsnAboveYRefLineColor = "red"           ; above ref line fill red
rts2@gsnBelowYRefLineColor = "blue"          ; below ref line fill blue

rts2@tmXBMode = "Explicit"
rts2@tmXBValues=(/"1950","1960","1970","1980","1990","2000","2010"/)
  rts2@tmXBLabels=(/"1950","1960","1970","1980","1990","2000","2010"/)
  ;rts2@tmYLMode="Explicit"
;rts2@tmYLValues=(/"-1.6","-1.2","-0.8","-0.4","0.0","0.4","0.8","1.2","1.6"/)
;rts2@tmYLLabels=(/"-1.6","-1.2","-0.8","-0.4","0.0","0.4","0.8","1.2","1.6"/)
;rts2@tiMainString    = "AP_Ts_02:1950-2014"
rts2@gsnLeftString  = "EOF "+(0+1)

rts2@tmXTOn =False
rts2@tmYROn=False
;rts2@gsnRightString = sprintf("%5.1f", eof@pcvar(1)) +"%"


plots = new(2,graphic)

  ;res2@tiMainString = "1950_2014hwf-std-eof1"
  res2@gsnLeftString  = "EOF1 "
  ;cmap2 = read_colormap_file("14")
;res2@cnFillColors = cmap2
  ;gsn_define_colormap(wks,"14")
res2@gsnRightString = sprintf("%5.1f", eof@pcvar(0)) +"%"
  plots(0) = gsn_csm_contour_map(wks,eof(0,:,:),res2)



plots(1)= gsn_csm_xy (wks,time,eof_ts(0,:),rts2)   






  pres = True
  ;pres@tiMainString = "1950_2014hwf_eof"
  pres@gsnMaximize = True
  pres@gsnPanelDebug=True
  pres@gsnPanelRowSpec = True
  pres@gsnPanelLabelBar = False
  pres@lbBoxLinesOn = False
  gsn_panel(wks,plots,(/1,2/),pres)
end






数据的经纬度格式

数据的经纬度格式

我出的图

我出的图
密码修改失败请联系微信:mofangbao
 楼主| 发表于 2017-10-5 06:56:21 | 显示全部楼层
自问自答了。。。。是因为经度设置是0.5-359.5,聪明的办法不会,就用笨办法。。。。将数组拆成两个再合并,重新赋值经纬度的坐标变量
begin
f1=addfile("NCEP1950-2014sstDJF.nc","r")

sst=f1->sst
;sst@missing_value=-1e+30
b=sst
copy_VarCoords(sst(0,:,:),b(0,:,:))
b!0="time"
time=ispan(1950,2014,1)
b&time=time
b&time@long_name="year"
b&time@units="year"
b@missing_value=-1e+30
lon=fspan(0.5,359.5,360)
lat=fspan(89.5,-89.5,180)
b&lon=lon
b&lat=lat
b&lon@long_name="lon"
b&lon@units="degrees_east"
b&lat@long_name="lat"
b&lat@units="degrees_north"
x=b
x1=b({lat|-25:25},{lon|295:360},{time|:})
x2=b({lat|-25:25},{lon|0:15},{time|:})
x3=new((/50,80,65/),"float")
x3(:,0:64,:)=x1
x3(:,65:79,:)=x2
x3!0="latitude"
x3!1="longitude"
x3!2="time"
x3&time=time
x3&time@long_name="year"
x3&time@units="year"
x3@missing_value=-1e+30
longitude=fspan(-65,15,80)
latitude=fspan(-25,25,50)
x3&longitude=longitude
x3&latitude=latitude
x3&longitude@long_name="longitude"
x3&longitude@units="degrees_east"
x3&latitude@long_name="latitude"
x3&latitude@units="degrees_north"

printVarSummary(x3)
printVarSummary(x1)
printVarSummary(x2)
;x=dim_standardize_n(b, 1, 0)
copy_VarCoords(b,x)
printVarSummary(x)
optEOF = False
optETS = False      
neval  =3
eof= eofunc_Wrap(x3({latitude|-25:25},{longitude|:},{time|:}),neval,optEOF)
eof_ts=eofunc_ts_Wrap (x3({latitude|-25:25},{longitude|:},{time|:}),eof, optETS)

如果有更好的办法求指教
密码修改失败请联系微信:mofangbao
发表于 2017-10-6 08:06:22 | 显示全部楼层
eof(0,:,:)你没有设置经纬度信息的原因不吧
密码修改失败请联系微信:mofangbao
发表于 2020-1-4 10:52:30 | 显示全部楼层
出现了相同问题,有没有更好的办法呢?
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

发表于 2020-7-8 18:02:41 | 显示全部楼层
我也想知道 咋办呀!
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

发表于 2021-5-8 21:08:59 | 显示全部楼层
songwei 发表于 2017-10-6 08:06
eof(0,:,:)你没有设置经纬度信息的原因不吧

hhh 在说啥呢hhhh
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

发表于 2021-5-10 09:47:02 | 显示全部楼层
lonflip
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表