积分 114
贡献
精华
在线时间 小时
注册时间 2019-2-25
最后登录 1970-1-1
4 金钱
本帖最后由 Helel 于 2019-5-31 11:18 编辑
改了几遍程序了,能添加的经纬度信息基本也都添加了,还是会报错,显示经纬度数据无单位,想请教各位一下原因代码如下
==========================================================
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"
begin
;yaogan
dir = "/home/data/"
year = "2003"
in1 = addfile(dir+"Ta_avg_monthly_"+year+"01.nc","r")
LST1=in1->LST
LST1!0 = "lat"
LST1!1 = "lon"
;printVarSummary(LST)
LST1&lat=fspan(25.01,45.00,2000)
LST1&lon=fspan(70.00,104.99,3500)
LST1&lat@units = "degrees_north"
LST1&lon@units = "degrees_east"
var = new((/12,2000,3500/),float)
do imo=1,12
smo = sprinti("%0.2i",imo)
in = addfile(dir+"Ta_avg_monthly_2003"+smo+".nc","r")
LST :=in->LST(:,:)
LST!0 = "lat"
LST!1 = "lon"
;printVarSummary(LST)
LST&lat=fspan(25.01,45.00,2000)
LST&lon=fspan(70.00,104.99,3500)
LST&lat@units = "degrees_north"
LST&lon@units = "degrees_east"
var(imo-1,:,:) = (/in->LST(:,:)/)
end do
var!0 = "time"
var!1 = "lat"
var!2 = "lon"
var&time = ispan(1,12,1)
var&lat=fspan(25.01,45.00,2000)
var&lon=fspan(70.00,104.99,3500)
varanl = dim_avg_n(var,0)
varanl!0 = "lat"
varanl!1 = "lon"
varanl&lat=fspan(25.01,45.00,2000)
varanl&lon=fspan(70.00,104.99,3500)
latgrd = fspan(28.25,36.25,45)
longrd = fspan(73.25,104.25,125)
vargrd=linint2(LST1&lon,LST1&lat,var,True,longrd,latgrd,0)
varanlgrd=linint2(LST1&lon,LST1&lat,varanl,True,longrd,latgrd,0)
;CN05
inu=addfile("/home/data/CN05.1_Tm_1961_2014_month_025x025.nc","r")
varinu=inu->tm(504:515,:,:)
latu=inu->lat
lonu=inu->lon
varu=new((/12,163,283/),float)
varu(0,:,:)=varinu(0,:,:)
varu(1,:,:)=varinu(1,:,:)
varu(2,:,:)=varinu(2,:,:)
varu(3,:,:)=varinu(3,:,:)
varu(4,:,:)=varinu(4,:,:)
varu(5,:,:)=varinu(5,:,:)
varu(6,:,:)=varinu(6,:,:)
varu(7,:,:)=varinu(7,:,:)
varu(8,:,:)=varinu(8,:,:)
varu(9,:,:)=varinu(9,:,:)
varu(10,:,:)=varinu(10,:,:)
varu(11,:,:)=varinu(11,:,:)
varu!0 = "time"
varu!1 = "lat"
varu!2 = "lon"
varu&time = ispan(1,12,1)
varu&lat = latu
varu&lon = lonu
varuanl = dim_avg_n(varu,0)
varuanl!0 = "lat"
varuanl!1 = "lon"
varuanl&lat = latu
varuanl&lon = lonu
diff = new((/12,163,283/),float)
diff = -9999.0
do kkd=0,11 ;month
do jjd=0,162 ;jini,jend
if ((latu(jjd).le. 36.25) .and. (latu(jjd).ge. 28.25)) then
jjw=doubletointeger((latu(jjd)-28.25)/0.25)
do iid= 0,282
if ( (lonu(iid).ge. 73.25) .and.(lonu(iid).le. 104.25)) then
iiw=doubletointeger((lonu(iid)-73.25)/0.25)
diff(kkd,jjd,iid)= vargrd(kkd,jjw,iiw)-varu(kkd,jjd,iid)
end if
end do
end if
end do
end do
diff!0 = "time"
diff!1 = "lat"
diff!2 = "lon"
diff&time = ispan(1,12,1)
diff&lat = latu
diff&lon = lonu
diff&lat@units = "degrees_north"
diff&lon@units = "degrees_east"
diff@_FillValue= -9999.0
wks = gsn_open_wks("ps","diff-TP-ATM-yaogan"+year)
setvalues NhlGetWorkspaceObjectId()
"wsMaximumSize" : 300000000
end setvalues
plotM = new(12,graphic)
plotU = new(12,graphic)
plotD = new(12,graphic)
gsn_define_colormap(wks,"BlueDarkRed18")
res = True
res@gsnDraw = False
res@gsnFrame = False
res@gsnSpreadColors = True
res@cnFillMode="RasterFill"
res@gsnLeftString = ""
res@gsnRightString = ""
res@gsnCenterString = ""
res@gsnLeftStringFontHeightF = 0.025
res@mpProjection = "LambertConformal" ; choose projection
res@mpLambertParallel1F = 30
res@mpLambertParallel2F = 60
res@mpLambertMeridianF = 100
res@mpLimitMode = "LatLon" ; choose range of map
res@mpMinLatF = 20 ;lat2d(0,0)
res@mpMinLonF = 60 ;lon2d(0,0)
res@mpMaxLatF = 50 ;lat2d(jlat-1,ilon-1)
res@mpMaxLonF = 110 ;lon2d(0,ilon-1)
res@mpGridAndLimbOn = True
res@mpFillOn = False
res@mpGridLineDashPattern =2
res@mpOutlineBoundarySets = "National"
res@gsnAddCyclic = False
res@pmTickMarkDisplayMode = "Always"
res@cnFillOn = True
res@cnLinesOn = False
res@cnLineLabelsOn = False
res@cnLevelSelectionMode = "ExplicitLevels" ;"ManualLevels"
res@cnLevels = (/-22.0,-20.0,-18.0,-16.0,-14.0,-12.0,-10.0,-8.0,-6.0,-4.0,-2.0,0.0,2.0,4.0,6.0,8.0,10.0,12.0,14.0,16.0,18.0,20.0,22.0/)
res@lbLabelBarOn = False
res@trGridType="TriangularMesh"
opt = True
;=================
pres = True
pres@gsnPanelLabelBar = True
pres@gsnPanelRowSpec = True
;==============================================================
do i=0,11
plotM(i) = gsn_csm_contour_map(wks,var(i,:,:),res)
end do
pres@txString = "yaoogan"
pres@gsnPanelFigureStrings =(/"Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"/)
pres@gsnPanelFigureStringsFontHeightF=0.01
pres@amJust = "TopLeft"
gsn_panel(wks,plotM,(/3,3,3,3/),pres)
;==============================================================
res@gsnLeftString = ""
do i=0,11
plotU(i) = gsn_csm_contour_map(wks,varu(i,:,:),res)
end do
pres@txString = "CN05"
pres@gsnPanelFigureStrings =(/"Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"/)
pres@gsnPanelFigureStringsFontHeightF=0.01
pres@amJust = "TopLeft"
gsn_panel(wks,plotU,(/3,3,3,3/),pres)
;======================================================
gsn_define_colormap(wks,"MPL_RdBu")
gsn_reverse_colormap(wks)
res@gsnLeftString = ""
delete(res@cnLevels)
;res@cnLevels = (/-16.0,-14.0,-12.0,-10.0,-8.0,-6.0,-4.0,-2.0,0.0,2.0,4.0,6.0,8.0,10.0,12.0,14.0,16.0/)
do i=0,11
plotD(i) = gsn_csm_contour_map(wks,diff(i,:,:),res)
end do
pres@txString = "Bias"
pres@gsnPanelFigureStrings =(/"Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"/)
pres@gsnPanelFigureStringsFontHeightF=0.01
pres@amJust = "TopLeft"
gsn_panel(wks,plotD,(/3,3,3,3/),pres)
end
=============================================================
报错信息如下:
(0) check_for_y_lat_coord: Warning: Data either does not contain
(0) a valid latitude coordinate array or doesn't contain one at all.
(0) A valid latitude coordinate array should have a 'units'
(0) attribute equal to one of the following values:
(0) 'degrees_north' 'degrees-north' 'degree_north' 'degrees north' 'degrees_N' 'Degrees_north' 'degree_N' 'degreeN' 'degreesN' 'deg north'
(0) check_for_lon_coord: Warning: Data either does not contain
(0) a valid longitude coordinate array or doesn't contain one at all.
(0) A valid longitude coordinate array should have a 'units'
(0) attribute equal to one of the following values:
(0) 'degrees_east' 'degrees-east' 'degree_east' 'degrees east' 'degrees_E' 'Degrees_east' 'degree_E' 'degreeE' 'degreesE' 'deg east'
我来回答