[PATCH] wrong handling of units in geotiff

Pierrick,

I tried this quickly with:

gdal_translate -a_srs EPSG:27572 openev/utm.tif out.tif

And got:

PROJCS["NTF (Paris) / Lambert zone II",
    GEOGCS["NTF (Paris)",
        DATUM["Nouvelle_Triangulation_Francaise_Paris",
            SPHEROID["Clarke 1880 (IGN)",6378249.2,293.4660212936269,
                AUTHORITY["EPSG","7011"]],
            AUTHORITY["EPSG","6807"]],
        PRIMEM["Paris",2.596921300000003],
        UNIT["grad",0.01570796326794895],
        AUTHORITY["EPSG","4807"]],
    PROJECTION["Lambert_Conformal_Conic_1SP"],
    PARAMETER["latitude_of_origin",57.7777777777779],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",0.99987742],
    PARAMETER["false_easting",600000],
    PARAMETER["false_northing",2200000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","27572"]]

on the resulting file. 

So I think the issue is either related to the version of the GDAL code, 
the way libgeotiff is integrated in your build, or whether all the right
files are being found.  

Could you attach a small example file you have produced with gdal_translate?
I would like to inspect the GeoTIFF file and see if it is correct coming
out of gdal_translate.  If so, then the issue is likely related to having
all the correct support files being found for gdalinfo.

Hi again,

Thanks for your answer. I will drive further tests tomorrow (5:30 PM here :-)...

But : see
http://docs.codehaus.org/display/GEOTOOLS/Installing+the+EPSG+database+in+PostgreSQL

It looks like PARAMETER["latitude_of_origin",57.77777777777791] is incorrect.

Well, working on Windows is not that easy :

I'm now working with OpenEV binaries that look more recent than the ones
provided with QGIS. 

However, OpenEV binaries don't include the "data" directory : that's why I'm now
working with the one provided with the 1.2.1 version at
http://dl.maptools.org/dl/gdal/.

For now, I still get the same results with a slighly more verbose output :

Corner Coordinates:
Upper Left  (   40000.000, 2420000.000) (  9d26'8.23"W, 59d38'50.79"N)
Lower Left  (   40000.000, 2410000.000) (  9d24'57.85"W, 59d32'53.95"N)
Upper Right (   50000.000, 2420000.000) (  9d16'7.01"W, 59d39'32.29"N)
Lower Right (   50000.000, 2410000.000) (  9d14'57.86"W, 59d33'35.36"N)
Center      (   45000.000, 2415000.000) (  9d20'32.74"W, 59d36'13.19"N)

Well, it should be around :  5°08'29" W / 48°26'58" N.

Anyway, stay tuned.

p.b.

Hi,

Test files available at http://perso.wanadoo.fr/pierrick.brihaye/tests.zip

FYI :
bin\gdal_translate --version
GDAL 1.2.0.0, released 2004/03/10
Usage: gdal_translate [--version]
       [-ot {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/
             CInt16/CInt32/CFloat32/CFloat64}] [-not_strict]
       [-of format] [-b band] [-outsize xsize[%] ysize[%]]
       [-scale [src_min src_max [dst_min dst_max]]]
       [-srcwin xoff yoff xsize ysize] [-a_srs srs_def]
       [-projwin ulx uly lrx lry] [-co "NAME=VALUE"]*
       [-gcp pixel line easting northing]*
       [-mo "META-TAG=VALUE"]* [-quiet]
       src_dataset dst_dataset

GDAL 1.2.0.0, released 2004/03/10

The following format drivers are configured and support output:
  VRT: Virtual Raster
  GTiff: GeoTIFF
  NITF: National Imagery Transmission Format
  HFA: Erdas Imagine Images (.img)
  ELAS: ELAS
  AAIGrid: Arc/Info ASCII Grid
  DTED: DTED Elevation Raster
  PNG: Portable Network Graphics
  MEM: In Memory Raster
  XPM: X11 PixMap Format
  BMP: MS Windows Device Independent Bitmap
  PCIDSK: PCIDSK Database File
  PNM: Portable Pixmap Format (netpbm)
  ENVI: ENVI .hdr Labelled
  EHdr: ESRI .hdr Labelled
  PAux: PCI .aux Labelled
  MFF: Atlantis MFF Raster
  MFF2: Atlantis MFF2 (HKV) Raster
  BT: VTP .bt (Binary Terrain) 1.3 Format
  FIT: FIT Image

It looks like the --version switch doesn't work properly...

Cheers,

p.b.

Hi,

Now working with the tools provided by the latest OpenEV 1.8.1 Windows release.
Their behaviour looks more compliant with the documentation :-) and I get the
same results than you.

However, I confirm that there is a problem with :

PARAMETER["latitude_of_origin",57.77777777777791],

It should be 52 *grades* (see : UNIT["grad",0.01570796326794895]). 

It is obvious here that 52 grades have been over-multiplicated by the
grade/deegree ratio which is equal to 400 / 360 or 1.111...

The problem looks to be the same with the longitude with an extra shift problem
; where does PARAMETER["central_meridian",-2.885468111111118] come from ?

Cheers,

p.b.

Pierrick,

OGR deliberately converts angular values from non degree units into degrees.
So, I don't feel that the translation of grads to degrees is wrong.  Is your
concern that the GEOGCS still lists grads as the angular unit and you feel 
this should apply to the angular parameters of the PROJCS declaration? 

Hmm, reading over my own document on WKT and how various items should be
interpreted I see that I suggest that if the GEOGCS is in gradians then 
the angular projection parameters have to be as well. 

  http://192.168.2.2/~warmerda/projects/opengis/wktproblems.html

OK ... so that is a bug I think. 

On the central meridian issue, I tried running testepsg.exe using 
OpenEV_FW 1.8.1 and got a central meridian of 0.  Note that I don't
think the OpenEV_FW 1.8.1 release included testepsg.exe, so I built mine
locally but used something that is essentially (but not exactly) the 
OpenEV_FW 1.8.1 release.  So, I am still unable to reproduce this problem. 

Hi again Frank,

Thank your for your feedback. I appreciate it very much.

>OGR deliberately converts angular values from non degree units into degrees.

No problem. However the csv file indicates 52.0 and I am *sure* we are speaking
about grades because this figure comes from what we can't bargain with : french
bureaucracy :-)

57.7777777777779 is the result of 52 grades * (400 grades/360 degrees) : we
definitely don't have degrees there IMHO.

Anyway, even if I find this PROJCS very strange (and buggy), my actual problem
is this bounding box :

Corner Coordinates:
Upper Left  (   40000.000, 2420000.000) (  9d26'8.23"W, 59d38'50.79"N)
Lower Left  (   40000.000, 2410000.000) (  9d24'57.85"W, 59d32'53.95"N)
Upper Right (   50000.000, 2420000.000) (  9d16'7.01"W, 59d39'32.29"N)
Lower Right (   50000.000, 2410000.000) (  9d14'57.86"W, 59d33'35.36"N)
Center      (   45000.000, 2415000.000) (  9d20'32.74"W, 59d36'13.19"N)

... a very particular view of interceltism :-) since the correct location is here :
http://gnswww.nga.mil/geonames/Gazetteer/Search/Results.jsp?Feature__Unique_Feature_ID=-1456387&Diacritics=Yes&reload=1

Regarding testepsg.exe, I have to check at my office if it is distributed with
OpenEV. But what is it used for ? I can't find any docs...

Furthermore, http://192.168.2.2/~warmerda/projects/opengis/wktproblems.html
can't be resolved.

Cheers,

p.b.

Bug still present in FWTools0.9.8.

C:\FWTools0.9.8>gdalinfo --version
GDAL 1.2.6.0, released 2005/03/13

C:\FWTools0.9.8>gdalinfo F004_026.TIF
Driver: GTiff/GeoTIFF
Size is 4000, 4000
Coordinate System is `'
Origin = (40000.000000,2420000.000000)
Pixel Size = (2.50000000,-2.50000000)
Metadata:
  TIFFTAG_SOFTWARE=Handmade Software, Inc. Image Alchemy v1.7.6
Corner Coordinates:
Upper Left  (   40000.000, 2420000.000) 
Lower Left  (   40000.000, 2410000.000) 
Upper Right (   50000.000, 2420000.000) 
Lower Right (   50000.000, 2410000.000) 
Center      (   45000.000, 2415000.000) 
Band 1 Block=4000x2 Type=Byte, ColorInterp=Palette
  Color Table (RGB with 256 entries)
    0: 0,0,0,255
    1: 50,0,0,255
    2: 100,0,0,255
    3: 150,0,0,255
[snip]

C:\FWTools0.9.8>gdal_translate --version
GDAL 1.2.6.0, released 2005/03/13

C:\FWTools0.9.8>gdal_translate -a_srs EPSG:27572 F004_026.TIF georeferenced.tif

Input file size is 4000, 4000
0...10...20...30...40...50...60...70...80...90...100 - done.

C:\FWTools0.9.8>gdalinfo georeferenced.tif

Driver: GTiff/GeoTIFF
Size is 4000, 4000
Coordinate System is:
PROJCS["NTF (Paris) / Lambert zone II",
    GEOGCS["NTF (Paris)",
        DATUM["Nouvelle_Triangulation_Francaise_Paris",
            SPHEROID["Clarke 1880 (IGN)",6378249.2,293.4660212936265,
                AUTHORITY["EPSG","7011"]],
            AUTHORITY["EPSG","6807"]],
        PRIMEM["Paris",2.33722917],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4807"]],
    PROJECTION["Lambert_Conformal_Conic_1SP"],
    PARAMETER["latitude_of_origin",46.8],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",0.99987742],
    PARAMETER["false_easting",600000],
    PARAMETER["false_northing",2200000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","27572"]]
Origin = (40000.000000,2420000.000000)
Pixel Size = (2.50000000,-2.50000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_SOFTWARE=Handmade Software, Inc. Image Alchemy v1.7.6
Corner Coordinates:
Upper Left  (   40000.000, 2420000.000) (  9d55'44.41"W, 48d32'6.79"N)
Lower Left  (   40000.000, 2410000.000) (  9d54'57.50"W, 48d26'44.66"N)
Upper Right (   50000.000, 2420000.000) (  9d47'39.32"W, 48d32'37.73"N)
Lower Right (   50000.000, 2410000.000) (  9d46'53.23"W, 48d27'15.54"N)
Center      (   45000.000, 2415000.000) (  9d51'18.62"W, 48d29'41.25"N)
Band 1 Block=4000x2 Type=Byte, ColorInterp=Palette
  Color Table (RGB with 256 entries)
    0: 0,0,0,255
    1: 50,0,0,255
    2: 100,0,0,255
    3: 150,0,0,255
[snip]

The PROJCS seems to be OK (no more grad/degree problems ; see below however) but
the Corner's longitudes are still wrong.

I expect them around 5d08'29" W, which is 5.1413 W (decimal). Notice that 5.1413
+ 2.3372 + 2.3372 = 9,8157, i.e. 9d48'56" W, inside the returned box.

It looks like there is a double shift for the prime meridian.

Cheers,

Same problem with :

GDAL 1.3.1.0, FWTools 1.0.0a3, released 2005/10/07

Same sign error (?) with :

GDAL 1.3.1.0, FWTools 1.0.0a5, released 2005/10/21

I encounter a problem in the reprojection of EPSG 27572 to EPSG 4326.
My reprojected layer are nearly 400 km on the left of the other layer.
I use gdal 1.3.1 and proj 4.4.9

When i manually change the original epsg entry 
(
<27572> +proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=-2.33722917
+k_0=0.99987742 +x_0=600000 +y_0=2200000 +a=6378249.2 +b=63565
15 +towgs84=-168,-60,320,0,0,0,0 +pm=paris +units=m +no_defs  <>
)
which comes with proj 4.4.9 for EPSG 27572 to:

<27572> +proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=2.33722917
+k_0=0.99987742 +x_0=600000 +y_0=2200000 +a=6378249.2 +b=6356515
+towgs84=-168,-60,320,0,0,0,0 +units=m +no_defs  <>

or

<27572> +proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 +k_0=0.99987742
+x_0=600000 +y_0=2200000 +a=6378249.2 +b=6356515
+towgs84=-168,-60,320,0,0,0,0 +pm=paris +units=m +no_defs  <>

the shift disappear.

I have no idea what happens behind the curtain but i assume that it has
something to do with the prime meridian. My be a sign error.

perhaps the bug is in the gt_wkt_srs.cpp file
in the function GTIFSetFromOGISDefn when Lambert_Conformal_Conic_1SP projection is recognized, ProjCoordTransGeoKey key is set with CT_LambertConfConic_2SP value. I think it should be set with CT_LambertConfConic_1SP value. Could you confirm me that it is the bug...
thank you

I think three bugs are mixed here.

The first one is related to a coordinate shift when prime meridians are used: it's bug #367.

The second one is related to a problem when writing tags in geotiff with Lambert_Conformal_Conic_1SP : it's bug #1315.

The third one is related to units' handling in geotiff: projection parameters are wrongly divided by (instead of multiplied by) the conversion ratio to degrees. This explains why one get 57.7777777777779° instead of 46.8°, the right value.

I think the correction is the following one:

In frmts/gtiff/gt_wkt_srs.cpp#L381 , replace:

adfParm[0] /= psDefn->UOMAngleInDegrees; adfParm[1] /= psDefn->UOMAngleInDegrees; adfParm[2] /= psDefn->UOMAngleInDegrees; adfParm[3] /= psDefn->UOMAngleInDegrees;

adfParm[5] /= psDefn->UOMLengthInMeters; adfParm[6] /= psDefn->UOMLengthInMeters;

with:

adfParm[0] *= psDefn->UOMAngleInDegrees; adfParm[1] *= psDefn->UOMAngleInDegrees; adfParm[2] *= psDefn->UOMAngleInDegrees; adfParm[3] *= psDefn->UOMAngleInDegrees;

adfParm[5] *= psDefn->UOMLengthInMeters; adfParm[6] *= psDefn->UOMLengthInMeters;

As far as I know, this only affects gdalinfo, but is not a problem for actual projection conversions.

proposed patch for fixing this bug

Small 10x10 GTiff image to check the bug fix

Metadata on the test image

gdalinfo on the test image without the patch applied

gdalinfo on the same test image once patch applied

patch against trunk