Opened 15 years ago

Closed 11 years ago

#601 closed defect (fixed)

[PATCH] wrong handling of units in geotiff

Reported by: pierrick.brihaye@… Owned by: warmerdam
Priority: normal Milestone: 1.5.3
Component: GDAL_Raster Version: 1.4.2
Severity: minor Keywords: geotiff grad
Cc: pierrick.brihaye@…, favard@…

Description (last modified by warmerdam)

Hi,

Trying to assign a french SRS to a tif.file :

bin\gdal_translate -a_srs EPSG:27572 29_test.tif 29_test_geo.tif

The gdalinfo result is :
Driver: GTiff/GeoTIFF
Size is 4000, 4000
Coordinate System is:

PROJCS["NTF (Paris) / Lambert zone II",
    GEOGCS["NTF (Paris)",
        DATUM["unknown",
            SPHEROID["unretrievable - using WGS84",6378137,298.257223563,
                AUTHORITY["EPSG","0"]],
            AUTHORITY["EPSG","6807"]],
        PRIMEM["Paris",2.596921300000003],
        UNIT["grad",0.01570796326794895],
        AUTHORITY["EPSG","4807"]],
    PROJECTION["Lambert_Conformal_Conic_1SP"],
    PARAMETER["latitude_of_origin",57.77777777777791],
    PARAMETER["central_meridian",-2.885468111111118],
    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)
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=Gray

Why these values for DATUM and SPHEROID ? It should be :

DATUM["Nouvelle Triangulation Francaise (Paris)",
      SPHEROID["Clarke 1880 (IGN)", 6378249.2, 293.4660212936269,
               AUTHORITY["EPSG","7011"]],
      AUTHORITY["EPSG","6807"]
],

Is it a problem with pcs.csv or with gdalinfo's interpretation of the SRS ?

Note : used today's http://qgis.org/GDALBuild.zip compiled distribution for Windows.

Cheers,

p.b.

Attachments (6)

gdal-1.5.0-bug-gt_wkt_srs.patch (1.1 KB) - added by dgrichard 11 years ago.
proposed patch for fixing this bug
test_gf.tif (810 bytes) - added by dgrichard 11 years ago.
Small 10x10 GTiff image to check the bug fix
test_gf.gtf (1.1 KB) - added by dgrichard 11 years ago.
Metadata on the test image
bug-601-gdalinfo-without-patch.txt (1.6 KB) - added by dgrichard 11 years ago.
gdalinfo on the test image without the patch applied
bug-601-gdalinfo-with-patch.txt (1.6 KB) - added by dgrichard 11 years ago.
gdalinfo on the same test image once patch applied
gdal-1.6.0-bug-gt_wkt_srs.patch (1.3 KB) - added by dgrichard 11 years ago.
patch against trunk

Download all attachments as: .zip

Change History (21)

comment:1 Changed 15 years ago by warmerdam

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.

comment:2 Changed 15 years ago by pierrick.brihaye@…

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.



comment:3 Changed 15 years ago by pierrick.brihaye@…

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.

comment:4 Changed 15 years ago by pierrick.brihaye@…

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.



comment:5 Changed 15 years ago by warmerdam

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. 



comment:6 Changed 15 years ago by pierrick.brihaye@…

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.

comment:7 Changed 14 years ago by pierrick.brihaye@…

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,


comment:8 Changed 14 years ago by pierrick.brihaye@…

Same problem with :

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

comment:9 Changed 14 years ago by pierrick.brihaye@…

Same sign error (?) with :

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

comment:10 Changed 13 years ago by mapserver@…

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.

comment:11 Changed 13 years ago by favard@…

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

comment:12 Changed 12 years ago by martinoty

Component: OGR_SRSGDAL_Raster
Priority: lownormal
Summary: pcs.csv accuracy ?wrong handling of units in geotiff
Version: unspecified1.4.2

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.

Changed 11 years ago by dgrichard

proposed patch for fixing this bug

Changed 11 years ago by dgrichard

Attachment: test_gf.tif added

Small 10x10 GTiff image to check the bug fix

Changed 11 years ago by dgrichard

Attachment: test_gf.gtf added

Metadata on the test image

Changed 11 years ago by dgrichard

gdalinfo on the test image without the patch applied

Changed 11 years ago by dgrichard

gdalinfo on the same test image once patch applied

comment:13 Changed 11 years ago by Even Rouault

Keywords: geotiff grad added
Milestone: 1.5.1

Changed 11 years ago by dgrichard

patch against trunk

comment:14 Changed 11 years ago by Even Rouault

Summary: wrong handling of units in geotiff[PATCH] wrong handling of units in geotiff

comment:15 Changed 11 years ago by warmerdam

Description: modified (diff)
Resolution: fixed
Status: assignedclosed

Didier,

Sorry for letting this languish. I have applied the portion of the patch that related to angular units. But the portion that related to linear units seems questionable. It certainly broke the autotest/gcore/gtiff_write.py test script:

Script: gcore/gtiff_write.py
  TEST: SetProjection: byte.tif ... fail
    Did not get expected projection reference.

Basically it reversed the conversion of feet to meters and I'm not convinced it should have. So I'd like you to file a new ticket specifically addressing linear units in GeoTIFF if you are convinced that it is still broken.

The reduced patch was applied in trunk (r15501) and 1.5 branch (r15502). I also added a test suite item for the angular units issue using your test file (r15503).

Note: See TracTickets for help on using tickets.