4

I am trying to transform some data (both raster and vector) from WGS84 to a Croatian Coordinate system defined by:

+proj=tmerc +lat_0=0 +lon_0=15 +k=0.9999 +x_0=5500000 +y_0=0 +ellps=bessel +towgs84=550.499,164.116,475.142,5.80967,2.07902,-11.62386,0.99999445824 +units=m +no_defs

This works great in CS2CS (Rel 4.8.0):

cs2cs +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +to +proj=tmerc +lat_0=0 +lon_0=15 +k=0.9999 +x_0=5500000 +y_0=0 +ellps=bessel +towgs84=550.499,164.116,475.142,5.80967,2.07902,-11.62386,0.99999445824 +units=m +no_defs -r

And also great in OGR2OGR (version 1.9.2):

ogr2ogr -s_srs "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs" -t_srs "+proj=tmerc +lat_0=0 +lon_0=15 +k=0.9999 +x_0=5500000 +y_0=0 +ellps=bessel +towgs84=550.499,164.116,475.142,5.80967,2.07902,-11.62386,0.99999445824 +units=m +no_defs" out.shp  in.shp 

But in GDALWARP (version 1.9.2), the data comes out c.400m offset to the West (I think it's possibly a rotation error based upon it's appearance):

gdalwarp -s_srs "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs" -t_srs "+proj=tmerc +lat_0=0 +lon_0=15 +k=0.9999 +x_0=5500000 +y_0=0 +ellps=bessel +towgs84=550.499,164.116,475.142,5.80967,2.07902,-11.62386,0.99999445824 +units=m +no_defs" in.tif out.tif

An example of this can be seen below, where I digitised a vector layer around existing raster data, transformed both using the above commands, and then loaded into QGIS:

example of offset error

I can't for the life of me figure out why this could be - it is worthy of note that I have tried this with multiple Raster datasets from a variety of sources, and they always seem to be offset to the West by the same amount.

I have tried upgrading to GDAL 1.10.1, and this has no effect upon the output.

Has anyone ever seen this before and have any idea why it could be?

10
  • PROJ.4 and GDAL version? Nov 8, 2013 at 12:44
  • I'd guess at maybe a precision issue with the SRS parameters.
    – til_b
    Nov 8, 2013 at 12:50
  • Apologies, I've added versions into the question, it's GDAL/OGR 9.1.2, and Proj is Rel 4.8.0
    – jonnyhuck
    Nov 8, 2013 at 12:52
  • 1
    @jonnyhuck if x is the ppm value, the actual equations use 1 + x/1000000 to calculate the scale factor.
    – mkennedy
    Nov 11, 2013 at 15:21
  • 1
    @jonnyhuck I certainly may have calculated it wrong! 0.99999445824 = 1 + x/1000000, x = -5.54176
    – mkennedy
    Nov 11, 2013 at 17:10

1 Answer 1

1

After loading the files to QGIS, check the layers with Rightclick -> Set CRS for layer.

Sometimes QGIS guesses a similar existing EPSG code when you do not want it.

Testing on QGIS 1.8.0 Lisboa, it thinks that EPSG:31275 is ment, with proj string:

+proj=tmerc +lat_0=0 +lon_0=15 +k=0.9999 +x_0=5500000 +y_0=0 +ellps=bessel +towgs84=577.326,90.129,463.919,5.137,1.474,5.297,2.4232 +units=m +no_defs

which is alomst the same, except for the +towgs84 values.

So, nothing wrong with Gdalwarp ;-)

1
  • A cracking idea, but I'm afraid not - it's using the resulting "user generated" CRS. To be sure I created a new CRS using the PROJ string, and set both layers and the project to that, but I'm afraid it made no difference. Thanks for the suggestion though!
    – jonnyhuck
    Nov 8, 2013 at 15:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.