我用cmorph融合的逐小时降水资料画的图很奇怪，请问是怎么回事？

pkn789

                               
登录/注册后可看大图

我用的是cmorph融合的逐小时降水资料，是grd格式的，想画105-120E，25-40N的降水，结果画出来的图很奇怪，下面是我的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 "./gsn_csm.ncl"
;============================================================================


; TEMPLATE TO PLOT A MAP OF THE PRECIPITABLE WATER.

begin


;SET RESOURCES:
res = True
res@gsnAddCyclic           =False
;res@mpGridAndLimbOn        =True
;res@mpGridLineColor        ="Black"
;res@mpGridLineDashPattern  =4
;res@mpLandFillColor        = "white"
;res@mpOceanFillColor       = "white"
;res@mpInlandWaterFillColor = "white"
res@cnFillOn               =True
res@cnLinesOn              ="False"
res@cnLineLabelsOn         ="False"
res@cnSmoothingOn          =True
res@cnLevelSelectionMode   = "ExplicitLevels"
res@cnLevels               =  (/ 6, 12, 25,50,70,100 /)
res@cnExplicitLabelBarLabelsOn = True
res@lbLabelStrings         = (/"6","12","25","50","70","100"/)
res@mpLimitMode = "LatLon"
res@mpDataBaseVersion       = "MediumRes"
res@mpDataSetName           = "Earth..4"
res@mpAreaMaskingOn         = True
res@mpMaskAreaSpecifiers    = (/"China:Henan"/)
res@mpMinLatF              =25
res@mpMaxLatF              =40
res@mpMinLonF              =105
res@mpMaxLonF              =120
res@mpOutlineBoundarySets   = "NoBoundaries"

res@mpLimitMode = "LatLon"
  res@pmTickMarkDisplayMode = "Always"
;reading files

                 tim = (/"SEVP_CLI_CHN_MERGE_CMP_PRE_HOUR_GRID_0.10-2016070818.grd",\
                       "SEVP_CLI_CHN_MERGE_CMP_PRE_HOUR_GRID_0.10-2016070819.grd",\
                             "SEVP_CLI_CHN_MERGE_CMP_PRE_HOUR_GRID_0.10-2016070820.grd",\
                       "SEVP_CLI_CHN_MERGE_CMP_PRE_HOUR_GRID_0.10-2016070821.grd",\
                    "SEVP_CLI_CHN_MERGE_CMP_PRE_HOUR_GRID_0.10-2016070822.grd",\
                   "SEVP_CLI_CHN_MERGE_CMP_PRE_HOUR_GRID_0.10-2016070823.grd",\
                 "SEVP_CLI_CHN_MERGE_CMP_PRE_HOUR_GRID_0.10-2016070900.grd",\
                 "SEVP_CLI_CHN_MERGE_CMP_PRE_HOUR_GRID_0.10-2016070901.grd",\
               "SEVP_CLI_CHN_MERGE_CMP_PRE_HOUR_GRID_0.10-2016070902.grd",\
               "SEVP_CLI_CHN_MERGE_CMP_PRE_HOUR_GRID_0.10-2016070903.grd",\
               "SEVP_CLI_CHN_MERGE_CMP_PRE_HOUR_GRID_0.10-2016070904.grd",\
                "SEVP_CLI_CHN_MERGE_CMP_PRE_HOUR_GRID_0.10-2016070905.grd",\
               "SEVP_CLI_CHN_MERGE_CMP_PRE_HOUR_GRID_0.10-2016070906.grd",\
                    "SEVP_CLI_CHN_MERGE_CMP_PRE_HOUR_GRID_0.10-2016070907.grd",\
                   "SEVP_CLI_CHN_MERGE_CMP_PRE_HOUR_GRID_0.10-2016070908.grd"/)


