+towgs84 parameter ignored using gdalwarp

[jachym]

When using gdalwarp on GeoTIFF raster file and comparing the result with r.proj module in GRASS GIS, it seems so, that gdalwarp does not use the +towgs84 parameter:

1) Using gdalwarp
==================

$ gdalinfo QBsubset.TIF
Driver: GTiff/GeoTIFF
Size is 317, 317
Coordinate System is `'
Origin = (746558.400000000023283,5496931.200000000186265)
Pixel Size = (2.400000000000000,-2.400000000000000)
Corner Coordinates:
Upper Left  (  746558.400, 5496931.200)
Lower Left  (  746558.400, 5496170.400)
Upper Right (  747319.200, 5496931.200)
Lower Right (  747319.200, 5496170.400)
Center      (  746938.800, 5496550.800)
Band 1 Block=317x12 Type=UInt16, ColorInterp=Gray
Band 2 Block=317x12 Type=UInt16, ColorInterp=Undefined
Band 3 Block=317x12 Type=UInt16, ColorInterp=Undefined
Band 4 Block=317x12 Type=UInt16, ColorInterp=Undefined

$ gdalwarp -s_srs "+ellipsoid=WGS84 +proj=utm +zone=33" -t_srs "epsg:2065" QBsubset.TIF QBsubset-jtsk.tif

$ gdalinfo QBsubset-jtsk.tif
Driver: GTiff/GeoTIFF
Size is 355, 355
Coordinate System is:
PROJCS["S-JTSK (Ferro) / Krovak",
   GEOGCS["S-JTSK (Ferro)",
       DATUM["S_JTSK_Ferro",
           SPHEROID["Bessel 1841",6377397.155,299.1528128000008,
               AUTHORITY["EPSG","7004"]],
           AUTHORITY["EPSG","6818"]],
       PRIMEM["Ferro",-17.66666666666667],
       UNIT["degree",0.0174532925199433],
       AUTHORITY["EPSG","4818"]],
   UNIT["metre",1,
       AUTHORITY["EPSG","9001"]],
   AUTHORITY["EPSG","2065"]]
Origin = (-464078.173189933761023,-1131377.828972819726914)
Pixel Size = (2.398827733050111,-2.398827733050111)
Metadata:
 AREA_OR_POINT=Area
Corner Coordinates:
Upper Left  ( -464078.173,-1131377.829)
Lower Left  ( -464078.173,-1132229.413)
Upper Right ( -463226.589,-1131377.829)
Lower Right ( -463226.589,-1132229.413)
Center      ( -463652.381,-1131803.621)
Band 1 Block=355x11 Type=UInt16, ColorInterp=Gray
Band 2 Block=355x11 Type=UInt16, ColorInterp=Undefined
Band 3 Block=355x11 Type=UInt16, ColorInterp=Undefined
Band 4 Block=355x11 Type=UInt16, ColorInterp=Undefined

with

gdalwarp -s_srs "+ellipsoid=WGS84 +proj=utm +zone=33" -t_srs "+proj=krovak +ellipsoid=bessel +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56" QBsubset.TIF QBsubset-jtsk.tif

the result is just the same (IMHO should not be(?))

2) Using r.proj (GRASS GIS 6.3)
===============================
created location UTM zone 33:

GRASS6.3.cvs utm33/jachym: tmp >g.proj -w
PROJCS["UTM Zone 33, Northern Hemisphere",
   GEOGCS["wgs84",
       DATUM["unknown",
           SPHEROID["WGS_1984",6378137,298.257223563]],
       PRIMEM["Greenwich",0],
       UNIT["degree",0.0174532925199433]],
   PROJECTION["Transverse_Mercator"],
   PARAMETER["latitude_of_origin",0],
   PARAMETER["central_meridian",15],
   PARAMETER["scale_factor",0.9996],
   PARAMETER["false_easting",500000],
   PARAMETER["false_northing",0],
   UNIT["meter",1]]

imported the data

GRASS6.3.cvs utm33/jachym: tmp > g.region -g
n=5496931.2
s=5496170.4
w=746558.4
e=747319.2
nsres=2.4
ewres=2.4
rows=317
cols=317
cells=100489

getting bbox coordinates for new Krovak location:

GRASS6.3.cvs utm33/jachym: tmp > echo "746558.4 5496170.4"|cs2cs
+proj=utm +ellipsoid=WGS84 +zone=33 +to +proj=krovak +ellipsoid=bessel
+towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56
-463951.58      -1132011.59 -42.36

GRASS6.3.cvs utm33/jachym: tmp > echo "747319.2 5496931.2"|cs2cs
+proj=utm +ellipsoid=WGS84 +zone=33 +to +proj=krovak +ellipsoid=bessel
+towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56
-463099.38      -1131355.44 -42.34

creating new krovak location
coord system: krovak
specify datum: yes
datum: hermannskogel
transformation parameters: 3

GRASS6.3.cvs jtskutm/jachym: tmp >  g.region n=-1131355.44
s=-1132011.59 e=-463099.38 w=-463951.58

zooming out, the rasters are a bit rotated and can not fit to such region

GRASS6.3.cvs jtskutm/jachym: tmp > g.proj -w
PROJCS["Krovak",
   GEOGCS["bessel",
       DATUM["Militar_Geographische_Institut",
           SPHEROID["Bessel_1841",6377397.155,299.1528128],
           TOWGS84[570.8,85.7,462.8,4.998,1.587,5.261,3.56]],
       PRIMEM["Greenwich",0],
       UNIT["degree",0.0174532925199433]],
   PROJECTION["Krovak"],
   PARAMETER["latitude_of_center",0],
   PARAMETER["longitude_of_center",0],
   PARAMETER["azimuth",0],
   PARAMETER["pseudo_standard_parallel_1",0],
   PARAMETER["scale_factor",1],
   PARAMETER["false_easting",0],
   PARAMETER["false_northing",0],
   UNIT["meter",1]]

reprojecting the data - qb is imported QBsubset.TIF file

GRASS6.3.cvs jtskutm/jachym: tmp > r.proj in=qb.4 out=qb.4 location=utm33

and exporting

GRASS6.3.cvs jtskutm/jachym: tmp > r.out.gdal in=qb.4 out=qb4.tif type=Int16
ERROR 6: SetColorTable() only supported for Byte or UInt16 bands in
TIFF format. 100%
r.out.gdal complete.
GRASS6.3.cvs jtskutm/jachym: tmp > gdalinfo qb4.tif
Driver: GTiff/GeoTIFF
Size is 1987, 1531
Coordinate System is `'
Origin = (-464519.016431930009276,-1130918.535061219939962)
Pixel Size = (1.000234741781585,-1.000228658536903)
Metadata:
 AREA_OR_POINT=Area
