[postgis-users] function to convert meters into degrees without reprojections

[THEETEN Franck]

Hi Everybody,
 
A few days ago I had a problem with projection units. 
I had a series of points in WGS 84 (with a projection which uses decimal degrees as unit)  and I needed to trace a circle around them, with a radius specified in meter.
The 'Buffer()' function from the postgis library allows the definition of a circle around a geometry, but it didn't work in this situation because the radius value must be expressed in the same unit as the reference point.
I found on the Postgis mailing list some solutions which used the definition of a constant value allowing conversions between degrees and meters.
But it was an unsatisfactory answer for my problem because the result of this conversion varies with the latitude.
I also wanted to avoid reprojections into custom projections because I do not master well this type of operation. 
 
I defined a function in PL\pgSQL which more or less solve this problem (but with some error gap because conversion is based on sheperical earth). 
It defines a starting point which has the same latitude as the latitude entered as input parameter and a longitude equals to 0.
It defines another point which will "move" along the same parallel as the starting point.
It calls the 'distance_sphere' function from the postgis library, which is based on a spherical earth, in order to calculate the distance in meters between the source point and the moving point.
We repeat the operation until the resulting value becomes bigger or equals to the meter value to convert.
The last longitude value known for the "moving" point is returned as the result for the conversion into degree.
This function works and gives consistent values,  but is rather slow, because it use a loop with many iterations.
 
 
You can find  the source code below for comments, critics or suggestions about improvements:
 
Franck Theeten
Database Manager
Royal Museum For Central Africa
Tervuren, Belgium
 
 
---------------------------
 
--function which convert a distans in meter into a distance in degree
--Principe: the result varies with the latitude of the conversion
--Input parameter are : -'latitudevalue' latitude of the point for the considered distance
   --'meter': distance to convert
   -- 'precisionoutput': the precision of the output value (in decimal degrees)
   --'srid' srid of the considered projection 
--process:
--The function defines a starting point with has the same latitude as the latitude entered as input parameter and a longitude equals to 0.
--It akso defines another point which will "move" along the same parallel as the starting point (is given by precisionoutput parameter).
--We call the 'distance_sphere' function from the postgis library, which is based on a spherical earth.
--to calculate the distance in meters between the source point and the moving point,
--We repeat the operation until the resulting becomes bigger or equals to the meter value to convert.
--The last longitude value known for the "moving" point is returned as the result ot the conversion into degree .
--
 

CREATE OR REPLACE FUNCTION rmca_meterToDegreelatitude(latitudevalue double precision, meter double precision, precisionoutput double precision, srid int)
RETURNS numeric 
AS $$
DECLARE
 variablelong double precision; 
 tmpresult double precision;
 geometrystart geometry;
 geometrymove geometry;
BEGIN

 variablelong:=0;
 tmpresult:=0;
 
 geometrystart:=GEOMETRYFROMTEXT('POINT(0 '|| latitudevalue ||')', srid);
  
 WHILE tmpresult<=meter LOOP
  geometrymove := GEOMETRYFROMTEXT('POINT('|| variablelong ::varchar||' '|| latitudevalue::varchar||')', srid);
  
  tmpresult=distance_sphere(geometrystart, geometrymove);
  variablelong :=variablelong+precisionoutput;  
 END LOOP;
 
 RETURN variablelong;
END $$ LANGUAGE plpgsql;
 
 
-------example of invocation: SELECT rmca_meterToDegreelatitude(27.9, 1324, 0.00005, 32733); 
--(latitude of the point is 27.9°, 
--distance to convert is 1324m, 
--accuracy is 0.000005°
--srid of the projection is 32733 (WGS84 and UTM zone 33 South)
--result:-> 0.0134899999999996°
 

 
###########################################

This message has been scanned by ICT - Africa Museum

________________________________________
30/3/2007 - Filtered through antispam by ICT

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20070330/27170b49/attachment.html>