nlat=80
nlon=70
ntim=15
lon = fspan(105,120,nlon)
lat = fspan(25,40,nlat)
lat!0 = "lat"
lat@long_name = "latitude"
lat@units = "degrees_north"
lon!0 = "lon"
lon@units = "degrees_east"
lon@long_name = "longitude"
olon = fspan(105,120,70)
olat = fspan(25,40,80)


r = new((/ntim,nlat,nlon/),"double")
do t = 0, ntim-1
rain = fbindirread(tim(t),0,-1,"double")
k=0
do i = 0, nlat-1
  do j = 0, nlon-1
  r(t,i,j)= rain(k)
   k = k+1
  end do
end do
end do

;r!0 = "time"
;r&time = tim
;r@units = "mm/hour"
;r@_FillValue = -999.0
;n=num(ismissing(r))
;print("n="+n)
;printVarSummary(r)

res@gsnAddCyclic = False

wks = gsn_open_wks("pdf","/home/wangl/pkn/station/sta")
cmap=(/(/1.,1.,1./),(/0.,0.,0./),(/1.,1.,1./),(/.26,.88,.95/),(/.153,.706,.733/),(/.08,.61,.26/),\
     (/.28,.90,.49/),(/.89,.88,.18/),(/.77,.75,.10/),(/.65,.45,.16/),(/.612,.125,.137/),\
     (/.8392,.2118,.2275/),(/.9059,.0745,.5098/),(/.65,.10,.86/),(/.24,.16,.45/)/)
gsn_define_colormap(wks,cmap)
;gsn_draw_colormap(wks)

res@gsnLeftString   = "rain"
res@gsnRightString  = "[mm~N~]"
res@tiMainString    = "2016-07-08-18:00 - 2016-07-09-08:00"

frain=r(0,:,:)
do t = 1,ntim-1
  srain = r(t,:,:)
  arain =srain + frain
  frain = arain
end do
print(frain)
filo="crain"
fbinrecwrite(filo,-1,frain)
nlat=80
nlon=70
ntim=15
lon = fspan(105,120,nlon)
lat = fspan(25,40,nlat)
frain!0 = "lat"
frain!1 = "lon"
frain&lat = lat
frain&lon = lon
frain@units = "mm/hour"
frain@_FillValue = -999.00
n=num(ismissing(frain))
print("n="+n)
printVarSummary(frain)
plot = gsn_csm_contour_map(wks,frain,res)
draw(plot)
   frame(wks)
end
;==================================================

werewolf

数据没有正确读进来。。。

pkn789

嗯嗯谢谢，我换了个脚本可以了

有趣的人

pkn789 发表于 2017-9-27 17:48
嗯嗯谢谢，我换了个脚本可以了

老师您好，我也在用cmorph融合的逐小时降水资料画图，但是NCL实在掌握太差，请问您最后用的什么脚本可以画出来啊？能分享下吗老师？

pkn789

;*******************************************************************************  
  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"
  load "/home/wangl/public/cnmap/cnmap.ncl"

;*************************************************
begin

raintot = conform_dims((/440,700/), 0.0 , -1)

do i = 0818,900
name=i
       if(isfilepresent("/home/wangl/pkn/station/SEVP_CLI_CHN_MERGE_CMP_PRE_HOUR_GRID_0.10-2016070"+name+".grd.nc")) then
       a=addfile("/home/wangl/pkn/station/SEVP_CLI_CHN_MERGE_CMP_PRE_HOUR_GRID_0.10-2016070"+name+".grd.nc","r")
  


rain =  a->GPCP  
copy_VarMeta(rain, raintot)
;print(dimsizes(x))  

raintot = rain + raintot
printVarSummary(raintot)

end if
end do



;********************************************************************************

    wks = gsn_open_wks("pdf","0818-0900")
   cmap = RGBtoCmap("15colorsrgb.txt")
gsn_define_colormap(wks,cmap)

res = True
  res@gsnDraw                = False
  res@gsnFrame               = False
; res@gsnAddCyclic           = False      ; not cycle
  res@gsnMaximize            = False     ; maximize

  res@tmXTOn = False
