Re: [Geotools-gt2-users] Problems with AXIS order in CRS

[Jody G. <jod...@gm...>]

Yeah it is a pain, for a detailed breakdown on how the expectation changes
depending on protocol see:
https://docs.geoserver.org/latest/en/user/services/wfs/axis_order.html

The take home is if you see EPSG:4326 (and not the longer urn formats) is
to not trust it and look at the data.

The specification (the EPSG database) actually defines in northing/easting
order.

Shapefile files or WKB in databases can store geometry in the order given
imported, the spatial reference system is your responsibility to get
correct :)

Both the CRS decode and fromWKT methods should not be summarized as you
have done. Both are susceptible to system property hints you provide on the
command line, or how you have configured the spatial reference system
factories.  You can use both methods to return your choice of axis order <—
you must choose and there is no “always”.

CRS:84 was provided by OGC as an apology on this madness and can reliably
give you x/y order.

If at all possible use the longer URN formats that can be much much more
specific.  Even down to what version of the EPSG database was used for the
definition.
—
Jody


On Tue, Apr 28, 2020 at 2:53 PM David Adler <dwa...@gm...> wrote:

> Jody & Ian,
>
> Thank you for enlightening me regarding EPSG:4326 and the link to the
> documentation related to this.
>
> In 30+ years in GIS, I've worked with WGS84 coordinates in longitude -
> latitude order in shapefile, WKT, WKB and spatial data constructors. I had
> assume that EPSG:4326 was just the embodiment of that and didn't realize
> that the specification defined a different order. It took quite a while to
> diagnose the problem by debugging the code.
>
> It appears that if you use CRS.fromWKT(wkt), the default is EAST-NORTH as
> opposed to the NORTH-EAST returned by CRS.decode("EPSG:4326"). DB2 is
> always going to return the data in WKB or WKT in longitude-latitude order.
>
> Is there any guidance on whether it is better to construct the CRS from
> WKT or set the hint to "longitude first" and use EPSG:4326?
>
> On Tue, Apr 28, 2020 at 2:52 PM Jody Garnett <jod...@gm...>
> wrote:
>
>> How does DB2 documentation have it?
>>
>>
>>
>> https://www.ibm.com/support/knowledgecenter/SSEPEK_11.0.0/pdf/dasz_03_spatlbook.pdf
>>
>> Looks like the have alter_cs and alter_crs so you can configure. Can you
>> check the definition in your database David?
>>
>> Note “EPSG:4326” is ambiguous on axis order intended (different databases
>> and protocols have different assumptions). You can use a longer form If you
>> need wish to have the exact order from EPSG, or provide a hint of you wish
>> to force XY order (required for some protocols).
>>
>> See docs for details
>> https://docs.geotools.org/latest/userguide/library/referencing/order.html
>>
>>
>> On Tue, Apr 28, 2020 at 7:37 AM Ian Turton <ijt...@gm...> wrote:
>>
>>> It sounds like the data is not actually in EPSG:4326 LAT/LON order then,
>>> thus the renderer is getting confused by the fact that the CRS claims to be
>>> a NORTH_EAST projection when the data isn't. In general avoid using
>>> EPSG:4326 unless you are very sure that is what you have and you are sure
>>> about the axis order (usually, to the extent of checking the axis order
>>> while writing the coordinates into the geometries). You never know what the
>>> user may have set the hints to until you get there.
>>>
>>> Ian
>>>
>>> On Tue, 28 Apr 2020 at 14:35, David Adler <dwa...@gm...> wrote:
>>>
>>>> A DB2 sample program that displays a map based on spatial data in a DB2
>>>> table used to work but no longer works when I tried it recently. It just
>>>> fails to display anything.
>>>>
>>>> It seems to be related to the AXIS ordering in the CRS which is
>>>> different depending on whether one uses CRS.decode("EPSG:4326) or
>>>> CRS.parseWKT(wkt) where wkt is the string representation of EPSG 4326. I
>>>> fixed the problem by changing the DB2 code to use CRS.parseWKT instead of
>>>> CRS.decode with "org" and "orgid". Shouldn't it work both ways?
>>>>
>>>> The following is my analysis of the problem:
>>>>
>>>> Using CRS.decode("EPSG:4326")
>>>> results in
>>>> GEOGCS["WGS 84",
>>>>   DATUM["World Geodetic System 1984",
>>>>     SPHEROID["WGS 84", 6378137.0, 298.257223563,
>>>> AUTHORITY["EPSG","7030"]],
>>>>     AUTHORITY["EPSG","6326"]],
>>>>   PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]],
>>>>   UNIT["degree", 0.017453292519943295],
>>>>   AXIS["Geodetic latitude", NORTH],
>>>>   AXIS["Geodetic longitude", EAST],
>>>>   AUTHORITY["EPSG","4326"]]
>>>>
>>>> In StreamingRenderer.paint
>>>> CRS.getAxisOrder(mapCRS) returns CRS.AxisOrder.NORTH_EAST
>>>> which executes the following code which flips the AffineTransform in
>>>> WorldToScreen which completely messes up the subsequent transformation
>>>> from world to screen coordinates, causing them to be outside the
>>>> clip area.
>>>>         if (CRS.getAxisOrder(mapCRS) == CRS.AxisOrder.NORTH_EAST) {
>>>>             try {
>>>>                 // sanitize, having flipped axis causes slowdowns, the
>>>> rendering
>>>>                 // subsystem has to go from data to rendering to screen
>>>> flipping axis order
>>>>                 // twice when advanced projection handling is enabled
>>>>                 Integer code = CRS.lookupEpsgCode(mapCRS, false);
>>>>                 if (code != null) {
>>>>                     String srs = "EPSG:" + code;
>>>>                     CoordinateReferenceSystem earthNorthCRS =
>>>> CRS.decode(srs, true);
>>>>                     mapArea =
>>>>                             new ReferencedEnvelope(
>>>>                                     mapArea.getMinY(),
>>>>                                     mapArea.getMaxY(),
>>>>                                     mapArea.getMinX(),
>>>>                                     mapArea.getMaxX(),
>>>>                                     earthNorthCRS);
>>>>                 }
>>>>
>>>>                 // flip world to screen too
>>>>                 worldToScreen =
>>>>                         new AffineTransform(
>>>>                                 worldToScreen.getShearX(),
>>>>                                 worldToScreen.getScaleX(),
>>>>                                 worldToScreen.getScaleY(),
>>>>                                 worldToScreen.getShearY(),
>>>>                                 worldToScreen.getTranslateX(),
>>>>                                 worldToScreen.getTranslateY());
>>>>             } catch (Exception e) {
>>>>                 LOGGER.log(
>>>>                         Level.FINER,
>>>>                         "Failed to turn the requested bbox in
>>>> east/north order, map rendering "
>>>>                                 + "should work anyways, but pay a
>>>> performance price");
>>>>             }
>>>>
>>>> Creating the CRS from WKT, everything works fine.
>>>> CRS from CRS.decodeWKT(wkt)
>>>> GEOGCS["GCS_WGS_1984",
>>>>   DATUM["D_WGS_1984",
>>>>     SPHEROID["WGS_1984", 6378137.0, 298.257223563]],
>>>>   PRIMEM["Greenwich", 0.0],
>>>>   UNIT["degree", 0.017453292519943295],
>>>>   AXIS["Longitude", EAST],
>>>>   AXIS["Latitude", NORTH]]
>>>>
>>>> _______________________________________________
>>>> GeoTools-GT2-Users mailing list
>>>> Geo...@li...
>>>> https://lists.sourceforge.net/lists/listinfo/geotools-gt2-users
>>>>
>>>
>>>
>>> --
>>> Ian Turton
>>> _______________________________________________
>>> GeoTools-GT2-Users mailing list
>>> Geo...@li...
>>> https://lists.sourceforge.net/lists/listinfo/geotools-gt2-users
>>>
>> --
>> --
>> Jody Garnett
>>
> --
--
Jody Garnett