CRS comparison and parameter precision

[aaime]

Hi,
some months ago Milton started a thread about CRS
comparisons and the issues in comparing the
WKT definitions of two different systems:
http://n2.nabble.com/Axis-problems-again-with-CRS-lookupEpsgCode-td3214651.html#a3225238

Today I've stumbled in pretty much the same problem,
with the different that this time I have a mandate to
actually spend solid time looking into the issue deeper (instead
of trying to cook and answer in 5 minutes between higher priority
issues) and... I guess I have to apologize with Milton as I've found
some of my statements to be not true.

First off, the x/easting y/northing comparison should be working
and there is code to equate x to easting and y to northing...
it's just buggy:
http://jira.codehaus.org/browse/GEOT-2744

Then comes the accuracy comparison issue. I looked into a huge list
of sample .prj files and found that after adding enough aliases
to get away from the datum name issues (if you have more aliases
to add I'm all ears btw) and found the case of mismatch due
to different precision in parameters to be common.

I'm open to add a tolerance to the parameters comparisons... I'm
just not exactly sure how.

For example, the following comparisong do fail:

PROJCS["NAD27(CGQ77) / SCoPQ zone 2",
   GEOGCS["NAD27(CGQ77)",
     DATUM["North American Datum 1927 (CGQ77)",
       SPHEROID["Clarke 1866", 6378206.4, 294.9786982138982,
AUTHORITY["EPSG","7008"]],
       AUTHORITY["EPSG","6609"]],
     PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]],
     UNIT["degree", 0.017453292519943295],
     AXIS["Geodetic longitude", EAST],
     AXIS["Geodetic latitude", NORTH],
     AUTHORITY["EPSG","4609"]],
   PROJECTION["Transverse Mercator", AUTHORITY["EPSG","9807"]],
   PARAMETER["central_meridian", -55.5],
   PARAMETER["latitude_of_origin", 0.0],
   PARAMETER["scale_factor", 0.9999],
   PARAMETER["false_easting", 304800.0],
   PARAMETER["false_northing", 0.0],
   UNIT["m", 1.0],
   AXIS["Easting", EAST],
   AXIS["Northing", NORTH],
   AUTHORITY["EPSG","2008"]]

PROJCS["NAD_1927_CGQ77_MTM_2_SCoPQ",
   GEOGCS["GCS_NAD_1927_CGQ77",
     DATUM["D_NAD_1927_CGQ77",
       SPHEROID["Clarke_1866", 6378206.4, 294.9786982]],
     PRIMEM["Greenwich", 0.0],
     UNIT["degree", 0.017453292519943295],
     AXIS["Longitude", EAST],
     AXIS["Latitude", NORTH]],
   PROJECTION["Transverse_Mercator"],
   PARAMETER["central_meridian", -55.5],
   PARAMETER["latitude_of_origin", 0.0],
   PARAMETER["scale_factor", 0.9999],
   PARAMETER["false_easting", 304800.0],
   PARAMETER["false_northing", 0.0],
   UNIT["m", 1.0],
   AXIS["x", EAST],
   AXIS["y", NORTH]]

The comparison fails because of the inverse flattening value
differ in accuracy, 294.9786982 vs 294.9786982138982.
Using an absolute difference in accuracy would not work,
other parameters have different scales. What about a
percentage one?
Assuming then that the comparison between x and y
passes if  y - y * percentage <= x <= y + y * percentage.

We could have the tolerance be defined as a hint in
GeoTools and have a default that's scrict enough
(eventually plain 0) and allow everybody to relax it
based on his needs?

Cheers
Andrea




--
Andrea Aime
OpenGeo - http://opengeo.org
Expert service straight from the developers.

------------------------------------------------------------------------------
Come build with us! The BlackBerry® Developer Conference in SF, CA
is the only developer event you need to attend this year. Jumpstart your
developing skills, take BlackBerry mobile applications to market and stay
ahead of the curve. Join us from November 9-12, 2009. Register now!
http://p.sf.net/sfu/devconf
_______________________________________________
Geotools-devel mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/geotools-devel

[aaime]

Andrea Aime ha scritto:

