Reprojection with ogr2ogr works fine on command line, but shifts the polygon to northeast when same command is run from inside python using os.system(...)

Any ideas why this could be?

I hard-coded the full ogr2ogr path to avoid hitting different installations from each environment.

PRJ files for both versions are identical so the difference must be in the "shp" file.

e.g. Here is the command:

/opt/local/bin/ogr2ogr -s_srs epsg:28992 -t_srs epsg:4326 -a_srs epsg:4326 /mnt/output_4326.shp /mnt/input_28992.shp
  • What is your Python command ? – gene Aug 29 '14 at 14:30
  • The command is exactly the same when run in OSX Terminal and using os.system(). i.e. os.system("/opt/local/bin/ogr2ogr -s_srs epsg:28992 -t_srs epsg:4326 -a_srs epsg:4326 /mnt/output_4326.shp /mnt/input_28992.shp") – python1981 Aug 29 '14 at 14:31
  • At least using -a_srs epsg:4326 is not needed when you have -t_srs. Take it away even it does not probably chance anything in your case. But think about situation with -t_srs epsg:4326 -a_srs epsg:3857. If the command works it would transform the data into epsg:4326 but write into metadata that the product is in epsg:3857. – user30184 Aug 29 '14 at 14:36
  • 1
    @user30184 The proj4 string of sr.org misses +towgs84 parameters. This is from current QGIS: +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.417,50.3319,465.552,-0.398957,0.343988,-1.8774,4.0725 +units=m +no_defs – AndreJ Aug 29 '14 at 14:43
  • 1
    Thus my first suspect is that ogr2ogr and python code are finding different Proj4 versions. For controlling what really happens use -s_srs "+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.417,50.3319,465.552,-0.398957,0.343988,-1.8774,4.0725 +units=m +no_defs" – user30184 Aug 29 '14 at 14:56

With GDAL and other Proj4 based programs it is handy to define the projection by its EPSG code shortcut like "epsg:28992". Shortcut is converted to actual Proj4 projection string by doing a lookup into proj installation files. Unfortunately there are lots of projections with several different projection definitions hanging aroud. Most typical case is that +towgs84 is either missing totally or different parameter sets have been delivered with different Proj4 versions.

For user it can be difficult to control what really happens when epsg code is user. There is a human readable "epsg" file which is stored, for example, to \proj\SHARE directory. It is possible to edit that file manually and fix the +proj string. Such fix works for example with Mapserver which really reads the string from this file. However, many other programs are generating the +proj string by utilizing other, more native files which are present in the proj installation directories. Actually, the epsg text file is a product from running the native methods for all the supported projections. In practise that means that editing the epsg text does not have any effect for QGIS not GDAL utilities like ogr2ogr.

The most safe way for the user is not to rely on the lookup system utilising the EPSG codes. By defining the projection directly with +proj string, in this case as "+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.417,50.3319,465.552,-0.398957,0.343988,-1.8774,4.0725 +units=m +no_defs" is explicit and will seldom lead to surprises caused by different Proj4 versions which are installed on the system. Some possibility to error still remains: it is possible that the mathematics behind for example "sterea" is different in different Proj4 versions because of a bug fix or other improvement.

When taking the explicit way it is better to use it on both sides of the conversion. Especially if +towgs84 is added, the target projection must have it defined as well either by +towgs84, +datum, or +nadgrids parameters because datum shift happens only if it is defined on both sides. For examle WGS84-UTM zones have it defined in current Proj4 versions as

# WGS 84 / UTM zone 32N
<32632> +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs  <>

whereas European GRS80 based UTM zones use +towgs84 parameters. Notice the use of a transform string "0,0,0,0,0,0,0" which seems to do nothing but it must still be there.

# ETRS89 / UTM zone 32N
<25832> +proj=utm +zone=32 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs  <>
  • Just another thing you missed: It is necessary to add towgs84=0,0,0,0,0,0,0 to projections based on WGS84. proj does the datum shift only if +towgs84 (or +datum or +nadgrids) is given on both sides. – AndreJ Aug 30 '14 at 18:06
  • Thanks, I edited my answer. I am not sure if +nadgrids alone is enough. I have thought that nadgrids just fine tune the 7-parameter transformation. – user30184 Aug 30 '14 at 18:54
  • No, it replaces them. – AndreJ Aug 30 '14 at 19:38
  • I rather lot trust you but I still made a small test. You are right, +nadgrids replaces +towgs84. – user30184 Aug 30 '14 at 20:12
  • The same technique also worked nicely to fix similar issues with cs2cs – python1981 Sep 3 '14 at 11:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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