;*************************************************************** ; cmorph_4.ncl ; ; Concepts illustrated: ; - Reading big endian binary files ; - Reading records written by a fortran *direct access* write ; - Reading CMORPH 8km data ; - Adding meta data (attributes and coordinates [time, lat, lon]) ; - Explicitly setting contour levels and colors ;*************************************************************** ;; ftp://ftp.cpc.ncep.noaa.gov/precip/global_CMORPH/README.cmorph.8km_30minute ;; ;; Each direct access record is a 4948 x 1649 CHARACTER*1 (use FORTRAN ichar ;; command to retrieve interger value for all parameters) array with grid ;; increments of 0.072756669 degrees of longitude and 0.072771377 of latitude, ;; which is apporoximately 8 km at the equator. The arrays are oriented from ;; North to South, beginning from latitude 59.963614N and from West to EAST from ;; longitude 0.036378335E. ;; ;; Missing data are denoted by values of 255. ;; ;; Note that the precipitation estimates have been scaled, and when multiplied by ;; "0.2", the data units are "mm/hour". ;; ;; For GrADS users, a descriptor ("ctl") file: CMORPH_8km-30-minute.ctl has been ;; provided. However, since the data are in CHARACTER*1 words (bytes) the parameters ;; after each variable (-1,40,1,-1 in our example "ctl file") are system dependent. ;; Our example is for an SGI system. ;; ;; Each file contains 6 records. The 1st 3 records pertain to the top half ;; of the hour (00-29 minus after the hour) and the last 3 records are for the ;; bottom half of the hour. Within each group: ;; ;; - the 1st record contains the CMORPH precipitation estimates ;; ;; - the 2nd record contains the time (in half hour units) since the most ;; recent microwave pass. Note that since we do both a forward & ;; backward interpolation in time, the nearest time may be prior to ;; the file time stamp or after it. ;; ;; - the 3rd record contains an ID that tells the satellite from which the last ;; microwave observation was made which can be interpretted by the following ;; table (as of the time of the last update of this documentation): ;; ;; 13 = DMSP-13 (SSM/I instrument) ;; 14 = DMSP-14 ( " " ) ;; 15 = DMSP-15 ( " " ) ;; 16 = DMSP-16 (SSMIS instrument, coming soon) ;; 17 = DMSP-17 ( " " ) ;; 18 = DMSP-18 ( " " ) ;; 115 = NOAA-15 (AMSU-B " ) ;; 116 = NOAA-16 ( " " ) ;; 117 = NOAA-17 ( " " ) ;; 118 = NOAA-18 (MHS ) ;; 119 = NOAA-19 ( " " ) ;; 151 = METOP-A ( " " ) ;; 201 = TRMM (TMI " ) ;; 211 = AQUA (AMSR-E " ) ;; ;; ;; Each direct access record is a 4948 x 1649 CHARACTER*1 (use FORTRAN ichar ;; command to retrieve interger value for all parameters) array with grid ;; increments of 0.072756669 degrees of longitude and 0.072771377 of latitude, ;; which is apporoximately 8 km at the equator. The arrays are oriented from ;; North to South, beginning from latitude 59.963614N and from West to EAST from ;; longitude 0.036378335E. ;; ;; Missing data are denoted by values of 255. ;; ;; Note that the precipitation estimates have been scaled, and when multiplied by ;; "0.2", the data units are "mm/hour". ;******************************************************************************** diri = "./" fili = "CMORPH_8KM-30MIN_2015111101" ; direct access (flat) pthi = diri+fili ;;setfileoption("bin","ReadByteOrder","BigEndian") ; *not* needed; input are type byte dlim = "_-" ; string delimiter nfld = str_fields_count(fili, dlim) ; nfld=4 ymdh = toint(str_get_field(fili, 4, dlim)) ; yyyymmddhh yyyy = ymdh/1000000 mdh = ymdh - yyyy*1000000 mm = mdh/10000 dh = mdh - mm*10000 dd = dh/100 hh = dh-dd*100 tunits= "hours since 2000-01-01 00:00:00" ; arbitrary date time = cd_inv_calendar(yyyy,mm,dd,hh, 0, 0,tunits, 0) time!0= "time" ntim = 1 nlat = 1649 nlon = 4948 lat = 59.963614d - ispan(0,nlat-1,1)*0.072771377d ; N->S lat!0 = "lat" lat@units = "degrees_north" printMinMax(lat,0) lon = 0.036378335d + ispan(0,nlon-1,1)*0.072756669d lon!0 = "lon" lon@units = "degrees_east" printMinMax(lon,0) ;---- Read 'top-half' of the hour (00-29 minus after the hour) ; This is the 1st record (recnum=0) prc_u = fbindirread(pthi,0,(/ntim,nlat,nlon/),"ubyte") prc_u@_FillValue = toubyte(255) printMinMax(prc_u,0) prc = where(ismissing(prc_u), 1e20, prc_u*0.2) prc@_FillValue = 1e20 delete(prc_u) ; no longer needed prc@long_name = "CMORPH 8km" prc@units = "mm/hr" prc!0 = "time" prc!1 = "lat" prc!2 = "lon" prc&time = time prc&lat = lat prc&lon = lon printVarSummary(prc) printMinMax(prc,0) PLOT = True if (PLOT) then pltType = "png" ; send graphics to PNG file pltDir = "./" pltName = "cmorph" wks = gsn_open_wks(pltType, pltDir+pltName) res = True ; plot mods desired res@gsnMaximize = True 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@cnLinesOn = False ; Turn off contour lines res@cnLineLabelsOn = False ; Turn off contour lines res@cnLevelSelectionMode = "ExplicitLevels" res@cnLevels = (/0.1,1,2.5,5,10,15,20,25,50,75/) ; "mm/day" res@cnFillPalette = (/"Snow","PaleTurquoise","PaleGreen","SeaGreen3" ,"Yellow" \ ; contour colors ,"Orange","HotPink","Red","Violet", "Purple", "Brown"/) ; one more color than contour levels res@cnMissingValFillPattern = 0 ; make 'missing' black res@cnMissingValFillColor = "black" res@lbLabelBarOn = False ; turn off individual cb's res@mpMinLatF = -60. ; CMORPH limits [approx] res@mpMaxLatF = 60. res@mpCenterLonF = 210. res@mpFillOn = False ;res@mpShapeMode = "FreeAspect" ;res@vpWidthF = 0.8 ;res@vpHeightF = 0.4 nt = 0 res@tiMainString = fili plot = gsn_csm_contour_map(wks,prc(nt,:,:), res) end if ; PLOT