> Hi,
> some months ago Milton started a thread about CRS
> comparisons and the issues in comparing the
> WKT definitions of two different systems:
> http://n2.nabble.com/Axis-problems-again-with-CRS-lookupEpsgCode-td3214651.html#a3225238
>
> Today I've stumbled in pretty much the same problem,
> with the different that this time I have a mandate to
> actually spend solid time looking into the issue deeper (instead
> of trying to cook and answer in 5 minutes between higher priority
> issues) and... I guess I have to apologize with Milton as I've found
> some of my statements to be not true.
>
> First off, the x/easting y/northing comparison should be working
> and there is code to equate x to easting and y to northing...
> it's just buggy:
> http://jira.codehaus.org/browse/GEOT-2744
>
> Then comes the accuracy comparison issue. I looked into a huge list
> of sample .prj files and found that after adding enough aliases
> to get away from the datum name issues (if you have more aliases
> to add I'm all ears btw) and found the case of mismatch due
> to different precision in parameters to be common.
>
> I'm open to add a tolerance to the parameters comparisons... I'm
> just not exactly sure how.
>
> For example, the following comparisong do fail:
>
> PROJCS["NAD27(CGQ77) / SCoPQ zone 2",
>    GEOGCS["NAD27(CGQ77)",
>      DATUM["North American Datum 1927 (CGQ77)",
>        SPHEROID["Clarke 1866", 6378206.4, 294.9786982138982,
> AUTHORITY["EPSG","7008"]],
>        AUTHORITY["EPSG","6609"]],
>      PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]],
>      UNIT["degree", 0.017453292519943295],
>      AXIS["Geodetic longitude", EAST],
>      AXIS["Geodetic latitude", NORTH],
>      AUTHORITY["EPSG","4609"]],
>    PROJECTION["Transverse Mercator", AUTHORITY["EPSG","9807"]],
>    PARAMETER["central_meridian", -55.5],
>    PARAMETER["latitude_of_origin", 0.0],
>    PARAMETER["scale_factor", 0.9999],
>    PARAMETER["false_easting", 304800.0],
>    PARAMETER["false_northing", 0.0],
>    UNIT["m", 1.0],
>    AXIS["Easting", EAST],
>    AXIS["Northing", NORTH],
>    AUTHORITY["EPSG","2008"]]
>
> PROJCS["NAD_1927_CGQ77_MTM_2_SCoPQ",
>    GEOGCS["GCS_NAD_1927_CGQ77",
>      DATUM["D_NAD_1927_CGQ77",
>        SPHEROID["Clarke_1866", 6378206.4, 294.9786982]],
>      PRIMEM["Greenwich", 0.0],
>      UNIT["degree", 0.017453292519943295],
>      AXIS["Longitude", EAST],
>      AXIS["Latitude", NORTH]],
>    PROJECTION["Transverse_Mercator"],
>    PARAMETER["central_meridian", -55.5],
>    PARAMETER["latitude_of_origin", 0.0],
>    PARAMETER["scale_factor", 0.9999],
>    PARAMETER["false_easting", 304800.0],
>    PARAMETER["false_northing", 0.0],
>    UNIT["m", 1.0],
>    AXIS["x", EAST],
>    AXIS["y", NORTH]]
>
> The comparison fails because of the inverse flattening value
> differ in accuracy, 294.9786982 vs 294.9786982138982.
> Using an absolute difference in accuracy would not work,
> other parameters have different scales. What about a
> percentage one?
> Assuming then that the comparison between x and y
> passes if  y - y * percentage <= x <= y + y * percentage.
>
> We could have the tolerance be defined as a hint in
> GeoTools and have a default that's scrict enough
> (eventually plain 0) and allow everybody to relax it
> based on his needs?

Btw, made a few tests, choosing percentage = 1^10-9
seems not to break the build (minus one test where the result
provided is different, but equivalent, just a different representation
of the same point on the ellipsoid) and still improve in a
significant way the number of matches against the sample WKT I do have.

Going towards a lower precision changes the kind of transformations
the referencing subsystem does. For example, if two datums are
close enough to match due to the relaxed precision the transformation
between them will be skipped, and so on: the referencing subsystem
users extensivedly identifiedObject.equalsIgnoreMetadata to decide
which steps of a general transformation actually use.

I would suggest to leave the default percentage to zero (complete
match) but to allow users to relax it a bit in case the ability
to match external wkt would be an important feature.

Cheers
Andrea

--
Andrea Aime
OpenGeo - http://opengeo.org
Expert service straight from the developers.

------------------------------------------------------------------------------
Come build with us! The BlackBerry® Developer Conference in SF, CA
is the only developer event you need to attend this year. Jumpstart your
developing skills, take BlackBerry mobile applications to market and stay
ahead of the curve. Join us from November 9-12, 2009. Register now!
http://p.sf.net/sfu/devconf
_______________________________________________
Geotools-devel mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/geotools-devel

[mbedward]

Hi Andrea,

> http://jira.codehaus.org/browse/GEOT-2744
>

I don't see any difference between the two snippets not matter how
hard I stare at them :)

> Using an absolute difference in accuracy would not work,
> other parameters have different scales.

Could you explain this... why wouldn't a per-parameter absolute
difference work with the example that you give ?

cheers
Michael