res@tmYROn = False
;res@gsnPaperOrientation = "Landscape"

res@pmTickMarkDisplayMode = "Always"

;######### SET FOR THE MAP  #########
  

   
mpres=res
mpres@mpProjection = "LambertConformal"
mpres@mpLambertMeridianF  = 113
mpres@mpLimitMode = "Corners"
mpres@mpLambertParallel1F = 30.
mpres@mpLambertParallel2F = 60.
mpres@mpLeftCornerLatF       = 33.6
   mpres@mpRightCornerLatF      = 38.3
   mpres@mpLeftCornerLonF       = 109.9
   mpres@mpRightCornerLonF      = 116.5


     
mpres@mpFillOn = False
mpres@mpOutlineOn = False
mpres@mpOutlineSpecifiers = "China:states"
mpres@mpGeophysicalLineThicknessF = 2
mpres@mpDataBaseVersion = "MediumRes"
mpres@mpDataSetName = "Earth..4"
mpres@mpProvincialLineThicknessF =4
mpres@tmXBLabelFontHeightF=0.022
mpres@tmXBLabelFontThicknessF=0.022
  
  
  

res@mpFillOn = False
res@mpOutlineOn = False
  res@cnFillOn = True     
  res@cnLinesOn =False  
  ;res@mpGridAndLimbOn   = False
res@mpFillOn = False
res@mpOutlineOn = False
  res@pmTickMarkDisplayMode = "Always"  
  res@cnLevelSelectionMode = "ExplicitLevels"
  res@cnLevels             = (/ 10,20,30,40,50,60,70,80,90,100,110,120,130,140,150/)
; res@cnFillColors         = (/"White","White","DarkOliveGreen1", \
                                    ;    "DarkOliveGreen3","Chartreuse", \
                                       ; "Chartreuse3","Green","ForestGreen", \
                                      ;  "Yellow","Orange","Red","Violet"/)

;res@cnMinLevelValF  = 0   
  ;res@cnMaxLevelValF  = 10     
  ;res@cnLevelSpacingF = 0.2she
;############   OTHERS   #############
  ;res@lbOrientation = "Vertical"
  ;res@pmLabelBarOrthogonalPosF = -0.01
  ;res@lbPerimOn            = False
  ; res@lbLabelPosition="left"
   res@gsnLeftString   = ""
   res@gsnRightString  = ""
   res@lbBoxLinesOn       = True

   res@tiMainString    = "0812-0818"
   res@tiMainFontHeightF  = 0.020              ; smaller title
res@tmBorderThicknessF = 3.0

   wrf_smooth_2d( raintot, 11 )
  
map = gsn_csm_contour_map(wks,raintot,res)
;######   ADD MAP OF CHINA   #############

  cnres           = True
  cnres@china     = True       ;draw china map or not
  cnres@river     = False     ;draw changjiang&huanghe or not
  cnres@province  = True       ;draw province boundary or not
  cnres@nanhai    = False       ;draw nanhai or not
  cnres@diqu      = False     ;draw diqujie or not   

  plot = add_china_map(wks,map,cnres)





;************************************************************************************

  draw(map)
  frame(wks)

  end
  
;**************************************************************************************
  

pkn789

有趣的人 发表于 2019-4-8 14:45
老师您好，我也在用cmorph融合的逐小时降水资料画图，但是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/wrf/WRFUserARW.ncl"
  load "cnmap.ncl"

;*************************************************
begin

raintot = conform_dims((/440,700/), 0.0 , -1)

do i = 0818,900
name=i
       if(isfilepresent("SEVP_CLI_CHN_MERGE_CMP_PRE_HOUR_GRID_0.10-2016070"+name+".grd.nc")) then
       a=addfile("SEVP_CLI_CHN_MERGE_CMP_PRE_HOUR_GRID_0.10-2016070"+name+".grd.nc","r")
  


