coordinate systemproj
I'm trying to convert a set of UTM32N coordinates to a local coordinate system with the same datum.
I have the four Helmert parameters (x,y,scale,theta) that defines the coordinate system.
How can I apply this transformation using proj4 strings?
I've tried the following one but it differs from the expected result:
proj +proj=helmert +ellps=WGS84 +x=-711066.2010726777 +y=-4532044.137231419 +s=0.999785911377552 +theta=-2.089700973339795
Your Gauss-Krueger projection uses +datum=potsdam. Up to 2012, this was hard coded in proj4 to a very unprecise value using a 3-parameter-transformation.
You find more exact values for 7-parameter transformations in this topic:
http://forum.openstreetmap.org/viewtopic.php?id=12723
There is an even better ntv2-grid transformation available here (take the binary), that has to be in the same folder as your application and data, unlsss you specify full pathnames.
To compare the different possible values, I made this test batch file:
echo epsg31467-epsg4326 >out.txt
cs2cs +init=epsg:31467 +to +init=epsg:4326 31467.txt >>out.txt
echo proj-Definition epsg >>out.txt
cs2cs +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs +to +init=epsg:4326 31467.txt >>out.txt
echo proj-definition Qgis >>out.txt
cs2cs +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel +towgs84=582,105,414,1.04,0.35,-3.08,8.3 +units=m +no_defs +to +init=epsg:4326 31467.txt >>out.txt
echo proj-definition nadgrid >>out.txt
cs2cs +proj=tmerc +lat_0=0 +lon_0=9 +x_0=3500000 +y_0=0 +k=1.000000 +ellps=bessel +units=m +nadgrids=./BETA2007.gsb +wktext +to +init=epsg:4326 31467.txt >>out.txt
echo epsg31467-epsg3785 >>out.txt
cs2cs +init=epsg:31467 +to +init=epsg:3785 31467.txt >>out.txt
echo proj-definition Qgis >>out.txt
cs2cs +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel +towgs84=582,105,414,1.04,0.35,-3.08,8.3 +units=m +no_defs +to +init=epsg:3785 31467.txt >>out.txt
echo proj-definition nadgrid >>out.txt
cs2cs +proj=tmerc +lat_0=0 +lon_0=9 +x_0=3500000 +y_0=0 +k=1.000000 +ellps=bessel +units=m +nadgrids=./BETA2007.gsb +wktext +to +init=epsg:3785 31467.txt >>out.txt
echo epsg31467-epsg3857 >>out.txt
cs2cs +init=epsg:31467 +to +init=epsg:3857 31467.txt >>out.txt
echo proj-definition Qgis >>out.txt
cs2cs +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel +towgs84=582,105,414,1.04,0.35,-3.08,8.3 +units=m +no_defs +to +init=epsg:3857 31467.txt >>out.txt
echo epsg31467-epsg900913 >>out.txt
cs2cs +init=epsg:31467 +to +init=epsg:900913 31467.txt >>out.txt
echo epsg31467-proj900913 >>out.txt
cs2cs +init=epsg:31467 +to +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +over +no_defs 31467.txt >>out.txt
echo proj31467-proj900913 >>out.txt
cs2cs +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel +towgs84=582,105,414,1.04,0.35,-3.08,8.3 +units=m +no_defs +to +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +over +no_defs 31467.txt >>out.txt
with any sample Gauss-Krüger lon-lat coordinate pair in 31467.txt
The UTM 23S grid aligns with a degree grid only on the center longitude -45°. Other X and Y lines get bended (not only rotated) against a degree grid:
The oblique mercator projection is rotated against true north in the center point. The mathematical formula rotates the local grid against the UTM grid, which does not align with true north in the center point. So the omerc projection can not give identical values as the formula.
To get a good precision, I took the given center point from the formulas as Corrego Allegre UTM 23S (EPSG:22523), and reprojected it to Corrego Allegre degrees (EPSG:4225) to get the center coordinates. With that, I created the following custom CRS:
+proj=omerc +y_0=8000 +x_0=6000 +alpha=-0.061 +lat_0=-17.219614491 +lonc=-46.91098576 +gamma=0 +ellps=intl +towgs84=-205.57,168.77,-4.12,0,0,0,0 +units=m +no_defs
+towgs84 is taken from Corrego Allegre, and y_0 and x_0 from the formulas.
From the formulas, I calculated a couple of points with LibreOffice Calc, and imported them as UTM 23S. The value +alpha is chosen empirically (for the reasons given above) to minimize the offset from the omerc grid to the points. It is about 0.4 meters for (10000 8000):
Note that the formula is only valid NorthEast of the center point. You see two of my points appearing at the "wrong" place.
Best Answer
The Helmert transformation that you're using is for converting between two geographic coordinate reference systems in 3D Cartesian, XYZ. It's not for applying a 2D Cartesian transformation to coordinates in planar / projected coordinate reference system.
I don't believe that PROJ.4 (or PROJ.5, just released in March 2018) has what you need.