[Thread Prev][Thread Next][Index]

Re: [ferret_users] current path length

Hi Leela,

This is not necessarily a better solution to your problem than use of transformations like @loc or @weq, but it illustrates how using Ferret as a digitizer, coupled with a Great Circle Distance calculator can be useful in a variety of situations. Apologies if this is already known to you but it may serve a wider audience.

First: Check that Ferret's "where" command is working for you:
    go ptest
When you move the cursor over the plot and click. Ferret should report the "mouse" coordinates in user units. This is the key to digitizing and a variety of other tasks one could imagine.

Next: Add the attached scripts to your directories (I would typically put the .jnl ones in my "ferret" directory and the shell script get_vertices [with chmod +x get_vertices] in "bin")

Demo: Here is an illustration. What is the length of the 25C isotherm of Levitus SST between Australia and South America?
    use levitus_climatology ; let sst=temp[k=1]
    fill/x=140:310/y=-30:10 sst ; cont/o/lev=(25)/x=140:310/y=-30:10 sst
    go polydef  ! prompts user to click on desired points, terminating by
! clicking offchart to the left. When this is done the coords of ! the point are loaded as a file with variables VX,VY
    sho data
   2> ./vertices.xy  (default)
 name     title                             I         J         K         L
VX VX ... ... 1:76 ... VY VY ... ... 1:76 ... ! As seen in the attached graphic I took 76 points to delineate the isotherm.

! In this case, where VX,VY are lon,lat values, what is needed is to sum the great
! circle distances along the curve. One way to do this is as follows:
    go greatcircle  ! definitions for great circle calculations
    let lon1=vx ; let lat1=vy ; let lon2=vx[k=@shf] ; let lat2=vy[k=@shf]
list gckm[k=1:75@sum] ! sum up the GC contributions, giving the answer (in kilometers)
Note that you sum over one less than the number of vertices so as not to run off the end. I did the above twice, because I forgot to save the plot, and the first time got 18042. This suggests the method is good enough for a certain class of problem. Perhaps a better demo would have been to use the hints in "greatcircle.jnl" to plot a portion of a great circle, then digitize it and compare the result with the known answer, but this one at least illustrates how it can deal with complex paths.

Hope it is of use,
On 5/19/11 11:06 AM, Leela Frankcombe wrote:
Dear Ferreters,

I've come across a question which is simple in theory but in practice I'm having a little trouble. So I'm wondering if anyone else has found a solution.

What I would like to do is calculate the path length of an ocean current. To define the path of the current I pick a particular sea surface height contour, now I would like to be able to calculate the length of that contour. I've been trying to use @loc or @weq to select the points at which SSH reaches the value that I've chosen but they only find the first instance along a given latitude or longitude so they miss parts where the current meanders (and also sometimes pick up eddies). Is there a way to use @loc or @weq to find every instance of a particular value? Or has someone got a better solution?


Leela Frankcombe
Post-doctoral researcher
Institute for Marine and Atmospheric research Utrecht
Utrecht University
The Netherlands

Attachment: demo.png
Description: PNG image

! polydef : use mouse click to define polygon vertices

can mode verify
let done=0 ; sp rm -f vertices.xy

say "****************************************************"
say "*                                                  *"
say "* Add polygon vertices by mouse clicks.  Terminate *"
say "*     by clicking to the left of the plot area.    *"
say "*                                                  *"
say "****************************************************"

! add new vertices to the file vertices.xy
repeat/range=1:1000 go add_vertex

! then read in the resulting file 
sp get_vertices ; go get_vertices

set mode verify

! add a new vertex to the vertex file while done=0

if `done eq 0` then 
  if `($XMOUSE) gt ($XAXIS_MIN)` then
    list/nohead/app/form=(2f12.6)/file=vertices.xy ($XMOUSE),($YMOUSE)
! if visible marks at the vertices are not desired drop the next line
    plot/o/nolab/vs/sym=1 ($XMOUSE),($YMOUSE)
    let done=1

! greatcircle : definitions for great circle calculations between
!               two locations lon1,lat1 and lon2,lat2 (in degrees)
let d2r=atan(1.)/45
let rlon1=d2r*lon1 ; let rlat1=d2r*lat1
let rlon2=d2r*lon2 ; let rlat2=d2r*lat2

! define great circle distances from lon1,lat1 to lon2,lat2 in radians ...
let/title="Great Circle Distance (radians)" \
! ... and kilometers
let/title="Great Circle Distance (km)" gckm=111.11*gcrad/d2r
let/title="Great Circle Distance (nm)" gcnm=60*gcrad/d2r

! define initial heading from lon1,lat1 to lon2,lat2 (clockwise from north)
let gcharg=acos((sin(rlat2)-sin(rlat1)*cos(gcrad))/(sin(gcrad)*cos(rlat1)))
! correct for near north-south pairings
let gchfix=gcharg[x=@fln]/d2r
let/title="Initial Heading (degrees)" \
   gchead=if(sin(rlon2-rlon1) gt 0)then gchfix else 360-gchfix

let/title="Secant Distance (3-D separation in km)" seckm=gckm*sin(gcrad/2)/(gcrad/2)

let lat=atan((tan(d2r*lat2)*sin(d2r*(lon-lon1))-tan(d2r*lat1)*sin(d2r*(lon-lon2)))/ \

! Usage: the results "gcrad", "gckm", and "gchead" are computed based on
!        existing variables lon1,lat1 representing the start point and
!        lon2,lat2 representing the destination.
! NOTE : To plot a greatcircle from lon1,lat1 to lon2,lat2 define a variable "lon"
!        that spans the interval between them with sufficient resolution, for example
!               def axis/x=`lon1`:`lon2`:0.1 xax ; let lon=x[gx=xax]
!               plot/o/vs/nolab/line=2 lon,lat
!        If a gridded file such as a topography is available it may be convenient to 
!        use that x-axis for "lon":
!                       let lon=x[g=...]
!               plot/o/vs/nolab/line=2/x=`lon1`:`lon2` lon,lat
#! /bin/csh -f

# get_vertices : create ferret script to access vertex data

set nrec=`wc -l vertices.xy | awk '{print $1}'`
echo 'def axis/z=1:'$nrec':1/mod zpoly ; def grid/z=zpoly gpoly' >! get_vertices.jnl
echo 'file/form=free/g=gpoly/var=vx,vy vertices.xy'              >> get_vertices.jnl

[Thread Prev][Thread Next][Index]
Contact Us
Dept of Commerce / NOAA / OAR / PMEL / Ferret

Privacy Policy | Disclaimer | Accessibility Statement