TRMM 数据读取出错，求指导

hhumuer

基于NCL官网提供的TRMM数据读取脚本，自己稍微改了一下，可是运行提示fatal:syntax error: line 65 in file hdfForm.ncl before or near

hhumuer

;*********** Load Libraries ************************************
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"
;**************************************************************
begin
;***************************************************************
; User Input
;***************************************************************
                                             ; INPUT
   diri   = "./"                             ; input directory
   fili   = "SATE_L3_TRM_MUTDS_MWB_3B42_GLB_V6-20030702-03.hdf"            ; binary uncompressed

   PLOT   = True                             ; generate plots

   if (PLOT) then
       pltDir  = "./"                        ; directory for plot output
       pltName = "SATE_L3_TRM_3B42_GLB_V6-20030702-03"      ; plot name root
       pltType = "pdf"                                          
   end if

;***************************************************************
; End User Input
;***************************************************************
; Miscellaneous: Parse the file name to extract strings/ date info
;***************************************************************
   filc        = stringtochar( fili )
   yy_str      = chartostring(filc(34: 37))     ; "yy"
   mon_str     = chartostring(filc(38: 39))     ; "mon"
   dd_str      = chartostring(filc(40:41))     ; "dd"
   hh_str      = chartostring(filc(43:44))    ; "hh"

   yyyy        = stringtoint( yy_str )        ; full year
   mm          = stringtoint( mon_str )       ; month
   dd          = stringtoint( dd_str )        ; day  of month
   hh          = stringtoint( hh_str )

;***************************************************************
; Read hdf
; The "scan" dimension corresponds to "time".
; Note the HDF dimension order is (time,lon,lat)
; We want (time,lat,lon) order to match the model
; Use NCL's dimension reordering syntax to accomplish this
;***************************************************************

   f      = addfile (diri+fili, "r")
   work   = f->precipitation                   ; (scan, longitude, latitude)
   prc    = work(scan|:,latitude|:,longitude|:)
   delete(work)


;***************************************************************
; Add meta data: See above README
;***************************************************************
   prc@_FillValue  = -9999.
   prc@units       = "mm/hr"
   ;prc@long_name   = "TRMM_3B42 precip"

;*****************************************************
; Create TRMM coordinate variables. See README
;*****************************************************
   nlat       = 400
   mlon       = 1440
   lat        = ispan(0,nlat-1,1)*0.25 -  49.875    ; 纬度换算
   lon        = ispan(0,mlon-1,1)*0.25 - 179.875    ； 经度换算

   lat@long_name = "latitude"
   lat@units  = "degrees_north"
   lat!0      = "lat"
   lat&lat    =  lat

   lon@long_name = "longitude"
   lon@units  = "degrees_east"
   lon!0      = "lon"
   lon&lon    =  lon
***************************************************************
; Associate the spatial coordinates with variables
; Rename the dimension to match the model.
; Convenience only, not required.
;***************************************************************
   prc!0      = "time"
   prc!1      = "lat"                   ; 1st ... name the dimensions
   prc!2      = "lon"
   prc&lat    =  lat                    ; create coordinate variable
   prc&lon    =  lon                    
;***************************************************************
; Simple data exploration:
;    Are there missing data?
;    Count the number of missing values in each variable
;
;    Calculate weighted areal averages: ignore missing grid points
;    Print results
;***************************************************************

   nMsg_prc = num(ismissing(prc )) ；ismissing returns true for every element of the input that contains a  missing value
                                   ; num, counts the number of True values in the input
   rad      = 4.*atan(1.0)/180.    ; atan, computes the inverse tangent of numeric types

   clat     = cos(lat*rad)    ; simple cosine weighting, cos, computes the cosine of numeric types

   prcAvg   = wgt_areaave(    prc, clat, 1.0, 0); wgt_areaave, calculates the area average of a quantity using weights.
   relerrAvg= wgt_areaave( relerr, clat, 1.0, 0); clat is the wgty containing the weights; 1.0, the wgtx

   print(" ")
   print("Number missing: nMsg_prc="+nMsg_prc)
   print(" ")
   print("TRMM_3B42: prcAvg="+prcAvg+" relerrAvg="+relerrAvg)
   print(" ")
   printMinMax(prc, False)      ; printMinMax, prints the minimum and maximum values of a variable
   printMinMax(relerr, True)
