[ncl-talk] Loop Over Latitude and Level Arrays
Dennis Shea
shea at ucar.edu
Tue Dec 31 13:43:20 MST 2019
[0]
ncl-talk would like to help but the *kf_filter.ncl* function and NCL driver
script are *user contributed *codes.
They were not developed by NCL core developers. Hence, we are not familiar
with the code.
[1]
You *must *be much clearer in your description of your input variable.
You state you have an array ordered:
*(latitude, longitude, level, time) *
Is this fortran [column major] ordering?
Typically*, *a 4D variable from a netCDF file has row-major ordering
*(time,level,lat,lon)*
What is the output from*:*
data = fin->$varIn$
printVarSummary(data)
[2]
As indicated in the function documentation the input array is* prototyped *as
2-dimensional: *inData[*][*]*
This means it will only allow two dimensional ( [*][*] )input.
;******************************************************************************
; kf_filter.ncl
; Carl Schreck (carl at atmos.albany.edu)
; January 2009
;******************************************************************************
undef("kf_filter")
function kf_filter( *inData[*][*]:*float, obsPerDay:numeric, \\
tMin:numeric, tMax:numeric, kMin:numeric, kMax:numeric,
\\
hMin:numeric, hMax:numeric, waveName:string )
;******************************************************************************
[snip]
; Input Variables:
*; inData: the data to be filtered. time should be 1st coordinate;
lon should be 2nd coordinate.*
;[snip]
; Return Value:
; * retVal[*][*]: the filtered data*
;******************************************************************************
=======================================================
[3]
In addition, the function documentation further states that the
two-dimensional input variable *must* be ordered: *(time,lon)*.
Any other 2D-ordering, will result in erroneous results.
[4]
If the data variable has row-major ordering
dim_data = *dimsizes*
<http://www.ncl.ucar.edu/Document/Functions/Built-in/dimsizes.shtml>(data)
print(dim_data)
ntim = dim_data(0)
klev = dim_data(1)
nlat = dim_data(2)
mlon = dim_data(3)
To loop
KF_SAVE = *new*
<http://www.ncl.ucar.edu/Document/Functions/Built-in/new.shtml>((/ntim,klev,nlat,mlon/),
"float")
do kl=0,klev-1
do nl=0,nlat-1
KF_SAVE( :,kl,nl,:) = kf_filter(data(:,kl,nl,:), ....) *;
data(:,kl,nl,:) reduces to 2D: (time,lon)*
end do ; nl
end do ; kl
KF_SAVE at long_name = "..."
copy_VarCoords(data,KF_SAVE)
printVarSummary(KF_SAVE)
printMinMax(KF_SAVE,0)
If your data are column major, you will have to make the appropriate
changes.
Good Luck
On Tue, Dec 31, 2019 at 5:56 AM Barry Lynn via ncl-talk <ncl-talk at ucar.edu>
wrote:
> Hi:
>
> You need to define your array as follows:
>
> Two_D_Array = Four_D_Array(i_lat,:,i_level,:) ; i_lat is the latitude
> array index and i_level is the level index
>
> On Tue, Dec 31, 2019 at 1:27 PM ANDIKA FAUZIAH HAPSARI, A.P <
> andika.hapsari at bmkg.go.id> wrote:
>
>> *Oh I'm so sorry, I clicked the reply button instead of repy all. Thank
>> you for remind me.*
>>
>> I think it's becasue kf_filter.ncl script only works on two dimension,
>> but when I tried to modify kf_filter, it couldn't work (but i think i did
>> it wrong too)
>>
>> I have tried 3d array, and the error message I got was
>> fatal:Number of subscripts do not match number of dimensions of
>> variable,(3) Subscripts used, (4) Subscripts expected
>> fatal:["Execute.c":8637]:Execute: Error occurred at or near line 229 in
>> file filter_vertical_waves_wrf2.ncl
>>
>> So I don't understand which part I got it wrong.
>>
>>
>> From: "Barry Lynn (barry.h.lynn at gmail.com)" <barry.h.lynn at gmail.com>
>> To: "ANDIKA FAUZIAH HAPSARI, A.P" <andika.hapsari at bmkg.go.id>, ncl-talk <
>> ncl-talk at ucar.edu>
>> Date: Tue, 31 Dec 2019 09:10:02 +0200
>> Subject: Re: [ncl-talk] Loop Over Latitude and Level Arrays
>>
>> Hi:
>>
>> Please also respond to the list. Someone on the list might be able to
>> help you faster than I can. Also, members of the list can learn useful
>> information from this exchange.
>>
>> It would be helpful if you would verify where the script fails, whether
>> it is the one you are using or a called routine.
>>
>> Lastly, the script is telling you that your variable is missing
>> information about what it is. You need to name it, as shown in the NCL
>> help pages.
>>
>> But, before you do that, try reading your 4d array into a 2d array and
>> see if you get the same error using the 2d array.
>>
>> Barry
>>
>> On Tue, Dec 31, 2019 at 8:59 AM ANDIKA FAUZIAH HAPSARI, A.P <
>> andika.hapsari at bmkg.go.id> wrote:
>>
>>>
>>> Thanks for your answer. But I'm not sure if I can pass the 4D variable
>>> directly. Because when I used this script below, I got this error message :
>>>
>>> warning:Argument 0 of the current function or procedure was coerced to
>>> the appropriate type and thus will not change if the function or procedure
>>> modifies its value
>>> fatal:["NclAtt.c":262]:Attribute assignment type mismatch.
>>> fatal:(lon) is not a named dimension in variable (inData).
>>> fatal:["Execute.c":8637]:Execute: Error occurred at or near line 36 in
>>> file $NCARG_ROOT/lib/ncarg/nclscripts/contrib/kf_filter.ncl
>>> fatal:["Execute.c":8637]:Execute: Error occurred at or near line 234 in
>>> file filter_vertical_waves_wrf.ncl
>>>
>>>
>>> Or maybe, is there something wrong with my script?
>>> Thanks in advance
>>>
>>>
>>> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>>> ; filter_waves.ncl
>>> ; Carl Schreck (carl at cicsnc.org)
>>> ; February 2011
>>> ; Updated for style October 2011
>>>
>>> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>>> ; Description: Filter data for each equatorial wave type
>>>
>>> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>>> 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/csm/shea_util.ncl"
>>> ;load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
>>> ;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/diagnostics_cam.ncl"
>>> load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/time_axis_labels.ncl"
>>> load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/kf_filter.ncl"
>>> load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/cd_string.ncl"
>>> load "$NCARG_ROOT/lib/ncarg/nclscripts/print_clock.ncl"
>>> ;load "~/lib/print_clock.ncl"
>>>
>>> begin
>>>
>>> print_clock( "Here we go! " )
>>> ; These are some parameters that could be useful to have up top
>>> mis = -999
>>> mis at _FillValue = -999
>>> obsPerDay = 4
>>> if( .not.isvar("varIn") ) then
>>> varIn = "u"
>>> end if
>>> if( .not.isvar("latitudeBand") ) then
>>> latitudeBand = "15N-15N"
>>> end if
>>> print( varIn + "." + latitudeBand )
>>> basePath = "/home/sarihmmm/datareanly/reanalysis/"
>>> pathIn = basePath + "anginU.nc"
>>> pathOut = basePath + "anginU" + varIn + "." + latitudeBand + ".
>>> vertkelvin.nc"
>>> wk99 = False ; True uses the original Wheeler-Kiladis filters,
>>> ; False uses slight modifications from recent papers
>>> calcHigh = False
>>> calcLow = False
>>> calcMrg = False
>>> calcMjo = False
>>> calcKelvin = True
>>> calcEr = False
>>> calcMtd = False
>>> calcTd = False
>>> makeNewFile = True
>>> ; Open the input files
>>> fin = addfile( pathIn, "r" )
>>> time = fin->time
>>> level = fin->level
>>> longitude = fin->longitude
>>> latitude = fin->latitude
>>> ; Open the output files
>>> setfileoption( "nc", "Format", "LargeFile" )
>>> if( makeNewFile ) then
>>> system( "rm " + pathOut )
>>> fout = addfile( pathOut, "c" )
>>> fout->time = time
>>> fout->level = level
>>> fout->longitude = longitude
>>> fout->latitude = latitude
>>> else
>>> fout = addfile( pathOut, "w" )
>>> end if
>>> print_clock( "Reading the input data... " )
>>> data = fin->$varIn$
>>> data&time at beginning_date = cd_string( data&time(0), "" )
>>> data&time at ending_date = cd_string(
>>> data&time(dimsizes(data&time)-1), "" )
>>> if( any( ismissing(data) ) ) then
>>> print_clock( "WARNING: Setting missing data to zero" )
>>> data = where( ismissing(data), 0, data )
>>> end if
>>> printVarSummary(data)
>>> if( calcHigh ) then
>>> print_clock( "Filtering High..." )
>>> high = data
>>> high at longitudeg_name = "High frequency"
>>> high at filter = "Highpass time filter"
>>> high at wavenumber = (/ -9999, 9999 /)
>>> high at period = (/ 0.01, 120 /)
>>> high at depth = (/ mis, mis /)
>>> do z = 0, ( dimsizes(level) - 1 )
>>> ; print_clock( z + " " + ( dimsizes(level) - 1 ) + " " )
>>> high(time|:,{level|level(z)},longitude|:) = (/ kf_filter( \\
>>> data(time|:,{level|level(z)},longitude|:), obsPerDay, \\
>>> high at period(0), high at period(1), \\
>>> high at wavenumber(0), high at wavenumber(1), \\
>>> high at depth(0), high at depth(1), \\
>>> "none" ) /)
>>> end do
>>> print_clock( "Writing high... " )
>>> fout->high = high
>>> delete(high)
>>> end if
>>>
>>> if( calcLow ) then
>>> print_clock( "Filtering Low..." )
>>> low = data
>>> low at longitudeg_name = "Low frequency"
>>> low at filter = "Lowpass time filter"
>>> low at wavenumber = (/ -9999, 9999 /)
>>> low at period = (/ 120, 9999 /)
>>> low at depth = (/ mis, mis /)
>>> do z = 0, ( dimsizes(level) - 1 )
>>> ; print_clock( z + " " + ( dimsizes(level) - 1 ) + " " )
>>> low(time|:,{level|level(z)},longitude|:,latitude|:) = (/
>>> kf_filter( \\
>>> data(time|:,{level|level(z)},longitude|:,latitude|:),
>>> obsPerDay, \\
>>> low at period(0), low at period(1), \\
>>> low at wavenumber(0), low at wavenumber(1), \\
>>> low at depth(0), low at depth(1), \\
>>> "none" ) /)
>>> end do
>>> print_clock( "Writing Low... " )
>>> fout->low = low
>>> delete(low)
>>> end if
>>>
>>> if( calcMrg ) then
>>> print_clock( "Filtering MRG..." )
>>> mrg = data
>>> mrg at longitudeg_name = "Mixed Rossby-Gravity Waves in " +
>>> str_upper(varIn)
>>> mrg at filter = "Wheeler & Kiladis (1999)"
>>> mrg at wavenumber = (/ -10, -1 /)
>>> mrg at period = (/ 3, 96 /)
>>> mrg at depth = (/ 8, 90 /)
>>> do z = 0, ( dimsizes(level) - 1 )
>>> ; print_clock( z + " " + ( dimsizes(level) - 1 ) + " " )
>>> mrg(time|:,{level|level(z)},longitude|:) = (/ kf_filter( \\
>>> data(time|:,{level|level(z)},longitude|:), obsPerDay, \\
>>> mrg at period(0), mrg at period(1), \\
>>> mrg at wavenumber(0), mrg at wavenumber(1), \\
>>> mrg at depth(0), mrg at depth(1), \\
>>> "mrg" ) /)
>>> end do
>>> print_clock( "Writing MRG... " )
>>> fout->mrg = mrg
>>> delete(mrg)
>>> end if
>>>
>>> if( calcMjo ) then
>>> print_clock( "Filtering MJO..." )
>>> mjo = data
>>> mjo at longitudeg_name = "Madden-Julian Oscillatitudeion in " +
>>> str_upper(varIn)
>>> if( wk99 ) then
>>> mjo at filter = "Wheeler & Kiladis (1999)"
>>> mjo at wavenumber = (/ 1, 5 /)
>>> ; mjo at wavenumber = (/ 0, 9 /)
>>> mjo at period = (/ 30, 96 /)
>>> mjo at depth = (/ mis, mis /)
>>> else
>>> mjo at filter = "Kiladis et al. (2005 JAS)"
>>> mjo at wavenumber = (/ 0, 9 /)
>>> mjo at period = (/ 20, 100 /)
>>> mjo at depth = (/ mis, mis /)
>>> end if
>>> do z = 0, ( dimsizes(level) - 1 )
>>> ; print_clock( z + " " + ( dimsizes(level) - 1 ) + " " )
>>> mjo(time|:,{level|level(z)},longitude|:) = (/ kf_filter( \\
>>> data(time|:,{level|level(z)},longitude|:), obsPerDay, \\
>>> mjo at period(0), mjo at period(1), \\
>>> mjo at wavenumber(0), mjo at wavenumber(1), \\
>>> mjo at depth(0), mjo at depth(1), \\
>>> "none" ) /)
>>> end do
>>> print_clock( "Writing MJO... " )
>>> fout->mjo = mjo
>>> delete(mjo)
>>> end if
>>>
>>> if( calcKelvin ) then
>>> print_clock( "Filtering KELVIN..." )
>>> kelvin = data
>>> kelvin at longitudeg_name = "Kelvin Waves in " + str_upper(varIn)
>>> if( wk99 ) then
>>> kelvin at filter = "Wheeler & Kiladis (1999)"
>>> kelvin at wavenumber = (/ 1, 14 /)
>>> kelvin at period = (/ 2.5, 30 /)
>>> kelvin at depth = (/ 8, 90 /)
>>> else
>>> kelvin at filter = "Straub & Kiladis (2002)"
>>> kelvin at wavenumber = (/ 1, 14 /)
>>> kelvin at period = (/ 2.5, 17 /)
>>> kelvin at depth = (/ 5, 90 /)
>>> end if
>>> do z = 0, ( dimsizes(level) - 1 )
>>> do y = 0, ( dimsizes(level) - 1 )
>>> ; print_clock( z + " " + ( dimsizes(level) - 1 ) + " " )
>>> kelvin(time|:,{level|level(z)},{latitude|latitude(y)},longitude|:)
>>> = (/ kf_filter( \\
>>>
>>> data(time|:,{level|level(z)},{latitude|latitude(y)},longitude|:),
>>> obsPerDay, \\
>>> kelvin at period(0), kelvin at period(1), \\
>>> kelvin at wavenumber(0), kelvin at wavenumber(1), \\
>>> kelvin at depth(0), kelvin at depth(1), \\
>>> "kelvin" ) /)
>>> end do
>>> end do
>>> print_clock( "Writing KELVIN... " )
>>> fout->kelvin = kelvin
>>> delete(kelvin)
>>> end if
>>>
>>> if( calcEr ) then
>>> print_clock( "Filtering ER..." )
>>> er = data
>>> er at longitudeg_name = "Equatorial Rossby Waves in " +
>>> str_upper(varIn)
>>> if( wk99 ) then
>>> er at filter = "Wheeler & Kiladis (1999)"
>>> er at wavenumber = (/ -10, -1 /)
>>> er at period = (/ 9.7, 48 /)
>>> er at depth = (/ 8, 90 /)
>>> else
>>> er at filter = "Kiladis et al. (2009 Rev. Geophys.)"
>>> er at wavenumber = (/ -10, -1 /)
>>> er at period = (/ 9.7, 48 /)
>>> er at depth = (/ 5, 90 /)
>>> end if
>>> do z = 0, ( dimsizes(level) - 1 )
>>> ; print_clock( z + " " + ( dimsizes(level) - 1 ) + " " )
>>> er(time|:,{level|level(z)},longitude|:) = (/ kf_filter( \\
>>> data(time|:,{level|level(z)},longitude|:), obsPerDay, \\
>>> er at period(0), er at period(1), \\
>>> er at wavenumber(0), er at wavenumber(1), \\
>>> er at depth(0), er at depth(1), \\
>>> "er" ) /)
>>> end do
>>> print_clock( "Writing ER... " )
>>> fout->er = er
>>> delete(er)
>>> end if
>>> if( calcMtd ) then
>>> print_clock( "Filtering MRG/TD..." )
>>> td = data
>>> td at longitudeg_name = "MRG/TD in " + str_upper(varIn)
>>> td at filter = "Frank and Roundy (2006), extended to
>>> wavenumber -20"
>>> td at wavenumber = (/ -20, -0 /)
>>> td at period = (/ 2.5, 10 /)
>>> td at depth = (/ mis, mis /)
>>> do z = 0, ( dimsizes(level) - 1 )
>>> ; print_clock( z + " " + ( dimsizes(level) - 1 ) + " " )
>>> td(time|:,{level|level(z)},longitude|:) = (/ kf_filter( \\
>>> data(time|:,{level|level(z)},longitude|:), obsPerDay, \\
>>> td at period(0), td at period(1), \\
>>> td at wavenumber(0), td at wavenumber(1), \\
>>> td at depth(0), td at depth(1), \\
>>> "none" ) /)
>>> end do
>>> print_clock( "Writing TD... " )
>>> fout->td = td
>>> delete(td)
>>> end if
>>>
>>> if( calcTd ) then
>>> print_clock( "Filtering TD..." )
>>> td = data
>>> td at longitudeg_name = "TD-type Disturbances in " +
>>> str_upper(varIn)
>>> td at filter = "Roundy and Frank (2004), extended to
>>> wavenumber -20"
>>> td at wavenumber = (/ -20, -6 /)
>>> td at period = (/ 2.5, 5 /)
>>> td at depth = (/ 90, mis /)
>>> do z = 0, ( dimsizes(level) - 1 )
>>> ; print_clock( z + " " + ( dimsizes(level) - 1 ) + " " )
>>> td(time|:,{level|level(z)},longitude|:) = (/ kf_filter( \\
>>> data(time|:,{level|level(z)},longitude|:), obsPerDay, \\
>>> td at period(0), td at period(1), \\
>>> td at wavenumber(0), td at wavenumber(1), \\
>>> td at depth(0), td at depth(1), \\
>>> "none" ) /)
>>> end do
>>> print_clock( "Writing TD... " )
>>> fout->td = td
>>> delete(td)
>>> end if
>>>
>>> print_clock( "Closing file..." )
>>> delete(fout)
>>> print_clock( "Thank you, come again. " )
>>> end ; filter_waves
>>>
>>>
>>> From: "Barry Lynn (barry.h.lynn at gmail.com)" <barry.h.lynn at gmail.com>
>>> To: "ANDIKA FAUZIAH HAPSARI, A.P" <andika.hapsari at bmkg.go.id>
>>> Cc: "ncl-talk at ucar.edu" <ncl-talk at ucar.edu>
>>> Date: Tue, 31 Dec 2019 07:48:15 +0200
>>> Subject: Re: [ncl-talk] Loop Over Latitude and Level Arrays
>>>
>>> Hi:
>>>
>>> When you pass the variable to the kf_filter, you can loop over both over
>>> latitude and level, as you pass the 4D variable directly into the routine.
>>>
>>> I am not sure if there is a faster way,
>>>
>>> Barry
>>>
>>> On Mon, Dec 30, 2019 at 8:17 PM ANDIKA FAUZIAH HAPSARI, A.P via ncl-talk
>>> <ncl-talk at ucar.edu> wrote:
>>>
>>>> Hello..
>>>>
>>>> I have 4 Dimension data (latitude, longitude, level, time) and I want
>>>> to filter kelvin waves vertically which use kf_filter function that only
>>>> has time and longitude. And I need to loop over both latitude and level so
>>>> that the array can pass kf_filter.
>>>>
>>>> Is there any instruction and script to do loop over latitude and level
>>>> arrays?
>>>>
>>>> Thanks in advance.
>>>> _______________________________________________
>>>> ncl-talk mailing list
>>>> ncl-talk at ucar.edu
>>>> List instructions, subscriber options, unsubscribe:
>>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>>
>>>
>>>
>>> --
>>> Barry H. Lynn, Ph.D
>>> Senior Associate Scientist, Lecturer,
>>> The Institute of the Earth Science,
>>> The Hebrew University of Jerusalem,
>>> Givat Ram, Jerusalem 91904, Israel
>>> Tel: 972 547 231 170
>>> Fax: (972)-25662581
>>>
>>> C.E.O, Weather It Is, LTD
>>> Weather and Climate Focus
>>> http://weather-it-is.com
>>> Jerusalem, Israel
>>> Local: 02 930 9525
>>> Cell: 054 7 231 170
>>> Int-IS: x972 2 930 9525
>>>
>>>
>>>
>>
>> --
>> Barry H. Lynn, Ph.D
>> Senior Associate Scientist, Lecturer,
>> The Institute of the Earth Science,
>> The Hebrew University of Jerusalem,
>> Givat Ram, Jerusalem 91904, Israel
>> Tel: 972 547 231 170
>> Fax: (972)-25662581
>>
>> C.E.O, Weather It Is, LTD
>> Weather and Climate Focus
>> http://weather-it-is.com
>> Jerusalem, Israel
>> Local: 02 930 9525
>> Cell: 054 7 231 170
>> Int-IS: x972 2 930 9525
>>
>>
>>
>
> --
> Barry H. Lynn, Ph.D
> Senior Associate Scientist, Lecturer,
> The Institute of the Earth Science,
> The Hebrew University of Jerusalem,
> Givat Ram, Jerusalem 91904, Israel
> Tel: 972 547 231 170
> Fax: (972)-25662581
>
> C.E.O, Weather It Is, LTD
> Weather and Climate Focus
> http://weather-it-is.com
> Jerusalem, Israel
> Local: 02 930 9525
> Cell: 054 7 231 170
> Int-IS: x972 2 930 9525
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20191231/67224c66/attachment.html>
More information about the ncl-talk
mailing list