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: [email protected] [mailto:[email protected]] On Behalf Of Guy Thomas Sent: Wednesday, May 06, 2009 2:49 AM To: [email protected] 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" <[email protected]> Subject: [postgis-users] RE: Conversion of Lambert72 to WGS84 and, vice-versa To: <[email protected]> Message-ID: <667f639a5162470e875576ad2cab4...@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:[email protected]] Sent: Tuesday, May 05, 2009 5:41 AM To: [email protected] Cc: [email protected] 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 [email protected] http://postgis.refractions.net/mailman/listinfo/postgis-users _______________________________________________ postgis-users mailing list [email protected] http://postgis.refractions.net/mailman/listinfo/postgis-users
