GeoTIFF Linear Units Wrong

[warmerdam]

Section 2.7.2 of the GeoTIFF specification makes it clear that the units of ProjFalseEastingGeoKey? and ProjFalseNorthingGeoKey? should be the value of ProjLinearUnit?. GDAL assumes these are always in meters when reading or writing.

Grr.

[gaopeng]

Frank,
Have you had chance to look at this?

Here is one approach that may work:
In GDAL 1.8:
- make a change in tif to remove the conversion of feet-to-meter and store easting/northing value as it is.
- store a version in tiff as, e.g. metadata, or citation key.

In GDAL 1.8, 1.6, 1.4, and etc:
- perform conversion on read only for existing tiff, i.e. if the version is not present.

I think we should fix this sooner so GDAL will stop creating bad tiffs.

[warmerdam]

Gao,

As I read your suggestion we would assume that any GeoTIFF not marked with a special marker would be assumed to be broken. I can't really pursue this approach in good conscience. It flies in the face of proper standards support.

I have discovered the problem with unit interpretation is two fold. First, libgeotiff's GTIFGetDefn() function was not properly normalizing things to meters as it claimed if the parameters were coming from geokeys as opposed to from epsg. Second, on write, the GDAL gt_wkt_srs.cpp code to write the geokeys was writing them in normalized (meters) form when it shouldn't have. I have corrected these issues in trunk (r22310) along with some tests for this scenario. There was also an upstream fix in libgeotiff (r1998).

Next I will try to construct a mechanism so those needing to can (temporarily) fallback to the previous behavior and look a porting into 1.8 and 1.8-esri branches.

[gaopeng]

Frank, My understanding is the fix will write correct tiffs, i.e. not normalize (to meters) when it shouldn't have. After the fix, there will be 3 types of tiffs: (1) old tiffs, with wrong unit, (2) new tiffs with correct unit, and (3) 3rd party tiffs, with correct unit. How do different versions of GDAL, before and after the fix handle these different tiffs? Is the following matrix correct:

           work with      doesnt work with

Before      1             2, 3
After       2, 3          1

[warmerdam]

Gao, your compatability matrix is accurate.

In trunk I have added logic to try and fix-up the linear units if we are using an older libgeotiff. Also added a configuration option, GTIFF_LINEAR_UNITS=BROKEN to read files with the wrong linear units interpretation (r22311).

It occurs to me that we can often identify if a GeoTIFF file was generated by GDAL based on the existance of GDAL specific metadata tags. If we used that in combination with your version id, it might be possible to mostly identify old broken GDAL generated files and handle them appropriately. I'm going to investigate this possibility.

[warmerdam]

I have reviewed a number of GeoTIFF files produced with GDAL and none had the gdal metadata tag which has brought me to realize that this is at best an occasional way of identifying old GDAL files. So I have not pursued this avenue.

I have back ported the changes so far to 1.8 branch (r22338) and 1.8-esri branch (r22339).

I'm not sure what additional steps would be appropriate at this time.

[gaopeng]

Frank, The current fix will create a situation where ArcGIS 10.1 will generate (new) tiff images, that early releases of ArcGIS can't read (can read but with wrong SR). This will cause problems. Could we make it work as the following by (1) introducing a version flag in the current GDAL version, together with the fix; (2) patching early GDAL versions to check the version info and apply unit conversion properly:

                             work with      doesnt work with

GDAL 1.4-esri, 1.6-esri      1, 2           3
GDAL 1.8-esri                1, 2           3

(1) Store a version# in newly created tiff

(2) Add logic to GDAL 1.4-esri, 1.6-esri, and 1.8-esri to check the version#. If present (and as expected), use the new logic for linear units; otherwise, assume it's an old tiff, and continue to use the old logic.

[warmerdam]

I have added a new (gdal specific) geokey "ProjLinearUnitsInterpCorrectGeoKey?" (3059) to mark files that are written by GDAL with non-meter linear units using the correct interpretation. If GTIFF_LINEAR_UNITS is set to BROKEN then these files are *not* assumed to be broken and adjusted.

This is in trunk (r22360, r22361, r22362) and 1.8 branch (r22365).

A similar change is in 1.8-esri but now 1.8-esri assumes GTIFF_LINEAR_UNITS BROKEN. So that means any files without this special geokey will be assumed to be old GDAL produced files, and fixed up accordingly. This means that files from other sources written properly will be incorrectly interpreted. (r22364).

Investigating appropriate backporting to 1.4-esri and 1.6-esri...

[warmerdam]

I have migrated the 1.8-esri change back to 1.6-esri (r22397) and 1.4-esri (r22398).

[rprinceley]

Frank,

We have hit issues with certain spatial references that have user defined linear units. In this case, when setting projection on oSRS object in function GTIFGetOGISDefn the false easting/northing values are converted to meters (via SetNormProjParm?).

E.g. UNIT["Foot_Clarke", 0.304797265]

[warmerdam]

Robin,

Can you provide more details on how to reproduce the problem?

[rprinceley]

Frank,

Attached is a list of SRSs that fail. I have also merged Xiuguang's workaround into 1.8-esri (r22710); will revert if there's a different fix.

[warmerdam]

The problem seems to occur when citation parsing causes the linear units to be set *before* the projection parameters are set. They are set using SetNormProjParm?() under the assumption that the linear units will still be meters.

The suggested solution seems a bit odd. What is the significance of the unitsCode being KvUserDefined?? why is that different than if it is a well known linear unit?

In my investigation I am also highly suspicious of this code which slipped in during a previous esri merge a couple years ago:

        int unitCode = 0;
        GTIFKeyGet(hGTIF, ProjLinearUnitsGeoKey, &unitCode, 0, 1  );
        if(unitCode != KvUserDefined)
        {
            adfParm[5] /= psDefn->UOMLengthInMeters;
            adfParm[6] /= psDefn->UOMLengthInMeters;
        }

This code also has an odd relationship to unitCode not being KvUserDefined?.

The suggested change will certainly introduce unnecessary false easting and northing parameters in SRSes where it does not normally occur which is going to be bad for some picky WKT processors. I also don't feel that it is really the right solution. I feel it would be better to either defer application of the linear units from the citation till near the end - in the same area we normally set the linear units, or we could temporarily "unapply" the linear units before doing all the projection parameters and then reapply it after.

I'm not sure I'm up to tackling this tonight but I wanted to get my notes down. Also, if you guys would prefer to take a crack at doing it you have all day Friday.

[xz]

Frank, The following code

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

was in GDAL at least before 9/28/2005. It was before esri's merge. This code has being made a lot of problems. “If (unitCode != KvUserDefined?)” was to fix bug for some tif datasets that the unitCode is KvUserDefined?. In this case UOMLengthInMeters was a strange number. The false easting/northig became wrong by dividing an undefined number. “ If (unitCode != KvUserDefined?) “ does not introduce any bad result but handles some special tif datasets. I don’t think the current issue relates to it.

My fix is a workaround for the current bug. It sets several conditions to correct back the wrong false easting/northing but not affect others. It fixes at least 14 issues.

The current issue is mostly in the function of OGRSpatialReference::SetNormProjParm?. The root problem is dfValue /= dfToMeter in the function. It sould be removed at least for false easting/northing. OGRSpatialReference should independ of format. The parameters should set what it is. The convertion is not necessary here. If sorm formats have such requirement, it should be done in the format itself.

You may also need to check GTIFGetProjTRFInfo to see if the following makes trouble:

adfProjParms[i] = atof(pszValue) * dfInMeters;

I think it is a big problem too.

Xiuguang Zhou