; ncl-talk
; I have latitude and longitude in cartesian coordinates (example of the text file is above). 
; That is the only information I have. I plotted the data before using gsn_add_polymarker 
; for the domain I want e.g. 15N 20N and 54W-58W. The center of this domain is 16.8N and 56W. 
; What the polymaker function did was if I had an observation at 16N and 53W a black dot 
; was drawn. What I need now is transform that 16N-53W (black dot) into cylindrical 
; coordinates using as center 16.8N and 56W. 
;
  diri = "./"
  fili = "station.txt" 
  pthi = diri+fili
  sdat = asciiread(pthi, -1, "string")
    
; sdat(0):  2010/08/28,04:00:00.356732, 11.7047, 133.9177,  4.8,  5
;           yyyy mm dd hh mn ----sc---   lat       lon      val   nnn
; field #:     1  2  3  4  5   6          7         8        9    10  

  delim = "/:,"
  nfld  = str_fields_count(sdat(0), delim)  ; use 1st line
  print(nfld)

  yyyy = toint(str_get_field(sdat,   1, delim))
  mm   = toint(str_get_field(sdat,   2, delim))
  dd   = toint(str_get_field(sdat,   3, delim))
  hh   = toint(str_get_field(sdat,   4, delim))
  mn   = toint(str_get_field(sdat,   5, delim))
  sc   = tofloat(str_get_field(sdat, 6, delim))
  lat  = tofloat(str_get_field(sdat, 7, delim))
  lon  = tofloat(str_get_field(sdat, 8, delim))
  val  = tofloat(str_get_field(sdat, 9, delim))
  nnn  = toint(str_get_field(sdat,  10, delim))

; User specified geograpcal areah

  clat =  16.8    ; center point
  clon = -56.0

  latS =  15.0    ; area of interest
  latN =  20.0
  lonL = -58.0
  lonR = -54.0

  print("clat="+clat+"  clon="+clon)
  print("---")

; Indices of data within the specified region 

  ii   = ind(lat.ge.latS .and. lat.le.latN .and. \
             lon.ge.lonL .and. lon.le.lonR)

  if (ismissing(ii(0))) then
      print("No observations in specified region") 
      exit
  else
       print("---")
       print("Obsverations  being used:")
       print("---")
       print(""+sdat(ii))
       print("---")
  end if

  nii  = dimsizes(ii)
  print("There are "+nii+" observations in specified region") 
  print("---")

; convenience: extract and overwrite using NCL reassignment operator   :=

  yyyy := yyyy(ii)
  mm   := mm(ii)
  dd   := dd(ii)
  hh   := hh(ii)
  mn   := mn(ii)
  sc   := sc(ii)
  lat  := lat(ii)
  lon  := lon(ii)
  val  := val(ii)
  nnn  := nnn(ii)

; Compute distance from central location to selected observations

  R     = gc_latlon(clat,clon,lat,lon,2,4)   ; great circle; polar radius   
  R@long_name = "distance (radius) from central location to each point"
  R@units     = "km"
  R@information = "great circle : ("+sprintf("%6.2f", clat)+","+sprintf("%6.2f", clon)+")"
  print(R)
  print("---")

  dlat  = (lat-clat)*110.567            ; degrees -> km
  dlon  = (lon-clon)*110.567            ; degrees -> km ... ok region of interest
  RC    = sqrt(dlon^2 + dlat^2)         ; hypotenuse:  km   
  RC@long_name = "Cartesian distance (radius) from central location to each point"
  RC@units     = "km"
  RC@information = "Cartesian approximation: ("+sprintf("%6.2f", clat)+","+sprintf("%6.2f", clon)+")"
  print(RC)
  print("---")

  theta = atan2(dlat,dlon)*57.29578     ; degrees
  theta@long_name = "theta: angle from central point"
  theta@units     = "degrees"
  print("R="+R+"  RC="+RC+"  ratio="+(RC/R)+"  theta="+theta)