[postgis-users] finding street address points

[Dave Blasby]

Here's some perl code that take a linestring (WKT), percentage along the
linestring, side of the linestring, and offset distance to produce the
given point.

You'll have to figure out if the address is to the left or right of the
linestring, and how far down the linestring to travel (make sure you're
going in the right direction!).

dave


#distance_along( percentage along the geometry, WKT of geometry, side)
# side:: 0 = left, 1= right, 2 = unknown
#
# returns a WKT point
# dont call with srid=-1 ; if you do, make sure the pattern matchers
handle the "-"
# geometry can be any linestring (multi or single segment).  
# 0% = 1st point, 100% = last point in the linestring
#
# $pt_offset - global variable with the offset distance from the input
line to output point 
sub distance_along
{
	my ($perc, $geom,$side) = @_;

	my
(@points, at points_x, at points_y,$p,$t,$d,$last,$dist_along,$wanted_dist,$result);
	my ($r_x,$r_y,$rr_x_unit,
$rr_y_unit,$r_xx,$r_yy,$c_x,$c_y,$c_x_unit,$c_y_unit);

	#parse geometry
	
	$geom =~ s/SRID=[0123456789]+;LINESTRING\(//; # not "-" safe
	$geom =~ s/\)//;

	@points = split ',' , $geom;

	@points_x = ();
	@points_y = ();

	
	foreach $p (@points)
	{
		$p =~ /(.+) (.+)/;
		push @points_x,$1;
		push @points_y,$2;
	}

	$d = 0;
	foreach $t (1..$#points)
	{
		$d += sqrt(	(  ($points_x[$t]-$points_x[$t-1]) *
($points_x[$t]-$points_x[$t-1])  ) +
			(  ($points_y[$t]-$points_y[$t-1]) * ($points_y[$t]-$points_y[$t-1]) 
)			);
				
	}

	#printf ("total len = $d\n");
	
	#have points, and percent along.  need to find distance along
	# calc total length

	$wanted_dist = $d*$perc;
	foreach $t (1..$#points)
	{
		$d = sqrt(	(  ($points_x[$t]-$points_x[$t-1]) *
($points_x[$t]-$points_x[$t-1])  ) +
			(  ($points_y[$t]-$points_y[$t-1]) * ($points_y[$t]-$points_y[$t-1]) 
)			);
		if ($wanted_dist <=$d)
		{
			$r_x = $points_x[$t-1]  + ($points_x[$t]- $points_x[$t-1] )*
$wanted_dist/$d;
			$r_y = $points_y[$t-1]  + ($points_y[$t]- $points_y[$t-1] )*
$wanted_dist/$d;


				#unit vector of points[$t-1] to points[$t] 

				$r_xx = $points_x[$t]-$points_x[$t-1];
				$r_yy = $points_y[$t]-$points_y[$t-1];


			$r_xx_unit = ($r_xx) / sqrt( $r_xx*$r_xx + $r_yy*$r_yy ); 
			$r_yy_unit = ($r_yy)/ sqrt( $r_xx*$r_xx + $r_yy*$r_yy ); 

			$c_x = $r_yy_unit;
			$c_y = -1*$r_xx_unit;

			$c_x_unit = ($c_x) / sqrt( $c_x*$c_x + $c_y*$c_y ); 
			$c_y_unit = ($c_y)/ sqrt( $c_x*$c_x + $c_y*$c_y ); 
			
			if ($side==0)	#left
			{
				$r_x = $r_x - $pt_offset*$c_x_unit;
				$r_y = $r_y - $pt_offset*$c_y_unit;
			}
			if ($side ==1)
			{
				$r_x = $r_x + $pt_offset*$c_x_unit;
				$r_y = $r_y + $pt_offset*$c_y_unit;
			}

			
				

			$result= "SRID=1;POINT(";
			$result .= $r_x;
			$result .=" ";
			$result .= $r_y;
			return $result.")";
		}
		$wanted_dist -= $d;				
	}
	#went off end - just return the last point
	$result = "SRID=1;POINT(".$points_x[$#points_x]." ".
$points_y[$#points_y].")";
	return $result;
}