Description: Flip N/E systems to E/N and geodetic systems to Lon/Lat, while leaving Polar systems as-is, references #4949, 3.1 branch Author: Paul Ramsey Origin: https://trac.osgeo.org/postgis/changeset/8baf0b07b26df12d246c82bdae8ecd77371f3d24/git Bug: https://trac.osgeo.org/postgis/ticket/4949 Bug-Debian: https://bugs.debian.org/1035921 --- a/liblwgeom/lwgeom_transform.c +++ b/liblwgeom/lwgeom_transform.c @@ -261,43 +261,63 @@ proj_cs_get_simplecs(const PJ *pj_crs) return NULL; } +#define STR_EQUALS(A, B) strcmp((A), (B)) == 0 +#define STR_IEQUALS(A, B) (strcasecmp((A), (B)) == 0) +#define STR_ISTARTS(A, B) (strncasecmp((A), (B), strlen((B))) == 0) + static uint8_t proj_crs_is_swapped(const PJ *pj_crs) { - PJ *pj_cs; - uint8_t rv = LW_FALSE; + int axis_count; + PJ *pj_cs = proj_cs_get_simplecs(pj_crs); + if (!pj_cs) + lwerror("%s: proj_cs_get_simplecs returned NULL", __func__); + + axis_count = proj_cs_get_axis_count(NULL, pj_cs); + if (axis_count >= 2) + { + const char *out_name1, *out_abbrev1, *out_direction1; + const char *out_name2, *out_abbrev2, *out_direction2; + /* Read first axis */ + proj_cs_get_axis_info(NULL, + pj_cs, 0, + &out_name1, &out_abbrev1, &out_direction1, + NULL, NULL, NULL, NULL); + /* Read second axis */ + proj_cs_get_axis_info(NULL, + pj_cs, 1, + &out_name2, &out_abbrev2, &out_direction2, + NULL, NULL, NULL, NULL); + + proj_destroy(pj_cs); + + /* Directions agree, this is a northing/easting CRS, so reverse it */ + if(out_direction1 && STR_IEQUALS(out_direction1, "north") && + out_direction2 && STR_IEQUALS(out_direction2, "east") ) + { + return LW_TRUE; + } + + /* Oddball case? Both axes north / both axes south, swap */ + if(out_direction1 && out_direction2 && + ((STR_IEQUALS(out_direction1, "north") && STR_IEQUALS(out_direction2, "north")) || + (STR_IEQUALS(out_direction1, "south") && STR_IEQUALS(out_direction2, "south"))) && + out_name1 && STR_ISTARTS(out_name1, "northing") && + out_name2 && STR_ISTARTS(out_name2, "easting")) + { + return LW_TRUE; + } + + /* Any lat/lon system with Lat in first axis gets swapped */ + if (STR_ISTARTS(out_abbrev1, "Lat")) + return LW_TRUE; + + return LW_FALSE; + } - pj_cs = proj_cs_get_simplecs(pj_crs); - if (!pj_cs) { - lwerror("%s: proj_cs_get_simplecs returned NULL", __func__); - } - int axis_count = proj_cs_get_axis_count(NULL, pj_cs); - if (axis_count > 0) - { - const char *out_name, *out_abbrev, *out_direction; - double out_unit_conv_factor; - const char *out_unit_name, *out_unit_auth_name, *out_unit_code; - /* Read only first axis */ - proj_cs_get_axis_info(NULL, - pj_cs, - 0, - &out_name, - &out_abbrev, - &out_direction, - &out_unit_conv_factor, - &out_unit_name, - &out_unit_auth_name, - &out_unit_code); - - /* Only swap Lat/Lon systems */ - /* Use whatever ordering planar systems default to */ - if (strcasecmp(out_abbrev, "Lat") == 0) - rv = LW_TRUE; - else - rv = LW_FALSE; - } - proj_destroy(pj_cs); - return rv; + /* Failed the axis count test, leave quietly */ + proj_destroy(pj_cs); + return LW_FALSE; } LWPROJ * --- a/regress/core/tickets.sql +++ b/regress/core/tickets.sql @@ -1318,3 +1318,25 @@ SELECT '#4727', _ST_DistanceTree('SRID=4 SELECT '#4796', st_astext(st_snaptogrid(st_normalize(st_simplifypreservetopology('MULTISURFACE(((178632.044 397744.007,178631.118 397743.786,178646.399 397679.574,178693.864 397690.889,178698.958 397669.487,178700.206 397669.784,178758.532 397683.689,178748.351 397726.468,178752.199 397727.384,178748.782 397741.904,178744.897 397740.98,178738.157 397769.303,178632.044 397744.007)))'::geometry,1)),1)); SELECT '#4812', st_srid('SRID=999999;POINT(1 1)'::geometry); + +-- New Zealand forward -- SRID=2193;POINT(1766289 5927325) +SELECT '#4949', 'NZ forward', ST_AsEWKT(ST_SnapToGrid(ST_Transform( + 'SRID=4326;POINT(174.863597538742 -36.785298415230315)'::geometry, 2193),0.1)); +--- New Zealand inverse (opposite EPSG order) -- SRID=4326;POINT(174.863598 -36.785298) +SELECT '#4949', 'NZ inverse', ST_AsEWKT(ST_SnapToGrid(ST_Transform( + 'SRID=2193;POINT(1766289 5927325)'::geometry, 4326),0.000001)); +-- British Columbia forward (respect EPSG order) -- SRID=3005;POINT(1286630.44 561883.98) +SELECT '#4949', 'BC forward', ST_AsEWKT(ST_SnapToGrid(ST_Transform( + 'SRID=4269;POINT(-122 50)'::geometry, 3005),0.01)); +-- British Columbia inverse (respect EPSG order) -- SRID=4269;POINT(-122 50) +SELECT '#4949', 'BC inverse', ST_AsEWKT(ST_SnapToGrid(ST_Transform( + 'SRID=3005;POINT(1286630.44 561883.98)'::geometry, 4269),0.000001)); +-- North Pole LAEA Europe inverse -- SRID=4326;POINT(19.4921659 69.7902258) +SELECT '#4949', 'North Pole LAEA inverse', ST_AsEWKT(ST_SnapToGrid(ST_Transform( + 'SRID=3575;POINT(370182 -2213980)'::geometry,4326),0.0000001)); +-- Polar Stereographic forward -- SRID=3413;POINT(2218082.1 -1409150) +SELECT '#4949', 'Arctic Stereographic forward', ST_AsEWKT(ST_SnapToGrid(ST_Transform( + 'SRID=4326;POINT(12.572160 66.081084)'::geometry,3413),0.1)); +-- Antarctic Polar Stereographic -- SRID=3031;POINT(-2399498.7 3213318.5) +SELECT '#4949', 'Antarctic Stereographic forward', ST_AsEWKT(ST_SnapToGrid(ST_Transform( + 'SRID=4326;POINT(-36.75 -54.25)'::geometry, 3031),0.1)); --- a/regress/core/tickets_expected +++ b/regress/core/tickets_expected @@ -440,3 +440,10 @@ ERROR: LWGEOM_addpoint: Invalid offset #4727|0 #4796|POLYGON((178632 397744,178738 397769,178745 397741,178749 397742,178752 397727,178748 397726,178759 397684,178699 397669,178694 397691,178646 397680,178632 397744)) #4812|999999 +#4949|NZ forward|SRID=2193;POINT(1766289 5927325) +#4949|NZ inverse|SRID=4326;POINT(174.863598 -36.785298) +#4949|BC forward|SRID=3005;POINT(1286630.44 561883.98) +#4949|BC inverse|SRID=4269;POINT(-122 50) +#4949|North Pole LAEA inverse|SRID=4326;POINT(19.4921659 69.7902258) +#4949|Arctic Stereographic forward|SRID=3413;POINT(2218082.1 -1409150) +#4949|Antarctic Stereographic forward|SRID=3031;POINT(-2399498.7 3213318.5)