[postgis-users] st_distance accuracy?

Sufficool, Stanley ssufficool at sbcounty.gov
Tue Oct 20 13:46:32 PDT 2009


I am attempting to get the last vertex index of a line for a given
measure using the following:

idx :=  max(np) from generate_series(1,st_npoints(m_geom)) as np where
st_distance(st_pointn(m_geom,1), st_pointn(m_geom,np)) <= m_distance;

However, some of my queries are returning the last point on the
linestring when "m_distance" is less than the st_length() of the given
geometry.

Is this an issue with st_distance? I don't get the same issue when using
Lat/Lon and st_distance_spheroid().

select postgis_version()
1.4 USE_GEOS=1 USE_PROJ=1 USE_STATS=1

select version()
PostgreSQL 8.4.0, compiled by Visual C++ build 1400, 32-bit


FULL FUNCTION TEXT
------------------------
CREATE OR REPLACE FUNCTION st_line_interpolate_point(m_geom geometry,
m_percent double precision, m_offset double precision)
  RETURNS geometry AS
$BODY$

declare idx int;
declare m_distance double precision;
declare p1 geometry;
declare p2 geometry;

BEGIN

IF (m_distance > 1 or m_distance < 0) THEN
 raise exception 'Percentage must be between 0 and 1';
END IF;

/* Convert percentage to actual distance measure */
m_distance := st_length(m_geom) * m_percent;

/* Get the last point number on the line that is less than the requested
distance */

/* For Feet (NAD1983) */
idx :=  max(np) from generate_series(1,st_npoints(m_geom)) as np where
st_distance(st_pointn(m_geom,1), st_pointn(m_geom,np)) <= m_distance;

/* For Lat/Long (TIGER Lines) */
--idx :=  max(np) from generate_series(1,st_npoints(m_geom)) as np where
st_distance_sphere(st_pointn(m_geom,1), st_pointn(m_geom,np)) <
m_distance;

/* Past the end of the line segment return exception (this should never
happen) */
if idx = st_npoints(m_geom) then 
 raise exception 'st_line_interpolate_point(geom, dbl, dbl): Attempt to
seek beyond end of line.';
end if;

/* Grab the 2 points of the matching line segment */
p1:=st_pointn(m_geom,idx);
p2:=st_pointn(m_geom,idx + 1);

/* get the delta of the line segment and offset the interpolated point
*/
/* Rotate deltas by 90 degrees (invert dx/dy ) */
/* spheroid projections will need the offset to vary based on lat/lon */
return st_translate(
 line_interpolate_point(m_geom, m_percent), 
 ((st_y(p2) - st_y(p1)) / st_distance(p1,p2)) * m_offset,
 ((st_x(p2) - st_x(p1)) / st_distance(p1,p2)) * (-m_offset)
);

END;

$BODY$
  LANGUAGE 'plpgsql' VOLATILE
  COST 100;
ALTER FUNCTION st_line_interpolate_point(geometry, double precision,
double precision) OWNER TO postgres;




More information about the postgis-users mailing list