- 积分
- 10656
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2013-10-10
- 最后登录
- 1970-1-1
|
发表于 2014-7-17 11:09:35
|
显示全部楼层
本帖最后由 longlivehj 于 2014-7-17 11:38 编辑
就在那个程序上改的。把diri改成bin文件所在目录。
diri = "./" ; input directory
两个输出，一个是nc，文件名与bin除了后缀，其它完全相同；另外一个是ps，提供预览。
; --------------------------3B42RT-----------------------------
; Header:
; Each file starts with a header that is one 2-byte-integer
; row in length, or 2880 bytes.
; Following the header, 4 data fields appear:
; precipitation (2-byte integer), precipitation_error (2-byte integer)
; source (1-byte integer), precipitation_uncalibrated (2-byte integer)
; All fields are 1440x480 grid boxes (0-360E,60N-S).
; The first grid box center is at (0.125E,59.875N)
; The binary file is 'big endian'
; -------------------------------------------------------------
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"
;***************************************************************
; User Input
;***************************************************************
begin
; INPUT
diri = "./" ; input directory
fili = systemfunc("ls " + diri + "/*.bin") ; binary uncompressed
; OUTPUT
netCDF = True ; generate netCDF file
PLOT = True ; generate plots
do i = 0, dimsizes(fili) - 1
sfx = get_file_suffix(fili(i), 0)
if (netCDF) then
ncDir = "./" ; directory for netCDF output
ncFil = sfx@fBase + ".nc" ; netCDF name output
end if
if (PLOT) then
pltDir = "./" ; directory for plot output
pltName = sfx@fBase ; plot name root
pltType = "ps"
end if
nlat = 480
mlon = 1440
lhead = 2880 ; characters (lhead/2=1440 type short)
lheads = lhead/2 ; header length as type 'short'
print("lhead="+lhead+" nlat="+nlat+" " \
+"mlon="+mlon+" nlat*mlon="+nlat*mlon)
prcScale = 0.01
errScale = 0.01
;***************************************************************
; End User Input
;***************************************************************
; Miscellaneous: Parse the file name to extract strings/ date info
;***************************************************************
filc = tochar( fili(i) )
yyyy_str = tostring(filc( 7:10)) ; "yyyy"
mon_str = tostring(filc(11:12)) ; "mon"
dd_str = tostring(filc(13:14)) ; "dd"
hh_str = tostring(filc(15:16)) ; "hr"
yyyy = toint( yyyy_str ) ; full year
mm = toint( mon_str ) ; month
dd = toint( dd_str ) ; day of month
hh = toint( hh_str )
print("yyyy="+yyyy+" mm="+mm+" dd="+dd+" hh="+hh)
;***************************************************************
; Read BigEndian binary file
;***************************************************************
; Read Header (2880 characters (character=byte in length)
;***************************************************************
setfileoption("bin","ReadByteOrder","BigEndian")
hdrc = fbindirread (diri+fili(i), 0, lhead , "character") ; header
hdrs = tostring(hdrc) ; create string
print(hdrs)
;;print(hdrc)
;***************************************************************
; Read precip (short); skip header; read type short; reshape & scale
; Read error (short); skip header and precip; read err; reshape & scale
;***************************************************************
works = fbindirread (diri+fili(i), 0, -1, "short")
works = where(works.le.toshort(-31998) .or. works.ge.toshort(31999)\
,toshort(-31999), works)
works@_FillValue = toshort(-31999)
;;print(works(lhead/2:lhead/2+100))
starts = lhead/2 ; length of header in terms of "short"
lasts = starts+nlat*mlon-1
prc = onedtond(works(starts:lasts)*prcScale, (/nlat,mlon/) )
starts = lasts+1
lasts = starts+nlat*mlon-1
err = onedtond(works(starts:lasts)*errScale, (/nlat,mlon/) )
;***************************************************************
; Read source (byte); skip header (char), precip and error (short)
;***************************************************************
startb = lhead + 4*nlat*mlon ; 4=2_bytes(short)*2_variables
lastb = startb+nlat*mlon-1
workb = fbindirread (diri+fili(i), 0, -1, "byte")
workb@_FillValue = tobyte(0)
src = onedtond(workb(startb:lastb),(/nlat,mlon/) )
delete(workb) ; no longer needed
;***************************************************************
; Read uncalibrated precip (short); skip header, precip, error, source
;***************************************************************
starts = lhead/2 + 2*nlat*mlon + (nlat*mlon)/2
lasts = starts+nlat*mlon-1
print("prcu: starts="+starts+" lasts="+lasts)
prcu = onedtond(works(starts:lasts)*prcScale, (/nlat,mlon/) )
delete(works) ; no longer needed
;***************************************************************
; Add meta data
;***************************************************************
prc@units = "mm/hr"
prc@long_name = "TRMM_3B42RT: Hourly Rain Rate"
;prc@_FillValue = ; ??
printMinMax(prc,0)
err@units = "mm/hr"
err@long_name = "TRMM_3B42RT: Hourly Rain Rate"
;err@_FillValue =
printMinMax(err,0)
src@long_name = "Data Source"
src@source = "0=no observation, 1=AMSU, 2=TMI, 3=AMSR, 4=SSMI, "+\
"5=SSMIS, 6=MHS, 30=AMSU&MHS avg, 31=conical avg, "+\
"50=IR, 1,2,3,4,5,6+100=Sparce Sample"
prcu@units = "mm/hr"
prcu@long_name = "TRMM_3B42RT: Uncalibrated Hourly Rain Rate"
prcu@information = "Multi-satellite precipitation before this calibration"
;prcu@_FillValue =
printMinMax(prcu,0)
;*****************************************************
; Create TRMM coordinate variables. See README
;*****************************************************
lat = 59.875 - ispan(0,nlat-1,1)*0.25
lon = ispan(0,mlon-1,1)*0.25
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
;***************************************************************
prc!0 = "lat" ; 1st ... name the dimensions
prc!1 = "lon"
prc&lat = lat ; create coordinate variable
prc&lon = lon
err!0 = "lat"
err!1 = "lon"
err&lat = lat
err&lon = lon
src!0 = "lat"
src!1 = "lon"
src&lat = lat
src&lon = lon
prcu!0 = "lat" ; 1st ... name the dimensions
prcu!1 = "lon"
prcu&lat = lat ; create coordinate variable
prcu&lon = lon
printVarSummary(prc)
printMinMax(prc, 0)
printVarSummary(err)
printMinMax(err, 0)
printVarSummary(src)
printMinMax(src, 0)
printVarSummary(prcu)
printMinMax(prcu, 0)
;************************************************
; Create plot ?
;************************************************
if (PLOT) then
wks = gsn_open_wks(pltType, pltName)
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 = "RasterFill" ; Raster Mode
res@cnLineLabelsOn = False ; Turn off contour lines
res@cnLevelSelectionMode = "ExplicitLevels"
res@cnMissingValFillPattern = 0
res@cnMissingValFillColor = "black"
res@lbOrientation = "vertical" ; vertical label barb's
res@lbLabelFontHeightF = 0.012 ; change font size
res@pmLabelBarOrthogonalPosF = -0.01 ; move a bit to left
res@pmLabelBarWidthF = 0.1
res@mpMinLatF = -50.
res@mpMaxLatF = 50.
res@mpCenterLonF = 210.
res@mpFillOn = False
res@mpOutlineOn = True
res@mpOutlineBoundarySets = "National" ; turn on country boundaries
;res@mpShapeMode = "FreeAspect"
;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@cnLevels = (/0.1,1,2,3,4,5,7.5,10,15,20/) ; "mm/3hr"
res@tiMainString = fili(i)
prc = where (prc.lt.0, prc@_FillValue, prc) ; bogus
plot = gsn_csm_contour_map_ce(wks,prc, res)
end if
;************************************************
; Create netCDF ?
; Recommend to always create a 'time' dimension
;************************************************
if (netCDF) then
ntim = 1
tunits = "hours since 1997-01-01 00:00:0.0"
time = (/cd_inv_calendar(yyyy,mm,dd,hh, 0,0d0,tunits, 0)/)
time!0 = "time"
date = yyyy*1000000 + mm*10000 + dd*100 + hh
date!0 = "time"
date@units = "yyyymmddhh"
nline = inttochar(10) ; new line character
globeAtt = 1
globeAtt@title = "TRMM_3B42"
globeAtt@ftp = nline + \
"http://disc.sci.gsfc.nasa.gov/data/datapool/TRMM_DP/01_Data_Products/02_Gridded/06_3-hour_Gpi_Cal_3B_42" +nline
globeAtt@ID = "3B42RT"
globeAtt@description = nline + \
"????" + nline
globeAtt@creation_date= systemfunc ("date" )
NCFILE = ncFil
system ("/bin/rm -f " + NCFILE) ; remove any pre-exist file
ncdf = addfile(NCFILE,"c")
;setfileoption(ncdf, "definemode", True)
fileattdef( ncdf, globeAtt ) ; create the global [file] attributes
dimNames = (/"time", "lat", "lon" /)
dimSizes = (/ ntim , nlat, mlon /)
dimUnlim = (/ True , False, False /)
filedimdef(ncdf, dimNames , dimSizes, dimUnlim )
filevardef (ncdf, "time" , typeof(time), getvardims(time) )
filevarattdef(ncdf, "time", time)
filevardef (ncdf, "lat", typeof(lat), getvardims(lat))
filevarattdef(ncdf, "lat", lat)
filevardef (ncdf, "lon", typeof(lon), getvardims(lon))
filevarattdef(ncdf, "lon", lon)
filevardef (ncdf, "date" , typeof(date), getvardims(date) )
filevarattdef(ncdf, "date", date)
filevardef (ncdf, "PRC" , typeof(prc) , (/ "time", "lat", "lon" /) )
filevarattdef (ncdf, "PRC" , prc)
filevardef(ncdf, "RELERR" , typeof(err), (/ "time", "lat", "lon" /) )
filevarattdef(ncdf, "RELERR", err)
ncdf->time = (/ time /)
ncdf->lat = (/ lat /)
ncdf->lon = (/ lon /)
ncdf->date = (/ date /)
ncdf->PRC(0:0,:,:) = (/ prc /)
ncdf->RELERR(0:0,:,:) = (/ err /)
end if ; netCDF
end do
end
|
|
窥视卡
雷达卡
发表于 2014-7-17 10:21:29
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
抢沙发
千斤顶
显身卡
楼主
发表于 2014-7-17 11:25:48