[ncl-talk] array indexing in ncl

[Dennis Shea]

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>