Two similar Projections don't match (but should they?)

[MattReimer]

NB: This is out of an email from @JamesSLC I've attached a zipfile containing two rasters with projections that are the same in all but name:

2003_DEM_Proj1.tif; 2016_DEM_Proj2.tif ProjectionTestData.zip

Projection 1 is:

Projection 2 is:

So, on the surface of it ... the two projections differ only in terms of their name. This is an archetypal case where we might want to allow the user to override the projection clash. Presumably this would involve (by default) applying the projection of the first dataset added to the project to any later clashing datasets?

Update

So datasets that were showing up as having the same coordinate system in GCD 6, are not in GCD 7. Watch here: https://youtu.be/GgakxZD-gIA

For first 5 minutes. You'll see this problem occuring on just the Feshie dataset, and in later 15 minutes I test with CHAMP data and it works.

[MattReimer]

I did a bit of investigation here. When you look at the Pretty WKT representation of the two projections you get:

2016_DEM_Proj2.tif

LOCAL_CS["NZGD_2000_New_Zealand_Transverse_Mercator",
    GEOGCS["NZGD2000",
        DATUM["New_Zealand_Geodetic_Datum_2000",
            SPHEROID["GRS 1980",6378137,298.2572221010042,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6167"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4167"]],
    AUTHORITY["EPSG","2193"],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]

2003_DEM_Proj1.tif

PROJCS["NZGD_2000_Transverse_Mercator",
    GEOGCS["NZGD2000",
        DATUM["New_Zealand_Geodetic_Datum_2000",
            SPHEROID["GRS 1980",6378137,298.2572221010042,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6167"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4167"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",173],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",1600000],
    PARAMETER["false_northing",10000000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]

The issue here is that one of these is a projected PROJCS and the other is a local, ungeoreferenced coordinate system LOCAL_CS

From Gdal's perspective these are fundamentally different and can't really be compared.

We can use GDAL to compare the GEOG_CS component of the projection and in this case we find them to be the same but I'm not sure that solves anything but this one case.

If we want to go there we can but I'm not a projection expert and way out of my depth here. If we want to be able to support this comparison I think we have two basic options:

1. Design and implement our own IsSame() comparison operation in C#. This will require some really careful thought about what "sameness" means and a bunch of carefully designed test cases to approximate the kinds of differences we expect users to have.
2. Show the user the difference in an intelligent way and allow them to reset the projection to one representative of the current project. This is the easier option and will simply require some intuitive dialogs and figuring out how to execute the workflow.

See issue #162 for more about this.

Regardless of whether we choose 1 or 2 this is a bigger conversation and I don't want to risk destabilizing our current ability to get rasters in so my suggestion for the workshop in April is to leave things as is and correct projections manually before importing them.

[jb10016]

Hi Matt,

This is interesting, and rather odd. That the two projections appear to contain the identical parameters in the ESRI rendering of the coordinate system, but not the OGC representation is rather difficult to understand! As you say, creating a validation routine that checks 'sameness' is not going to be straightforward, as we would need to be extremely comprehensive to make this work, and that is likely to be a rabbit hole I, for one, don't fancy running down. Perhaps the best option here is to think through some guidance that pops up when we get a clash of projections, but which ultimately makes this the users problem and requires them to ensure consistency by maintaining a consistent projection file across the project. I suspect in the long run this might be the more sustainable option, as your bit of digging seems to suggest that this is deeper than 'name only'. Perhaps @joewheaton will have thoughts here ....

[philipbaileynar]

@jb10016 I suspect that the local versus projected representation of the same coordinate system originates from how the datasets were processed. Perhaps one came through CAD or some survey bespoke instrument software that was smart enough to pick the correct coordinate system, but not smart enough to position it on the Earth properly. This is one of the fundamental differences between GIS and CAD.

I believe that while the projections themselves are the same, the differences occur because one of them is not positioned on the Earth (it's merely a local representation of the coordinate system). See how the projected version above possesses meridian information for how to apply the coordinate system to the earth.

So... how were these two datasets generated?

[jb10016]

@philipbaileynar I'm not sure of the details regarding the generation of these two datasets ... the projections were attached to the files I inherited so I'm not in a position to assess this directly. However, what still seems odd to me nonetheless, is that the ERSI rendering in Matt's first post clearly possess all of the relevant projection parameters ... down to the central meridian/datums etc., but that this is not somehow correctly transferred into OGC metadata?

[MattReimer]

@jb10016 that does seem odd. Maybe ESRI guesses at them if they aren't there or maybe there's some aspect to this that I've missed completely. It's tricky with ESRI being closed source and with gdal's sparse documentation it might as well be.

I think the most important thing here is that every time this comes up we add the offending raster to a bucket of exceptions and eventually we'll know all the most common differences that can exist and try to account for as many of them as we can.

We can supplement that with instructions on how to manually add a new projection to andataset using arcmap for the cases we can't handle.