CDO

Hi Andreas,

in my opinion, there is no way to retrieve the indices of given lat/lon values.

But there is a workaround using NCL. Here, a ksh script is used to create the NCL script, commit the file name,lat and lon values to the NCL script which returns the indices:

#!/bin/ksh

  infile="infile.nc" 

  set -A lats
  set -A lons

  lats+=( 31.0 17.31 24.05 )
  lons+=( -86.45 -101.00 -92.46 )

  cat <<EOF > ncl_script
begin

  if(.not. isvar("infile")) then
     print("ERROR: infile not given")
     exit
  end if
  if(.not. isvar("lat_array")) then
     print("ERROR: lat_array variable not set")
     exit
  else
     lat_array := tofloat(str_split(lat_array," "))
  end if
  if(.not. isvar("lon_array")) then
     print("ERROR: lon_array variable not set")
     exit
  else
     lon_array := tofloat(str_split(lon_array," "))
  end if

  f     =  addfile(infile,"r")                  ;-- open file
  lat2d =  f->lat                               ;-- get 2D latitude array
  lon2d =  f->lon                               ;-- get 2D longitude array

  ndims   = dimsizes(lat_array)
  indices = getind_latlon2d (lat2d,lon2d, lat_array, lon_array)     ;-- retrieve the indices

  lat_indices = indices(:,0)                    ;-- latitude indices
  lon_indices = indices(:,1)                    ;-- longitude indices

  do i=0,ndims-1
     print(sprinti("%5i",lat_indices(i))+"  "+sprinti("%5i",lon_indices(i))) ;-- print ("lat index" "lon index")
  end do
end

EOF

  set -A ret $(ncl -nQ infile=\""${infile}"\" lat_array=\""${lats[*]}"\" lon_array=\""${lons[*]}"\" ncl_script)

  echo "latid lonid" 
  echo ${ret[0]} ${ret[1]}
  echo ${ret[2]} ${ret[3]}
  echo ${ret[4]} ${ret[5]}

exit

Hope this will help.
-Karin

Thanks a lot, Karin! I'll have to look into using NCL. (Actually, I'm building a Python application which already uses CDO for regridding and other things, so CDO will be nice. I'll look into NCL Python bindings ...)