------------------------------------------------------------------------------
Come build with us! The BlackBerry® Developer Conference in SF, CA
is the only developer event you need to attend this year. Jumpstart your
developing skills, take BlackBerry mobile applications to market and stay
ahead of the curve. Join us from November 9-12, 2009. Register now!
http://p.sf.net/sfu/devconf
_______________________________________________
Geotools-devel mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/geotools-devel

[aaime]

Michael Bedward ha scritto:
> Hi Andrea,
>
>> http://jira.codehaus.org/browse/GEOT-2744
>>
>
> I don't see any difference between the two snippets not matter how
> hard I stare at them :)

SPHEROID["Clarke 1866", 6378206.4, 294.9786982138982,
AUTHORITY["EPSG","7008"]],

vs

DATUM["D_NAD_1927_CGQ77",
        SPHEROID["Clarke_1866", 6378206.4, 294.9786982]],

that is: 294.9786982138982 vs
          294.9786982

This is enough for the current referencing subsystem to state
the two projections are different (whilst by all intents and purposes
they are the same one, but using a different accuracy level for their
parameters).

>> Using an absolute difference in accuracy would not work,
>> other parameters have different scales.
>
> Could you explain this... why wouldn't a per-parameter absolute
> difference work with the example that you give ?

If we could have a per parameter accuracy it could be an absolute
difference, but then we'd end up with soo many absolute
differences.

Consider for example:
294.9786982138982
304800.0

Using an absolute difference of 10e-6 would make sense for the
first, but would not make much sense for the latter, which is
big. If we use a tolerance based on the size of the value we
get something that adapts to the order of magnitude of the
value itself. If we apply 10-6 as a percentage for the first value
we'd tolerate differences of 2.94 * 10e-4 and for the second
of 0.3 instead. Makes sense?

Cheers
Andrea

--
Andrea Aime
OpenGeo - http://opengeo.org
Expert service straight from the developers.

------------------------------------------------------------------------------
Come build with us! The BlackBerry® Developer Conference in SF, CA
is the only developer event you need to attend this year. Jumpstart your
developing skills, take BlackBerry mobile applications to market and stay
ahead of the curve. Join us from November 9-12, 2009. Register now!
http://p.sf.net/sfu/devconf
_______________________________________________
Geotools-devel mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/geotools-devel

[mbedward]

Hi Andrea

>> I don't see any difference between the two snippets not matter how
>> hard I stare at them :)
>
> SPHEROID["Clarke 1866", 6378206.4, 294.9786982138982,
> AUTHORITY["EPSG","7008"]],
>
> vs
>
> DATUM["D_NAD_1927_CGQ77",
>       SPHEROID["Clarke_1866", 6378206.4, 294.9786982]],
>

No - I mean't the code snippets on your Jira issue where you have:

The code:
xy = xy.trim();
        if (xy.length() == 0) {
            // try match
        }
should read:

xy = xy.trim();
        if (xy.length() == 0) {
            // try match
        }

Thanks for explaining your idea of adaptive tolerance !

Michael

------------------------------------------------------------------------------
Come build with us! The BlackBerry® Developer Conference in SF, CA
is the only developer event you need to attend this year. Jumpstart your
developing skills, take BlackBerry mobile applications to market and stay
ahead of the curve. Join us from November 9-12, 2009. Register now!
http://p.sf.net/sfu/devconf
_______________________________________________
Geotools-devel mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/geotools-devel

[aaime]

Michael Bedward ha scritto:

> Hi Andrea
>>> I don't see any difference between the two snippets not matter how
>>> hard I stare at them :)
>> SPHEROID["Clarke 1866", 6378206.4, 294.9786982138982,
>> AUTHORITY["EPSG","7008"]],
>>
>> vs
>>
>> DATUM["D_NAD_1927_CGQ77",
>>       SPHEROID["Clarke_1866", 6378206.4, 294.9786982]],
>>
>
> No - I mean't the code snippets on your Jira issue where you have:
>
> The code:
> xy = xy.trim();
>         if (xy.length() == 0) {
>             // try match
>         }
> should read:

Dang, you're right the second should be:

> xy = xy.trim();
>         if (xy.length() == 1) {
>             // try match
>         }

Fixed the jira too

Cheers
Andrea


--
Andrea Aime
OpenGeo - http://opengeo.org
Expert service straight from the developers.

------------------------------------------------------------------------------
Come build with us! The BlackBerry® Developer Conference in SF, CA
is the only developer event you need to attend this year. Jumpstart your
developing skills, take BlackBerry mobile applications to market and stay
ahead of the curve. Join us from November 9-12, 2009. Register now!
http://p.sf.net/sfu/devconf
_______________________________________________
Geotools-devel mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/geotools-devel