Opened 20 years ago
Closed 15 years ago
#601 closed defect (fixed)
[PATCH] wrong handling of units in geotiff
Reported by: | 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 )
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)
Change History (21)
comment:2 by , 20 years ago
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 by , 20 years ago
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 by , 20 years ago
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 by , 20 years ago
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 by , 20 years ago
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 by , 19 years ago
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:9 by , 18 years ago
Same sign error (?) with : GDAL 1.3.1.0, FWTools 1.0.0a5, released 2005/10/21
comment:10 by , 18 years ago
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 by , 17 years ago
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 by , 17 years ago
Component: | OGR_SRS → GDAL_Raster |
---|---|
Priority: | low → normal |
Summary: | pcs.csv accuracy ? → wrong handling of units in geotiff |
Version: | unspecified → 1.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.
by , 16 years ago
Attachment: | gdal-1.5.0-bug-gt_wkt_srs.patch added |
---|
proposed patch for fixing this bug
by , 16 years ago
Attachment: | bug-601-gdalinfo-without-patch.txt added |
---|
gdalinfo on the test image without the patch applied
by , 16 years ago
Attachment: | bug-601-gdalinfo-with-patch.txt added |
---|
gdalinfo on the same test image once patch applied
comment:13 by , 16 years ago
Keywords: | geotiff grad added |
---|---|
Milestone: | → 1.5.1 |
comment:14 by , 16 years ago
Summary: | wrong handling of units in geotiff → [PATCH] wrong handling of units in geotiff |
---|
comment:15 by , 15 years ago
Description: | modified (diff) |
---|---|
Resolution: | → fixed |
Status: | assigned → closed |
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).