[ncl-talk] Advice on extracting point locations from array

[Will Hobbs]

Hi Dennis

Thanks for the reply; I’m already saving the profiles to file for later work, but I need to do this for 3 model experiments and ~15 CMIP5 models. The script is slow but low on memory and CPU burden, so maybe I just need to dip my toe into parallel processing…

>I'm not sure how the underlying netCDF software would implement something like


>Tprof = fi->T(pInd(:,0),:,pInd(:,1),pInd(:,2))


For future reference, it returns a HUGE array, because each time index will have every profile’s lat and every profile’s lon, so the array is (nprofile x nlev nprofile x nprofile). With a small number of profiles (i.e. less than the time, lat, and lon dimensions) it would be fine, but I have ~50,000.

Will

From: Dennis Shea <shea at ucar.edu>
Date: Wednesday, 17 October 2018 at 12:45 PM
To: Will Hobbs <will.hobbs at utas.edu.au>
Cc: Ncl-talk <ncl-talk at ucar.edu>
Subject: Re: [ncl-talk] Advice on extracting point locations from array

'Brain Trust' ... "great idea"  ===> Nope ... Sorry
---
IMHO:
I don't think the loop can be replaced by array operations. I'm not sure how the underlying netCDF software would implement something like

Tprof = fi->T(pInd(:,0),:,pInd(:,1),pInd(:,2))
---

  nprofile   = ...
  nlev        = ...
  Tprof      = new((/nprofile,nlev/),float)   ; reorder from your code
  begTime = get_cpu_time()
  do i = 0, nprofile-1
     Tprof(i,:) = fi->T(pInd(i,0),:,pInd(i,1),pInd(i,2))    ;pInd is an array with coord triplets
  end do

  printVarSummary(Tprof)
  printMinMax(Tprof,0)
  print("Tprof generation time: " + (get_cpu_time() - begTime)+" seconds; nprofile="+nprofile)
  print("---")
---

If it takes a long time, I'd suggest saving the created variable(s) in a netCDF file. Untested:

;------------------
;  netCDF
;-----------------
  netCDF = True
  if (netCDF) then
      profile   = ispan(0,nprofile-1,1)   ; or, ispan(1,nprofile,1)
      profile at long_name = "Profile Number"
      profile!0 = "profile"
      profile&profile = profile

      TIME   = fi->time(pInd(:,0))
      LAT     = fi->lat(pInd(:,1))
      LON    = fi->lon(pInd(:,2))

      Tprof!0= "profile"
      TIME!0 = "profile"
      LAT!0   = "profile"
      LON!0  = "profile"
      Tprof&profile  = profile
      TIME&profile  = profile
      LAT&profile    = profile
      LON&profile   = profile

      dir_nc  = "./"
      fil_nc   = "FOO.nc<http://FOO.nc>"
      pth_nc = dir_nc + fil_nc
      print(pth_nc)

      system("/bin/rm -f "+pth_nc)           ; remove any pre-existing file
      ncdf = addfile(pth_nc ,"c")            ; open output netCDF file

;===================================================================
; Create global attributes of the file (optional)
;===================================================================
       fAtt               = True            ; assign file attributes
       fAtt at title         = "Selected Profiles"
      ;fAtt at source_info   =  "..."
       fAtt at Conventions   = "None"
       fAtt at creation_date = systemfunc ("date")
       fileattdef( ncdf, fAtt )            ; copy file attributes

;===================================================================
; Make profile an UNLIMITED dimension; recommended  for most applications
; Allows for future profiles to be appended
;===================================================================
       filedimdeff(ncdf,"profile",-1,True)

;===================================================================
; output variable(s) directly; NCL will call appropriate functions
;===================================================================
       ncdf->TIME  = TIME
       ncdf->LAT   = LAT
       ncdf->LON   = LON
       ncdf->TPROF = Tprof
  end if   ; netCDF



On Mon, Oct 15, 2018 at 11:53 PM Will Hobbs <will.hobbs at utas.edu.au<mailto:will.hobbs at utas.edu.au>> wrote:
Hi all

This one is for the NCL ‘Brains Trust’. I’m extracting a large number of ocean profiles from a 4-d array (dimensioned time,depth,lat,lon), using coordinate triplets based on time,lat and lon, (i.e. a 2-d array of indices dimensioned nprofile, 3),  but it’s a slow process and I wonder if anyone has any smart ideas on speeding it up.

At the moment (partly to manage memory), I’m reading each profile in a loop, viz:


Tprof = new((/nlev,nprofile/),float)
do i = 0, nprofile-1
     Tprof(:,i) = fi->T(pInd(i,0),:,pInd(i,1),pInd(i,2))    ;pInd is an array with coord triplets
end do

Obviously one obvious way of speeding things up would be to read the entire array of input data and extract the locations, to avoid multiple i/o calls. What I really want to do though is somehow get rid of the loop entirely. Is there any way of extracting coordinate pairs (or in this case triplets) without looping through each profile?

The only way I can think of is to turn the input array (T in my example above) into a 1-d array, and somehow convert ‘pInd’ into the elements of that 1-d array – this seems fraught with opportunity for error though. Catastrophic error I can deal with – it’s the sneaky, inobvious error that scares me….

Hoping someone has a great idea….

Will


University of Tasmania Electronic Communications Policy (December, 2014).
This email is confidential, and is for the intended recipient only. Access, disclosure, copying, distribution, or reliance on any of it by anyone outside the intended recipient organisation is prohibited and may be a criminal offence. Please delete if obtained in error and email confirmation to the sender. The views expressed in this email are not necessarily the views of the University of Tasmania, unless clearly intended otherwise.
_______________________________________________
ncl-talk mailing list
ncl-talk at ucar.edu<mailto:ncl-talk at ucar.edu>
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk<http://mailman.ucar.edu/mailman/listinfo/ncl-talk>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20181017/6b5251bb/attachment.html>