[ncl-talk] array indexing in ncl
Dennis Shea
shea at ucar.edu
Tue Nov 28 21:05:27 MST 2017
As noted ... NCL syntax operates differently. C'est la vie!
Two elaborations of DaveA's response.
---
T = 100
N = 20
M = 30
x = random_normal(0,10, (/T,N,M/))
n = (/1,2/)
m = (/2,9/)
kpts = dimsizes(n)
z = new( (/T,kpts/), typeof(x))
do k=0,kpts-1
z(:,k) = x(:,n(k),m(k))
end do
printVarSummary(z)
===============================================
undef("extract_loc")
function extract_loc(x:numeric, n:integer, m:integer)
local k, kpts, z, dimx, T
begin
dimx = dimsizes(x)
T = dimx(0)
kpts = dimsizes(n)
z = new( (/T,kpts/), typeof(x))
do k=0,kpts-1
z(:,k) = x(:,n(k),m(k))
end do
return(z)
end
;-----------------------------------
; MAIN
;-----------------------------------
T = 100
N = 20
M = 30
x = random_normal(0,10, (/T,N,M/))
n = (/1,2/)
m = (/2,9/)
Z = extract_loc(x,n,m) ; Z = extract_loc(x,(/1,2/),(/2,9/))
printVarSummary(Z)
------------------------------------------------------------------------------
D
On Tue, Nov 28, 2017 at 8:44 PM, Dave Allured - NOAA Affiliate <
dave.allured at noaa.gov> wrote:
> Jatin,
>
> In NCL, you can do this with only one loop over your two-column list of
> point indices. I think this is reasonably obvious and efficient. Here is
> a compact example to match your example, untested:
>
> pre-define hfx2 as (103,2)
> hfx2(:,0) = hfx(:,1,2)
> hfx2(:,1) = hfx(:,2,4)
>
> NCL also provides coordinate subscripting, which may add some value when
> picking out grid points like this. Then you don't need to separately
> convert coordinates to indices. Note the curly brackets:
>
> hfx2(:,0) = hfx(: , {lats(0)} , {lons(0)} )
> etc.
>
> My editorial opinion is that is a rather fancy multidimensional indexing
> scheme that python provides, perhaps not altogether intuitive. Well, um,
> NCL's curly brackets are also a bit unusual, but you get used to them.
>
> --Dave
>
>
> On Tue, Nov 28, 2017 at 7:59 PM, Jatin Kala <jatin.kala.jk at gmail.com>
> wrote:
>
>> Hi,
>>
>> Suppose we have a 3-D array, size (103,lat,lon), or whatever...
>>
>> If i want to extract a time series at two points, say (:,1,2) and
>> (:,2,4), then in Python you can do it quite simply, and you end up with an
>> array which is size (103,2), which one would expect:
>>
>> import netCDF4 as nc
>> f = nc.Dataset('wrfout_d03_cntl')
>> hfx = f.variables['HFX'][:]
>> hfx2 = hfx[:,[1,2],[2,4]]
>> print(hfx2.shape)
>>
>> But if i do the same thing in ncl, i actually get an array (103,2,2),
>> which is not what one would expect:
>>
>> f = addfile("wrfout_d03_cntl","r")
>> hfx = f->HFX
>> hfx2 = hfx(:,(/1,2/),(/2,4/))
>> print(dimsizes(hfx2))
>>
>> The only way to get the right answer in ncl, is to pre-define hfx as
>> (103,2), then loop through the indices, and put in the array within the
>> loop.
>>
>> Does anyone else find this rather counter-intuitive? I.e., the way ncl
>> does array indexing in such cases? Is there a way to achieve the same
>> result in ncl without having to loop?
>>
>> Cheers,
>>
>> Jatin
>>
>
> _______________________________________________
> 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/20171128/abc85717/attachment.html>
More information about the ncl-talk
mailing list