# 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
where
```
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)
```    18033.
```
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,
Mick
----------------
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?

Thanks!
Leela.

-----------------------------------------------------------------------------------
Leela Frankcombe
Post-doctoral researcher
Institute for Marine and Atmospheric research Utrecht
Utrecht University
The Netherlands
www.phys.uu.nl/~frankcmb
l.m.frankcombe@xxxxx
-----------------------------------------------------------------------------------
```
```
```

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

! 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
where
if `(\$XMOUSE) gt (\$XAXIS_MIN)` then
! if visible marks at the vertices are not desired drop the next line
plot/o/nolab/vs/sym=1 (\$XMOUSE),(\$YMOUSE)
else
let done=1
endif
endif

```
```! 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)
! 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)))/ \
sin(d2r*(lon2-lon1)))/d2r

! 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
```