Corner Coordinates:
Upper Left  ( -464519.016,-1130918.535)
Lower Left  ( -464519.016,-1132449.885)
Upper Right ( -462531.550,-1130918.535)
Lower Right ( -462531.550,-1132449.885)
Center      ( -463525.283,-1131684.210)
Band 1 Block=1987x2 Type=Int16, ColorInterp=Gray
 NoData Value=-32768
 Metadata:
   COLOR_TABLE_RULES_COUNT=1
   COLOR_TABLE_RULE_RGB_0=2.250000e+02 2.007000e+03 0 0 0 255 255 255

NOTE: the file does not contain datum transformation parameters, but this is realted to bug #1511

3) Final test
=============

qgis qb4.tif QBsubset.TIF

screenshot of this composition is available at ​http://les-ejk.cz/tmp/krovak-datumshift.png test data area available at ​http://les-ejk.cz/tmp/QBsubset.TIF

[jachym]

Original data

[warmerdam]

Jachym,

I don't know if it is related to your problem, but there is no +ellipsoid keyword. You are thinking +ellps butin this situation you should likely be using +datum=WGS84 instead.

Second, the TIFF format does not provide a mechanism to store TOWGS84 values. However, gdalwarp should use the towgs84 specified via -t_srs for this warp so that shouldn't matter.

I will try to reproduce the problem.

[jachym]

Replying to warmerdam:

I don't know if it is related to your problem, but there is no +ellipsoid keyword. You are thinking +ellps butin this situation you should likely be using +datum=WGS84 instead.

oh, thanks. I just knew, that some program from PROJ/GDAL family told me ones, that I used wrong parameter and because of no error occured, I thought, everything was all right

Second, the TIFF format does not provide a mechanism to store TOWGS84 values. However, gdalwarp should use the towgs84 specified via -t_srs for this warp so that shouldn't matter.

ok for me.

I will try to reproduce the problem.

Thanks a lot. I hope, ones will be the "krovak problem" really solved. We can get correct negative coordinates now, but the datum transformation parameters seems to be a bit messed up (cs2cs seems to work, gdalwarp not so well).

