[postgis-users] Spherical expand

[Rob Young]

Hey guys,

A while ago I mentioned that I had written a postgres function that
acted as a expand function that worked with lat/long data.  I thought
I'd pass it on to the list, in case other people might find it useful.

Included are two functions: expand_sphere and distance_from.
distance_from is a helper function for expand_sphere.  expand_sphere
creates a polygon around the lat and long coordinates passed to it,
the size of which is determined by the radius value you give it (in
metres.)  It uses a spherical earth approximation for the underlying
math.  Currently it doesn't handle wraparound.  So, it could in theory
create a box that bounds outside of the normal lat/long bounds.  (I'm
on vacation, I'll get around to writing the boundery cases when I get
back home.  Also, more commenting then...)  :)

Cheers!

Rob Young

Ps: hopefully the attached code is preserved in some useable form..  :)


---- See => http://williams.best.vwh.net/avform.htm#LL
----      lat=asin(sin(lat1)*cos(d)+cos(lat1)*sin(d)*cos(tc))
----      IF (cos(lat)=0)
----         lon=lon1      // endpoint a pole
----      ELSE
----         lon=mod(lon1-asin(sin(tc)*sin(d)/cos(lat))+pi,2*pi)-pi
----      ENDIF

CREATE OR REPLACE FUNCTION distance_from( IN lat real, -- Latitude in degrees
					IN long real, -- Longitude in degrees
					IN course real, -- course in degrees
					IN dist real ) -- distance in metres
RETURNS geometry AS $$
DECLARE
	radius_km real; -- CONSTANT radius of the earth
	latRad real;
	longRad real;
	courseRad real;
	distRad real;

	newLatRad real;
	newLongRad real;

	newLat real;
	newLong real;
BEGIN
	radius_km = 6378.137; -- radius of the earth (ROUGHLY, CHECK!)

	latRad := ( pi() / 180.0 ) * lat;
	longRad := ( pi() / 180.0 ) * long;
	courseRad := ( pi() / 180.0 ) * course;
	distRad := ( dist / 1000.0 ) / radius_km;

	newLatRad := asin(sin(latRad)*cos(distRad)+cos(latRad)*sin(distRad)*cos(courseRad));
	newLongRad := 0.0;


	IF cos( newLatRad ) = 0 THEN
		newLongRad = longRad;
	ELSE
		--newLongRad =
mod(longRad-asin(sin(courseRad)*sin(distRad)/cos(lat))+pi(), 2 * pi())
- pi();
		newLongRad = cast(longRad-asin(sin(courseRad)*sin(distRad)/cos(latRad))+pi()
as numeric ) % cast( 2 * pi() as numeric ) - pi();
	END IF;

	newLat := ( 180.0 / pi() ) * newLatRad;
	newLong := ( 180.0 / pi() ) * newLongRad;

	RETURN SetSRID( MakePoint( newLong, newLat ), 4326 );
END;
$$ LANGUAGE plpgsql IMMUTABLE;

CREATE OR REPLACE FUNCTION expand_sphere( IN lat real, -- Latitude in degrees
					IN long real, -- Longitude in degrees
					IN dist real ) -- distance in metres
RETURNS geometry AS $$
DECLARE
	a geometry;
	b geometry;
	c geometry;
	d geometry;

	llPoint geometry;
	urPoint geometry;

	sphereBox box2D;
BEGIN
	a := distance_from( lat, long, 0.0, dist );
	b := distance_from( lat, long, 90.0, dist );
	c := distance_from( lat, long, 180.0, dist );
	d := distance_from( lat, long, 270.0, dist );

	llPoint := MakePoint( X( d ), Y( c ) );
	urPoint := MakePoint( X( b ), Y( a ) );

	sphereBox := MakeBox2D( llPoint, urPoint );

	RETURN SetSRID( sphereBox, 4326 );
END;
$$ LANGUAGE plpgsql IMMUTABLE;
-------------- next part --------------
---- See => http://williams.best.vwh.net/avform.htm#LL
----      lat=asin(sin(lat1)*cos(d)+cos(lat1)*sin(d)*cos(tc))
----      IF (cos(lat)=0)
----         lon=lon1      // endpoint a pole
----      ELSE
----         lon=mod(lon1-asin(sin(tc)*sin(d)/cos(lat))+pi,2*pi)-pi
----      ENDIF

CREATE OR REPLACE FUNCTION distance_from( IN lat real, -- Latitude in degrees
					IN long real, -- Longitude in degrees
					IN course real, -- course in degrees
					IN dist real ) -- distance in metres
RETURNS geometry AS $$
DECLARE 
	radius_km real; -- CONSTANT radius of the earth
	latRad real;
	longRad real;
	courseRad real;
	distRad real;

	newLatRad real;
	newLongRad real;

	newLat real;
	newLong real;
BEGIN 
	radius_km = 6378.137; -- radius of the earth (ROUGHLY, CHECK!)

	latRad := ( pi() / 180.0 ) * lat;
	longRad := ( pi() / 180.0 ) * long;
	courseRad := ( pi() / 180.0 ) * course;
	distRad := ( dist / 1000.0 ) / radius_km;

	----      lat=asin(sin(lat1)*cos(d)+cos(lat1)*sin(d)*cos(tc))
	newLatRad := asin(sin(latRad)*cos(distRad)+cos(latRad)*sin(distRad)*cos(courseRad));
	newLongRad := 0.0;

----      IF (cos(lat)=0)
----         lon=lon1      // endpoint a pole
----      ELSE
----         lon=mod(lon1-asin(sin(tc)*sin(d)/cos(lat))+pi,2*pi)-pi
----      ENDIF

	IF cos( newLatRad ) = 0 THEN
		newLongRad = longRad;
	ELSE
		--newLongRad = mod(longRad-asin(sin(courseRad)*sin(distRad)/cos(lat))+pi(), 2 * pi()) - pi();
		newLongRad = cast(longRad-asin(sin(courseRad)*sin(distRad)/cos(latRad))+pi() as numeric ) % cast( 2 * pi() as numeric ) - pi();
	END IF;

	newLat := ( 180.0 / pi() ) * newLatRad;
	newLong := ( 180.0 / pi() ) * newLongRad;

	RETURN SetSRID( MakePoint( newLong, newLat ), 4326 );
END;
$$ LANGUAGE plpgsql IMMUTABLE;

CREATE OR REPLACE FUNCTION expand_sphere( IN lat real, -- Latitude in degrees
					IN long real, -- Longitude in degrees
					IN dist real ) -- distance in metres
RETURNS geometry AS $$
DECLARE 
	a geometry;
	b geometry;
	c geometry;
	d geometry;

	llPoint geometry;
	urPoint geometry;

	sphereBox box2D;
BEGIN 
	a := distance_from( lat, long, 0.0, dist );
	b := distance_from( lat, long, 90.0, dist );
	c := distance_from( lat, long, 180.0, dist );
	d := distance_from( lat, long, 270.0, dist );

	llPoint := MakePoint( X( d ), Y( c ) );
	urPoint := MakePoint( X( b ), Y( a ) );

	sphereBox := MakeBox2D( llPoint, urPoint );

	RETURN SetSRID( sphereBox, 4326 );
END;
$$ LANGUAGE plpgsql IMMUTABLE;