rain =  a->GPCP  
copy_VarMeta(rain, raintot)
;print(dimsizes(x))  

raintot = rain + raintot
printVarSummary(raintot)

end if
end do


;********************************************************************************

    wks = gsn_open_wks("pdf",0818-0900")
   cmap = RGBtoCmap("15colorsrgb.txt")
gsn_define_colormap(wks,cmap)

res = True
  res@gsnDraw                = False
  res@gsnFrame               = False
; res@gsnAddCyclic           = False      ; not cycle
  res@gsnMaximize            = False     ; maximize

  res@tmXTOn = False
res@tmYROn = False
;res@gsnPaperOrientation = "Landscape"

res@pmTickMarkDisplayMode = "Always"

;######### SET FOR THE MAP  #########
  

   
mpres=res
mpres@mpProjection = "LambertConformal"
mpres@mpLambertMeridianF  = 113
mpres@mpLimitMode = "Corners"
mpres@mpLambertParallel1F = 30.
mpres@mpLambertParallel2F = 60.
mpres@mpLeftCornerLatF       = 33.6
   mpres@mpRightCornerLatF      = 38.3
   mpres@mpLeftCornerLonF       = 109.9
   mpres@mpRightCornerLonF      = 116.5


     
mpres@mpFillOn = False
mpres@mpOutlineOn = False
mpres@mpOutlineSpecifiers = "China:states"
mpres@mpGeophysicalLineThicknessF = 2
mpres@mpDataBaseVersion = "MediumRes"
mpres@mpDataSetName = "Earth..4"
mpres@mpProvincialLineThicknessF =4
mpres@tmXBLabelFontHeightF=0.022
mpres@tmXBLabelFontThicknessF=0.022
  
  
  

res@mpFillOn = False
res@mpOutlineOn = False
  res@cnFillOn = True     
  res@cnLinesOn =False  
  ;res@mpGridAndLimbOn   = False
res@mpFillOn = False
res@mpOutlineOn = False
  res@pmTickMarkDisplayMode = "Always"  
  res@cnLevelSelectionMode = "ExplicitLevels"
  res@cnLevels             = (/ 10,20,30,40,50,60,70,80,90,100,110,120,130,140,150/)
; res@cnFillColors         = (/"White","White","DarkOliveGreen1", \
                                    ;    "DarkOliveGreen3","Chartreuse", \
                                       ; "Chartreuse3","Green","ForestGreen", \
                                      ;  "Yellow","Orange","Red","Violet"/)

;res@cnMinLevelValF  = 0   
  ;res@cnMaxLevelValF  = 10     
  ;res@cnLevelSpacingF = 0.2she
;############   OTHERS   #############
  ;res@lbOrientation = "Vertical"
  ;res@pmLabelBarOrthogonalPosF = -0.01
  ;res@lbPerimOn            = False
  ; res@lbLabelPosition="left"
   res@gsnLeftString   = ""
   res@gsnRightString  = ""
   res@lbBoxLinesOn       = True

   res@tiMainString    = "0812-0818"
   res@tiMainFontHeightF  = 0.020              ; smaller title
res@tmBorderThicknessF = 3.0

   wrf_smooth_2d( raintot, 11 )
  
map = gsn_csm_contour_map(wks,raintot,res)
;######   ADD MAP OF CHINA   #############

  cnres           = True
  cnres@china     = True       ;draw china map or not
  cnres@river     = False     ;draw changjiang&huanghe or not
  cnres@province  = True       ;draw province boundary or not
  cnres@nanhai    = False       ;draw nanhai or not
  cnres@diqu      = False     ;draw diqujie or not   

  plot = add_china_map(wks,map,cnres)





;************************************************************************************

  draw(map)
  frame(wks)

  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'
warning:ContourPlotInitialize: scalar field is constant; no contour lines will appear; use cnConstFEnableFill to enable fill
warning:ContourPlotSetValues: Data values out of range of levels set by EXPLICITLEVELS mode