I tryed this: gdalwarp -s_srs "+datum=WGS84 +proj=utm +zone=33" -t_srs "+proj=krovak +ellps=bessel +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56" QBsubset.TIF QBsubset-jtsk.tif

result seems to be the same

[warmerdam]

Jachym,

I tried the following two warps, and they produce noticably different results, one with towgs84 and one without. I notice that your example with the +towgs84 parameter above did not include all the projection details.

gdalwarp -s_srs '+proj=utm +zone=33 +datum=WGS84' -t_srs '+proj=krovak +lat_0=49.5 +lon_0=42.5 +alpha=30.28813972222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +pm=ferro +units=m +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56 +no_defs' QBsubset.TIF out.tif

gdalwarp -s_srs '+proj=utm +zone=33 +datum=WGS84' -t_srs '+proj=krovak +lat_0=49.5 +lon_0=42.5 +alpha=30.28813972222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +pm=ferro +units=m +no_defs' QBsubset.TIF out.tif

[jachym]

Replying to warmerdam:

Jachym,

I tried the following two warps, and they produce noticably different results, one with towgs84 and one without. I notice that your example with the +towgs84 parameter above did not include all the projection details.

gdalwarp -s_srs '+proj=utm +zone=33 +datum=WGS84' -t_srs '+proj=krovak +lat_0=49.5 +lon_0=42.5 +alpha=30.28813972222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +pm=ferro +units=m +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56 +no_defs' QBsubset.TIF out.tif

gdalwarp -s_srs '+proj=utm +zone=33 +datum=WGS84' -t_srs '+proj=krovak +lat_0=49.5 +lon_0=42.5 +alpha=30.28813972222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +pm=ferro +units=m +no_defs' QBsubset.TIF out.tif

Are these parameters taken from epsg code 2065? It uses Ferro meridian and so the map is shifted to Ukraine. To make the map usable, I would need to run something like

I took parameter produced by g.proj and got also different results: all wrong (because of the East coordinate should be negative):

gdalwarp -s_srs '+proj=utm +zone=33 +datum=WGS84' -t_srs "`g.proj -wef`" QBsubset.TIF out-wef.tif

gdalwarp -s_srs '+proj=utm +zone=33 +datum=WGS84' -t_srs "`g.proj -jf`" QBsubset.TIF out-wef.tif

The results are different - I would expect, they are the same... and .. Easting is non-negative :-/

Anyway, I would expect, that gdalwarp will produce same result, as r.proj does

Is it a bug or feature? How to run set -s_srs so it produces map in S-JTSK coordinate system (+proj=krovak +ellps=bessel +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56), so that all coordinates in Czech and Slovak region are negative (and so correct)?

Thanks a lot and sorry, if I do not undertsand something right or if I'm messgin things together

[warmerdam]

Replying to jachym:

gdalwarp -s_srs '+proj=utm +zone=33 +datum=WGS84' -t_srs '+proj=krovak +lat_0=49.5 +lon_0=42.5 +alpha=30.28813972222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +pm=ferro +units=m +no_defs' QBsubset.TIF out.tif

Are these parameters taken from epsg code 2065? It uses Ferro meridian and so the map is shifted to Ukraine. To make the map usable, I would need to run something like

The above parameters are what GDAL/OGR expands EPSG:2065 too, which I noticed you used once earlier on. I just focused on the point you raised that the towgs84 parameter was having no effect, and tried to test this. With the above settings the towgs84 was having an effect.

I took parameter produced by g.proj and got also different results: all wrong (because of the East coordinate should be negative):

gdalwarp -s_srs '+proj=utm +zone=33 +datum=WGS84' -t_srs "`g.proj -wef`" QBsubset.TIF out-wef.tif

gdalwarp -s_srs '+proj=utm +zone=33 +datum=WGS84' -t_srs "`g.proj -jf`" QBsubset.TIF out-wef.tif

The results are different - I would expect, they are the same... and .. Easting is non-negative :-/

As presented above, I have no idea what parameters actually got used.

Anyway, I would expect, that gdalwarp will produce same result, as r.proj does

If we can get to the point that r.proj and gdalwarp are using the same underlying PROJ.4 definitions then I would expect the results to be similar.

Is it a bug or feature? How to run set -s_srs so it produces map in S-JTSK coordinate system (+proj=krovak +ellps=bessel +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56), so that all coordinates in Czech and Slovak region are negative (and so correct)?

