[postgis-users] RE: Conversion of Lambert72 to WGS84 and vice-versa
Paragon Corporation
lr at pcorp.us
Wed May 6 00:26:44 PDT 2009
Guy,
Well your proj4 install looks fairly new which is the important thing for
transformation.
The fastest thing to do might be to just update your spatial_ref_sys
proj4text for srid = 31370 with what I have. Looks like yours is missing
the towgs84 etc. which I suspect might do the trick for you.
"+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 "
If it is not an issue -- yes upgrading is better since your postgis and geos
install are a little behind.
Leo
-----Original Message-----
From: postgis-users-bounces at postgis.refractions.net
[mailto:postgis-users-bounces at postgis.refractions.net] On Behalf Of Guy
Thomas
Sent: Wednesday, May 06, 2009 2:49 AM
To: postgis-users at postgis.refractions.net
Subject: [postgis-users] RE: Conversion of Lambert72 to WGS84 and vice-versa
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;
_______________________________________________
postgis-users mailing list
postgis-users at postgis.refractions.net
http://postgis.refractions.net/mailman/listinfo/postgis-users
More information about the postgis-users
mailing list