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

Paragon Corporation lr at pcorp.us
Tue May 5 10:46:32 PDT 2009


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;

-- 
Guy Thomas                                    guy.thomas at fks.be
fks bvba - Formal and Knowledge Systems       http://www.fks.be/
Stationsstraat 108                            Tel:  ++32-(0)11-21 49 11
B-3570 ALKEN                                  Fax:  ++32-(0)11-22 04 19






More information about the postgis-users mailing list