I'm concerned that we have wandered through several related questions, but I have trouble unless I focus on one thing.

Are you satisfied that towgs84 parameters are working fine? If so, and if you feel that something else isn't working properly, it would likely be best to close this ticket, and start a new one focused on the new issue.

[jachym]

Replying to warmerdam:

Are you satisfied that towgs84 parameters are working fine? If so, and if you feel that something else isn't working properly, it would likely be best to close this ticket, and start a new one focused on the new issue.

yes, I can confirm. IMHO it was caused by the "+ellipsoid" parameters

{{{gdalwarp -s_srs '+proj=utm +zone=33 +datum=WGS84' -t_srs '+proj=krovak +ellps=bessel +units=m' QBsubset.TIF jtsk2.tif }}} and gdalwarp -s_srs '+proj=utm +zone=33 +datum=WGS84' -t_srs '+proj=krovak +ellps=bessel +units=m +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56 +no_defs' QBsubset.TIF jtsk1.tif

are working nice, like expected (the output is on wrong place though - but that is another issue).

I'm closing this bug, thanks for explanation and help.

[jachym]

however:

gdalwarp -sgdalwarp -s_srs '+proj=utm +zone=33 +datum=WGS84' -t_srs "+init=epsg:2065" QBsubset.TIF out1.tif

and

gdalwarp -s_srs '+proj=utm +zone=33 +datum=WGS84' -t_srs "+init=epsg:2065 +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56" QBsubset.TIF out1.tif

are producing the same result:

gdalinfo out1.tif
Driver: GTiff/GeoTIFF
Size is 374, 374
Coordinate System is `'
Origin = (1317941.395995946368203,-987858.858327466063201)
Pixel Size = (2.401169904767594,-2.401169904767594)
Metadata:
  AREA_OR_POINT=Area
Corner Coordinates:
Upper Left  ( 1317941.396, -987858.858)
Lower Left  ( 1317941.396, -988756.896)
Upper Right ( 1318839.434, -987858.858)
Lower Right ( 1318839.434, -988756.896)
Center      ( 1318390.415, -988307.877)
Band 1 Block=374x10 Type=UInt16, ColorInterp=Gray
Band 2 Block=374x10 Type=UInt16, ColorInterp=Undefined
Band 3 Block=374x10 Type=UInt16, ColorInterp=Undefined
Band 4 Block=374x10 Type=UInt16, ColorInterp=Undefined

Seems like if GDAL could not understand the -t_srs parameters, it just transforms the map in some strange system..

BTW, trying this:

gdalwarp -s_srs '+proj=utm +zone=33 +datum=WGS84' QBsubset.TIF jtsk4.tif

So running gdalwarp withou any -t_srs parameters produces the same result, like using above mentioned parameters.

The same result was obtained, when using the "+ellipsoid" parameter.

Is this a bug or feature?

[Markus Neteler]

I have the same problem (GDAL 1.5.1):

Attempt to fix missing proj information in GeoTIFF:

# note: output manually broken for trac convenience:
cat /geodata2_originals_raid5/gaussboaga1.prj
PROJCS["Transverse Mercator",GEOGCS["international",DATUM["Monte_Mario",SPHEROID["International_1924",6378388,297],
TOWGS84[-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68]],PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",9],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",1500000],
PARAMETER["false_northing",0],UNIT["Meter",1]]


gdal_translate -a_srs /geodata2_originals_raid5/gaussboaga1.prj originals/102010.tif 102010.tif

gdalinfo 024120.tif
Driver: GTiff/GeoTIFF
Files: 024120.tif
Size is 13200, 11600
Coordinate System is:
PROJCS["Transverse Mercator",
    GEOGCS["international",
        DATUM["Monte_Mario",
            SPHEROID["International 1924",6378388,297.000000000005,
                AUTHORITY["EPSG","7022"]],
            AUTHORITY["EPSG","6265"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",9],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",1500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (1621400.000000000000000,5150920.000000000000000)
Pixel Size = (0.500000000000000,-0.500000000000000)
Metadata:
  AREA_OR_POINT=Area
...

The towgs84 part is omitted. No apparent way to include it (I also tried the full PROJ4 string on command line).

Could it be the same problem as reported in bug #2047?

Markus

[Markus Neteler]

Just learned that GeoTIFF has no mechanism to store TOWGS84 parameters.