;************************************************
; Create plot ?
;************************************************
   if (PLOT) then
       plot    = new ( 2, "graphic"); new,creates an NCL variable.2 is the dimension size; "graphic"is the vartype.
   
       wks    = gsn_open_wks(pltType, pltDir+pltName)  ; open a workstation on which to draw graphics, type and name.
       colors = (/"white","black", "Snow"    \
               ,"PaleTurquoise","PaleGreen","SeaGreen3" ,"Yellow"  \
               ,"Orange","HotPink","Red","Violet", "Purple", "Brown"/)

       gsn_define_colormap(wks, colors)               ; generate new color map
   
       res                      = True     ; plot mods desired
       res@gsnDraw              = False    ; let gsn_panel do plotting
       res@gsnFrame             = False
       res@gsnMaximize          = True     ; make ps/eps/pdf large
     
       res@cnFillOn             = True     ; turn on color fill
       res@cnLinesOn            = False    ; turn of contour lines
       res@cnFillMode           = "CellFill"           ; Cell Mode
       res@cnFillMode           = "RasterFill"         ; Raster Mode
       res@cnLineLabelsOn       =  False       ; Turn off contour lines
       res@cnLevelSelectionMode = "ExplicitLevels"              
       res@cnMissingValFillPattern = 0
       res@cnMissingValFillColor   = "black"

       res@lbOrientation        = "horizontal"   ; vertical label barb's
       res@lbLabelFontHeightF   = 0.02        ; change font size,default, 0.02
       res@pmLabelBarOrthogonalPosF = 0.1     ; move a bit to left,default,0.02
       res@pmLabelBarWidthF     =  0.1       ; default, 0.15,0.6
       res@mpMinLatF            = 24.         ; CMORPH limits [approx]-50
       res@mpMaxLatF            =  40.         ; 50
      ;res@mpCenterLonF         = 115.         ; 210
       res@mpMinLonF            =  104.
       res@mpMaxLonF            =  126.
       res@mpFillOn             = False
       res@mpOutlineOn          = True
       res@mpOutlineBoundarySets  = "National"   ; turn on country boundaries
       res@mpDataBaseVersion="MediumRes"    ;加入中国省界边界线
       res@mpDataSetName="Earth..4"         ;加入中国省界边界线
       res@mpOutlineSpecifiers = "China:states"   ;加入中国省界边界线
       ;res@mpGeophysicalLineColor = "Black"
       ;res@mpProjection = "LambertConformal"  ; or LambertEqualArea

      ;res@mpShapeMode          = "FreeAspect" ; not preserve the map projection aspect ratio.
      ;res@vpWidthF             = 0.8
      ;res@vpHeightF            = 0.4
   
       res@cnLevels             = (/0.1,1,2.5,5,10,15,20,25,50,75/) ; "mm/3hr"
       res@gsnCenterString = "Areal Mean="+sprintf("%4.2f",  prcAvg)
      plot(0)    = gsn_csm_contour_map_ce(wks,prc(0,:,:), res)
       ;plot   = gsn_csm_contour_map_ce(wks,prc(0,:,:), res)

       resP = True
       resP@gsnMaximize         = True                ; make ps/eps/pdf large [no effect x11]
     ;;resP@gsnPaperOrientation = "Portrait"          ; force portrait
       resP@txString            = "TRMM_L3_3B42_"+":  "+yyyy+"-"+mm+"-"+dd+"_"+hh_str+":00:00"
       gsn_panel(wks,plot,(/2,1/),resP)               ; now draw as one plot
   end if

end
   

红尘_海

你的数据是V6还是V7？官网的程序可以运行V6数据，对V7数据不可以。

yrovl

一大清早就看到有一排“TRMM数据读取出错”，额。。。要不要一次发这么多篇呀！

hhumuer

红尘_海 发表于 2014-11-5 07:40
你的数据是V6还是V7？官网的程序可以运行V6数据，对V7数据不可以。

应该是V6的，因为官网的程序是可以读取，这个只是我自己改了一下，跑不起来，想知道是不是哪个地方没写对~~~

hhumuer

yrovl 发表于 2014-11-5 09:04
一大清早就看到有一排“TRMM数据读取出错”，额。。。要不要一次发这么多篇呀！

第一次发贴，囧，见谅啦~~~

longlivehj

lon        = ispan(0,mlon-1,1)*0.25 - 179.875    ； 经度换算
这句最后的注释用的是中文的分号！

hhumuer

longlivehj 发表于 2014-11-6 10:10
lon        = ispan(0,mlon-1,1)*0.25 - 179.875    ； 经度换算
这句最后的注释用的是中文的分号！

脚本的确有好多中文分号的原因，谢谢！。但是它自己原本加的一些注释后来运行时也提示错误，楼上的有碰到过这种情况吗？

longlivehj

hhumuer 发表于 2014-11-6 10:44
脚本的确有好多中文分号的原因，谢谢！。但是它自己原本加的一些注释后来运行时也提示错误，楼 ...

脚本我没有运行过。
如果有什么其它类型的错误，贴出来看看。

hhumuer

longlivehj 发表于 2014-11-6 10:54
脚本我没有运行过。
如果有什么其它类型的错误，贴出来看看。

就是不断提示注释的地方有问题，所以我把注释都去掉，现在可以跑了