Re: problem to use lancos filter with NCL

Bonjour,

The response function of the 241-weight Lancos filter
[Duschon, 1979] may be seen at:
             http://www.cgd.ucar.edu/web/public_html/shea/lancos_241.gif

I suggest that you create a netCDF file containing the band passed
filtered data. This way you create the band passed data one time.

An untested sample script is attached. It
     [1] generates the weights
     [2] reads/reorders the OLR
     [3] removes the 12 harmonics like the Matthews reference [below]
     [4] applies the Lancos weights to create band passed data
     [5] creates a netCDF file

Bon Chance

; generate Lancos weights
     nwgt = 241
     fca = 1./200. ; 200 day cutoff
     fcb = 1./20. ; 20 day cutoff
     ihp = 2 ; bandpass
     sigma = 1.
     wgts = filwgts_lancos (nwgt, ihp, fca, fcb, sigma)

> On May 16, 2007, at 3:49 AM, temeching gladice joelle wrote:
>
>> I want to use ncl program of lancos filter to filter OLR satellite data
>> In fact, the purpose of my work, is to find the result of a
>> publication of
>>
>> ADRIAN J. MATTHEWS : Intraseasonal Variability over Tropical Africa
>> during Northern Summer.er. Journal of Climate, volume 17, 2427-2440.
>> A Lanczos filter used in this paper has 241-weights, and a
>> frequency between 20-200days.
>>
>> I am trying to filter out the band-pass-frequency signal and
>> using "wgt_runave" to compute this with lancos weights
>> But there are some mistake because i use 241 weights.
>>
>> if anyone can tell me how to filter out it in using NCL
>> build-in function or has any other method for the band-pass filter
>>
>> ------------------------------------------------------------------------
>

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

begin
;----------------------------------------------------------
; Specify (1) start and end dates [yyyymmdd]
; (2) region
;----------------------------------------------------------

     dStrt = 20010101 ; start date (yyyymmdd)
     dLast = 20061231 ; last date

     latS = -30
     latN = 30
     lonL = 0
     lonR = 360

;----------------------------------------------------------
; generate Lancos weights
;----------------------------------------------------------
     nwgt = 241 ; number of weights used by Matthews
     fca = 1./200. ; 200 day cutoff
     fcb = 1./20. ; 20 day cutoff
     ihp = 2 ; bandpass
     sigma = 1.
     wgts = filwgts_lancos (nwgt, ihp, fca, fcb, sigma)

;----------------------------------------------------------
; Open file
; Assuming you are looking at a specific period
; read the time [units days since...] and convert to yyyymmdd
; Create indicies for the selected period.
; If data are of type short convert to float.
;----------------------------------------------------------

     f = addfile ("olr.day.mean.nc", "r")

     time = f->time
     date = ut_calendar(time, -2) ; date: yyyymmdd
     delete(time)

     iStrt = ind(date.eq.dStrt) ; index of start date/time
     iLast = ind(date.eq.dLast) ; last
     delete(date)
     
     OLR = short2flt( f->olr(iStrt:iLast,{latS:latN},{lonL:lonR}) ) ; (time,lat,lon)
     printVarSummary( OLR )
     printMinMax( OLR, True )

;----------------------------------------------------------
; OLR must be reordered for next set of operations
;----------------------------------------------------------
     olr = OLR(lat|:,lon|:,time|:)
     delete(OLR)

;----------------------------------------------------------
; As per Matthews paper:
; Remove first 12 harmonics: taper 1st to minimize leakage
;----------------------------------------------------------
     olr = taper( olr, 0.1, 0)

     cf = ezfftf( olr ) ; (2,nlat,mlon,ntim/2)
     printVarSummary( cf )
     cf(:,:,:,0:12) = 0.0 ; set fourier coef to 0.0
     olr = ezfftb( cf , 0.0 ) ; reconstruct
     printVarSummary( olr )
     printMinMax( olr, True )

;----------------------------------------------------------
; Apply Lancos weights to the reconstructed data
;----------------------------------------------------------
     olr = wgt_runave(olr, wgts, 0)
     printVarSummary( olr )
     printMinMax( olr, True )
     
;----------------------------------------------------------
; For convenience ... save to netCDF
; Let NCL do the tedious work. This may be slow if you
; have many variables or large records.
;----------------------------------------------------------
; Eliminate attributes associated with the original data
; that are no longer applicable
;----------------------------------------------------------
     delete(olr_at_precision)
     delete(olr_at_valid_range)
     delete(olr_at_actual_range)
     olr_at_statistic = "tapered; 1st 12 harmonics removed; 241 Lancos filter applied"
     olr@_FillValue = 1.e20 ; change for convenience
     olr_at_missing_value = olr@_FillValue

;----------------------------------------------------------
; create some global attributes
;----------------------------------------------------------
     nline = inttochar(10) ; new line character

     globeAtt = 1
     globeAtt_at_title = "OLR: Band Passed with 241-wgt Lancos Filter: 20-200 day"
     globeAtt@data_source = "http://www.cdc.noaa.gov/ ; Search PSD Gridded Climate Data"
     globeAtt_at_reference = nline + \
          "ADRIAN J. MATTHEWS "+nline+\
          "Intraseasonal Variability over Tropical Africa during Northern Summer"+nline+\
          "Journal of Climate, 2004: 17, 2427-2440"+nline
     globeAtt_at_creation_date= systemfunc ("date" )

;----------------------------------------------------------
; create a date variable for convenience
;----------------------------------------------------------

     date = ut_calendar(olr&time, -2) ; date: yyyymmdd
     date!0= "time"
     date_at_units = "yyyymmdd"

;----------------------------------------------------------
; create netCDF
;----------------------------------------------------------

     diro = "./" ; output directory
     filo = "OLR.matthews.nc"
     NCFILE = diro+filo
     system ("/bin/rm -f " + NCFILE) ; remove any pre-exist file

     ncdf = addfile(NCFILE,"c")
     fileattdef( ncdf, globeAtt ) ; assign globe file attributes
     ncdf->date = date
     ncdf->OLR = olr ; write data (lat,lon,time)
end

_______________________________________________
ncl-talk mailing list
ncl-talk_at_ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Wed May 16 2007 - 13:09:01 MDT