[Thread Prev][Thread Next][Index]
Re: [ferret_users] Sample a field along a set of line segments?
Hi Billy,
Initially i got a method to solve it using Fortran. But
it was a lengthy one. My friend Bharath suggested a method based
on GCD (almost same as the solution i got but a neat one). Here
is the Ferret implimentation of the solution as a "GO" file
(please see attached jnl file grid_pts_along_line_segments.jnl).
Comments/suggestions are most welcome. If you face any problems
with this GO file, I will be happy to look into it.
If anybody is interested to have a look at the Fortran stuff,
please send me a mail.
Hope this helps,
Regards,
Jaison
On Tue, 3 Oct 2006, William S. Kessler wrote:
> I want to sample a field along a line defined by a set of (x,y)
> points. Clearly I could use SAMPLEXY to extract the values at the
> given points, but I want to, in effect, draw straight lines between
> those points and sample all along those lines.
>
> It is necessary to say how densely I want to sample, that is, how
> many intermediate points to add in between the given ones. Suppose I
> want to add an intermediate point at each value of the x-grid of the
> field to be sampled. As an example, suppose I have two points (155,7)
> and (158,13), and my background grid is 1 by 1, on integers. The
> points I want to add would be (156,9) and (157,11). How can I find
> those points and augment my original set to include them? Then I can
> use SAMPLEXY to sample the field.
>
> Billy K
> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
> William S. Kessler
> NOAA / Pacific Marine Environmental Laboratory
> 7600 Sand Point Way NE
> Seattle WA 98115 USA
>
> william.s.kessler@noaa.gov
> Tel: 206-526-6221
> Fax: 206-526-6744
> Home page: http://www.pmel.noaa.gov/people/kessler/
>
>
\ cancel mode verify
!
! Go file to find the grid points along specified line segments. In other
! words to find the grid points on which a cruise track will pass through.
!
! Assumptions : 1. Each Cruise track station (x,y) points exists in the gridded
! dataset.
! 2. Gridded data have "uniform" dx and dy
!
! Usage : GO grid_pts_along_line_segments stn_x, stn_y, dx, dy, points
!
! Arguments
!
! $1 : mandatory : X(longitude) points of cruise stations, on X-grid : INPUT
! $2 : " : Y(latitude) points of cruise stations, on X-grid : "
! $3 : " : dX, spacing along X(Lon) for gridded data : "
! $4 : " : dY, spacing along Y(Lat) for gridded data : "
! $5 : optional : total number of grid points, including new ones : OUPUT
!
! This GO file will write all the (X,Y) pairs of new as well as old grid
! points to a text file grid_pts.txt, which can be loaded into Ferret
! in order to sample the gridded data along the cruise track (without
! any interpolation error) using SAMPLEXY external function.
!
! NOTE : All the internal variables of this script have a prefix "gpl"
! except "npts".
!
! Example :
! use levitus_climatology ! d=1
! define axis/x=0:360:1/units=longitudes xlon ! for the sake of
! define axis/y=-90:90:1/units=latitudes ylat ! example
! let sst = temp[d=1,k=1,gx=xlon,gy=ylat] ! d=1 is important as we
! let dx = x[gx=sst,i=2]-x[gx=sst,i=1] ! will load a new
! let dy = y[gy=sst,j=2]-y[gy=sst,j=1] ! file soon
! !--Cruise track : station lat-lon points : defined on X-Axis (A MUST)
! let stn_x = XSEQUENCE({155,158,158,160,157,160,162,165,168,165})
! let stn_y = XSEQUENCE({ 7, 10, 14, 16, 19, 22, 24, 24, 21, 18})
!
! GO grid_pts_along_line_segments stn_x, stn_y, dx, dy, points
!
! define axis/x=1:`points`:1 xfile
! define grid/x=xfile gfile
! FILE/grid=gfile/var="xpts, ypts" grid_pts.txt ! d=2
! fill/x=120:200/y=0:30/pal=inverse_greyscale/lev=(14,30,0.5) sst
! plot/vs/ov/nolab/size=0.12/symb=1/color=8 xpts[d=2], ypts[d=2] ! new grid points in RED
! plot/vs/ov/nolab/size=0.12/symb=1/color=11 stn_x, stn_y
! go land
! pause
! let sst_cruise = SAMPLEXY(sst,xpts[d=2],ypts[d=2])
! plot sst_cruise
!
!---------------------------------------------------------------------------------------
!
! Created By : Jaison Kurian (jaison@caos.iisc.ernet.in) &
! Bharath Raj G. N. (ozone@caos.iisc.ernet.in)
!
! Created on : 06/Oct/2006
!
! Modifications : None
!----------------------------------------------------------------------------------------
query/ignore $1%<Use: GO grid_pts_along_line_segments xpts, ypts, dx, dy, npts%
query/ignore $2%<Use: GO grid_pts_along_line_segments xpts, ypts, dx, dy, npts%
query/ignore $3%<Use: GO grid_pts_along_line_segments xpts, ypts, dx, dy, npts%
query/ignore $4%<Use: GO grid_pts_along_line_segments xpts, ypts, dx, dy, npts%
query/ignore $1%<Use: GO grid_pts_along_line_segments xpts, ypts, dx, dy, npts%
let gpl_stn_x = $1
let gpl_stn_y = $2
let gpl_dx = $3
let gpl_dy = $4
!--do some consistency checks
def sym x_shape = `gpl_stn_x,r=shape`
def sym y_shape = `gpl_stn_y,r=shape`
IF `"($x_shape)" NE "X"` THEN
SAY " ERROR : 1st argument (xpts) to grid_pts_along_line_segments should be "
SAY " defined on X-axis (XSEQUENCE will make it easy) "
EXIT
ENDIF
IF `"($y_shape)" NE "X"` THEN
SAY " ERROR : 2nd argument (ypts) to grid_pts_along_line_segments should be "
SAY " defined on X-axis (XSEQUENCE will make it easy) "
EXIT
ENDIF
let gpl_stns = `gpl_stn_x,r=iend`
let gpl_stns1 = `gpl_stn_y,r=iend`
IF `gpl_stns NE gpl_stns1` THEN
SAY " ERROR : From grid_pts_along_line_segments"
SAY " Number of xpts and ypts (1st & 2nd Arguments) are not same"
EXIT
ENDIF
IF `gpl_stns LT 2` THEN
SAY " ERROR : From grid_pts_along_line_segments"
SAY " Number of xpts and ypts (1st & 2nd Arguments) should be > 1."
EXIT
ENDIF
IF `gpl_dx LE 0 OR gpl_dy LE 0` THEN
SAY " ERROR : From grid_pts_along_line_segments"
SAY " dx and dy (3rd & 4th Arguments) should be > 0."
EXIT
ENDIF
!--if everything is ok, then do it
sp rm -f grid_pts.txt
define symbol gpl_qual = quiet/nohead/file=grid_pts.txt/format=(2x,2(f9.4,3x))/APPEND
let gpl_points = 1
list/($gpl_qual) gpl_stn_x[i=1], gpl_stn_y[i=1]
REPEAT/RANGE=1:`gpl_stns-1`:1/NAME=gpl_st (;\
let gpl_xs = gpl_stn_x[i=`gpl_st`] ; let gpl_xe = gpl_stn_x[i=`gpl_st+1`] ;\
let gpl_ys = gpl_stn_y[i=`gpl_st`] ; let gpl_ye = gpl_stn_y[i=`gpl_st+1`] ;\
let gpl_xn = ((gpl_xe - gpl_xs)^2)^0.5 / gpl_dx ;\
let gpl_yn = ((gpl_ye - gpl_ys)^2)^0.5 / gpl_dy ;\
let gpl_xo = (gpl_xe - gpl_xs) * gpl_dx/MAX(gpl_xn,1) ;\
let gpl_yo = (gpl_ye - gpl_ys) * gpl_dy/MAX(gpl_yn,1) ;\
IF `gpl_xn EQ 0 OR gpl_yn EQ 0` THEN ;\
let gpl_gcd = MAX(gpl_xn,gpl_yn) ;\
ELSE ;\
define axis/x=1:`MIN(gpl_xn,gpl_yn)`:1 gpl_xax ;\
let gpl_d = x[gx=gpl_xax] ;\
let gpl_devisors_x = IF gpl_xn/gpl_d EQ INT(gpl_xn/gpl_d) THEN gpl_d ;\
let gpl_devisors_y = IF gpl_yn/gpl_d EQ INT(gpl_yn/gpl_d) THEN gpl_d ;\
let gpl_devisors = gpl_devisors_x + gpl_devisors_y*0 ;\
let gpl_gcd = gpl_devisors[i=@MAX] ;\
ENDIF ;\
let gpl_xinc = gpl_xn/gpl_gcd ; let gpl_yinc = gpl_yn/gpl_gcd ;\
REPEAT/RANGE=1:`gpl_gcd`:1/NAME=gpl_in ( ;\
let gpl_x2 = gpl_xs + gpl_xinc * gpl_in * gpl_dx * gpl_xo ;\
let gpl_y2 = gpl_ys + gpl_yinc * gpl_in * gpl_dy * gpl_yo ;\
list/($gpl_qual) gpl_x2, gpl_y2 ;\
let gpl_points = `gpl_points` + 1 ;\
);\
)
define symbol gpl_npts = `gpl_points`
let $5"npts" = `($gpl_npts)`
cancel var gpl_*
SAY
SAY Output is avialable in the text file "grid_pts.txt".
SAY
[Thread Prev][Thread Next][Index]
Dept of Commerce /
NOAA /
OAR /
PMEL /
TMAP
Contact Us | Privacy Policy | Disclaimer | Accessibility Statement