Re: area averaging across 0E.

Cathy,

Here is a general solution. This method allows the requested bounds
to be independent from the actual range of coordinates. For example,
data set coordinates from 0 to 360, requested bounds -120 to -60.
Split longitudes are also handled correctly, e.g -20 to +20.

mlons = modulo (lons(:) - minlon, 360) + minlon
                ; lons modular within requested range
xind = ind (minlon .le. mlons .and. mlons .le. maxlon)
                      ; get indices of lons within range; this
                   ; method takes care of split lon subsets
nxind = num (.not. ismissing (xind))

; Make 3-D subset array for current region.
; Special handling needed to prevent dimension reduction.

if (nxind .eq. 1) then
   subset = orig(:, ilat1:ilat2, xind:xind) ; single longitude
else
   subset = orig(:, ilat1:ilat2, xind) ; two or more, vector subscripting
end if

Note the custom "modulo" math function. For this formula, It is
important to use modulo rather than NCL's mod function or % (mod)
operator. Otherwise you can get incorrect results when negative
coordinates are involved.

http://en.wikipedia.org/wiki/Modulo_operation

You can download modulo.ncl here:
http://www.esrl.noaa.gov/psd/people/dave.allured/data/ncl/lib/modulo.ncl

You might want to test for nxind = 0. This means that the specified
range of longitudes is so narrow that it is in between adjacent
longitudes in the data set. Special handling or an error message is
needed.

This method is inefficient for very large file variables, because of
the vector subscripting. Either read the whole array into memory
before subsetting, or modify the method to read the array in two parts
with range subscripting.

--Dave

On Fri, Nov 9, 2012 at 11:58 AM, Cathy Smith (NOAA Affiliate)
<cathy.smith@noaa.gov> wrote:
> All
>
> I have gridded data that generally goes from 0 to 360 in longitude
> (sometimes -180 to 180 though I can easily switch that to 0 to 360). I
> would like to produce area averages from this data.
> Sometimes, the area to be averaged will cross 0E (350E to 10E, say). I
> looked and could not find any routine that does this. I thought
> wgt_areaave_Wrap did but it doesn't seem to (and it doesn't give an error).
> How can I easily do this? A few ways that occurred to me are to append
> arrays in longitude (0 to 720) which will create very big arrays as I'm
> looking at time series. Another alternative is to split up
> the calculation at 0 but I would need to combine the results back
> together correctly. I could flip the longitude depending on the desired
> area but that wouldn't handle -10 to 200, for example though it would
> handle -10 to 120.
>
> I looked through examples and didn't see one though I would think it
> would be a common task. So, maybe I'm just misunderstanding something.
>
> Cathy Smith
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Fri Nov 9 13:37:28 2012