[postgis-users] st_distance accuracy?

[Sufficool, Stanley]

Correction inline

> -----Original Message-----
> From: postgis-users-bounces at postgis.refractions.net 
> [mailto:postgis-users-bounces at postgis.refractions.net] On 
> Behalf Of Sufficool, Stanley
> Sent: Tuesday, October 20, 2009 1:47 PM
> To: PostGIS Users Discussion
> Subject: [postgis-users] st_distance accuracy?
> 
> 
> 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;

Should read:

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

But it still errors at "Attempt to seek beyond..."

> 
> /* 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;
> 
> _______________________________________________
> postgis-users mailing list postgis-users at postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users
>