[postgis-users] RE: Conversion of Lambert72 to WGS84 and vice-versa

Guy Thomas gt at fks.be
Tue May 5 23:49:24 PDT 2009


Hmm, I'm definitely using an older version of PostGIS. Have to try the 
newer version. Thank you for pointing this out to me.

SELECT postgis_full_version();
                                postgis_full_version 

----------------------------------------------------------------------------------
  POSTGIS="1.2.1" GEOS="3.0.0-CAPI-1.4.1" PROJ="Rel. 4.6.0, 21 Dec 2007" 
USE_STATS


SELECT proj4text from spatial_ref_sys where srid = 31370;
 
  proj4text 

---------------------------------------------------------------------------------------------------------------------------------------------------------
  +proj=lcc +lat_1=51.16666723333333 +lat_2=49.8333339 +lat_0=90 
+lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438 +ellps=intl 
+units=m +no_defs



-----

Date: Tue, 5 May 2009 13:46:32 -0400
From: "Paragon Corporation" <lr at pcorp.us>
Subject: [postgis-users] 	RE: Conversion of Lambert72 to WGS84 and,
	vice-versa
To: <postgis-users at postgis.refractions.net>
Message-ID: <667F639A5162470E875576AD2CAB4519 at b>
Content-Type: text/plain;	charset="us-ascii"

Guy,

Which version of proj are you running in PostGIS.

I tried the below and it comes closer to the Geotools answer you have.

SELECT ST_AsText(ST_Transform('SRID=31370;POINT(148378.77 172011.96)',4326))

Gives:

POINT(4.34572624388819 50.8584956187369)


I'm running:
SELECT postgis_full_version();

"POSTGIS="1.3.5" GEOS="3.1.0-CAPI-1.5.0" PROJ="Rel. 4.6.1, 21 August 2008"
USE_STATS"

So could be your proj for PostGIS is different or the entry in
spatial_ref_sys

Mine returns
SELECT proj4text from spatial_ref_sys where srid = 31370;

"+proj=lcc +lat_1=51.16666723333333 +lat_2=49.8333339 +lat_0=90
+lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438 +ellps=intl
+units=m
+towgs84=-106.868628,52.297783,-103.723893,0.336570,-0.456955,1.842183,-1
+no_defs "

Leo

-----Original Message-----
From: Guy Thomas [mailto:gt at fks.be]
Sent: Tuesday, May 05, 2009 5:41 AM
To: postgis-users at postgis.refractions.net
Cc: lr at pcorp.us
Subject: Conversion of Lambert72 to WGS84 and, vice-versa

Regina,

These are the outputs of the three different tools I used:

The input coordinate (long, lat) in Lambert72 (EPSG:31370): 148378.77
172011.96

In WGS84: ogr2ogr: 4.343195763754792  50.857974812545223

            via Geotools: 4.345726237889486  50.85849524233313

            via postGIS: 4.34446096709109 50.859039939172

[sql-statement: select node_id, AsText (node_geom),
AsText(st_transform(node_geom, 4326)) from nodes;

   node_id |           astext           |                 astext

---------+----------------------------+---------------------------------
---------+----------------------------+--------
       120 | POINT(148378.77 172011.96) | POINT(4.34446096709109
50.859039939172)
]

Best regards,

   Guy

Further spec:

1. ogr2ogr: I had to convert from csv to csv; as this didn't seem to work, I
converted first to gml and then back to csv:

ogr2ogr -overwrite -f GML -t_srs 'EPSG:4326' -s_srs 'EPSG:31370'
havenbrusselgegevens.gml havenbrusselgegevens.vrt

ogr2ogr --debug on -overwrite -f CSV -t_srs 'EPSG:4326' -s_srs 'EPSG:4326'
out havenbrusselgegevens.gml -sql 'select *,OGR_GEOM_WKT from
havenbrusselgegevens'

The second command is only there to get the WGS83 coordinates (generated in
the first step) in the csv file. I checked: the second step simply copies
the original and the converted coordinates in the output csv file.

The havenbrusselgegevens.vrt file:

<OGRVRTDataSource>
      <OGRVRTLayer name="havenbrusselgegevens">
          <SrcDataSource
relativeToVRT="0">havenbrusselgegevens.csv</SrcDataSource>
          <SrcLayer>havenbrusselgegevens</SrcLayer>
          <GeometryType>wkbPoint</GeometryType>
          <GeometryField encoding="PointFromColumns" x="Longitude"
y="Latitude" />
          <LayerSRS>epsg:31370</LayerSRS>
      </OGRVRTLayer>
</OGRVRTDataSource>

The first lines of the havenbrusselgegevens.csv file:

Longitude,Latitude,Location
148378.77,172011.96,"P.Place Sainctelette"
148389.51,171988.13,"(leeg)"
148416.83,172094.99,"Saictelettebrug"
148439.62,172070.48,"(leeg)"
...

The first lines of havenbrusselgegevens.gml:

<?xml version="1.0" encoding="utf-8" ?>
<ogr:FeatureCollection
       xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
       xsi:schemaLocation="http://ogr.maptools.org/havenbrusselgegevens.xsd"
       xmlns:ogr="http://ogr.maptools.org/"
       xmlns:gml="http://www.opengis.net/gml">
    <gml:boundedBy>
      <gml:Box>

<gml:coord><gml:X>4.296037795179485</gml:X><gml:Y>50.81251473822174</gml:Y><
/gml:coord>

<gml:coord><gml:X>4.411269648572022</gml:X><gml:Y>50.92352164264287</gml:Y><
/gml:coord>
      </gml:Box>
    </gml:boundedBy>
    <gml:featureMember>
      <ogr:havenbrusselgegevens fid="F0">
        <ogr:geometryProperty><gml:Point
srsName="EPSG:4326"><gml:coordinates>4.343195763754792,50.857974812545223,34
7.641519372361927</gml:coordinates></gml:Point></ogr:geometryProperty>
        <ogr:Longitude>148378.77</ogr:Longitude>
        <ogr:Latitude>172011.96</ogr:Latitude>
        <ogr:Location>P.Place Sainctelette</ogr:Location>
      </ogr:havenbrusselgegevens>
    </gml:featureMember>
...

The first lines of the final output file:
Longitude,Latitude,Location,OGR_GEOM_WKT
148378.77,172011.96,P.Place Sainctelette,POINT (4.343195763754792
50.857974812545223 347.641519372361927)
148389.51,171988.13,(leeg),POINT (4.343348401648505 50.857760633655566
347.641805360402032)
148416.83,172094.99,Saictelettebrug,POINT (4.34373592745083
50.858721283251469 347.63981045788023)
...


2. My Geotools wrapper code:

[In Geotools the JTS library is used for the transformation and I do not
know how good or bad the conversion code is.]

try
        {
           Hints hints = new
Hints(Hints.FORCE_LONGITUDE_FIRST_AXIS_ORDER, Boolean.TRUE);
           CRSAuthorityFactory factory =
ReferencingFactoryFinder.getCRSAuthorityFactory("EPSG", hints);
           CoordinateReferenceSystem crsWGS84 =
factory.createCoordinateReferenceSystem("EPSG:4326");
           //CoordinateReferenceSystem crsWGS84 = CRS.decode("EPSG:4326");
           CoordinateReferenceSystem crsLambert72 = 
CRS.decode("EPSG:31370");

           System.out.println("WGS84: " + crsWGS84.toWKT());
           System.out.println("Lambert72: " + crsLambert72.toWKT());

           // transforming from LAMBERT72 to WGS84
           MathTransform transformLambert72ToWGS84 =
CRS.findMathTransform(crsLambert72, crsWGS84);

           Coordinate c10 = new Coordinate(148378.77, 172011.96);
           Coordinate c11 = new Coordinate();
           JTS.transform(c10, c11, transformLambert72ToWGS84);
           System.out.println("c10 : " + c10);
           System.out.println("c11 : " + c11);
        }
        catch (NoSuchAuthorityCodeException nsae)
        {
           nsae.printStackTrace();
        }
        catch (FactoryException fe)
        {
           fe.printStackTrace();
        }
        catch (TransformException te)
        {
           te.printStackTrace();
        }

3. Postgres-PostGIS:

create table nodes (
     node_id    integer,
     lc_country text
);

select AddGeometryColumn('', 'nodes', 'node_geom', 31370, 'POINT', 2);

insert into nodes (node_id, node_geom, lc_country) values (120,
GeomFromText('POINT(148378.77  172011.96)', 31370), 'BE');

select node_id, AsText (node_geom), AsText(st_transform(node_geom,
4326)) from nodes;




More information about the postgis-users mailing list