whiff on implementing axis ordering for filters.reprojection

filters/ReprojectionFilter.cpp

@@ -66,12 +66,19 @@ void ReprojectionFilter::addArgs(ProgramArgs& args)
 {
     args.add("out_srs", "Output spatial reference", m_outSRS).setPositional();
     args.add("in_srs", "Input spatial reference", m_inSRS);
+    args.add("in_axis_ordering", "Axis ordering override for in_srs", m_inAxisOrdering, {} );
+    args.add("out_axis_ordering", "Axis ordering override for out_srs", m_outAxisOrdering, {} );
 }
 
 
 void ReprojectionFilter::initialize()
 {
     m_inferInputSRS = m_inSRS.empty();
+    if (m_inAxisOrdering.size())
+        m_inSRS.setAxisOrdering(m_inAxisOrdering);
+
+    if (m_outAxisOrdering.size())
+        m_outSRS.setAxisOrdering(m_outAxisOrdering);
     setSpatialReference(m_outSRS);
 }
 

filters/ReprojectionFilter.hpp

@@ -68,6 +68,8 @@ class PDAL_DLL ReprojectionFilter : public Filter, public Streamable
     SpatialReference m_outSRS;
     bool m_inferInputSRS;
     std::unique_ptr<SrsTransform> m_transform;
+    std::vector<int> m_inAxisOrdering;
+    std::vector<int> m_outAxisOrdering;
 };
 
 } // namespace pdal

pdal/SpatialReference.cpp

@@ -81,6 +81,12 @@ SpatialReference::SpatialReference(const std::string& s)
     set(s);
 }
 
+SpatialReference::SpatialReference(const std::string& s,
+                                   const std::vector<int>& ordering)
+{
+    m_axisOrdering = ordering;
+    set(s);
+}
 
 bool SpatialReference::empty() const
 {
@@ -139,6 +145,18 @@ void SpatialReference::parse(const std::string& s, std::string::size_type& pos)
 }
 
 
+/// Set axis mapping strategy
+void SpatialReference::setAxisOrdering(const std::vector<int>& v)
+{
+    m_axisOrdering = v;
+}
+
+const std::vector<int>& SpatialReference::getAxisOrdering() const
+{
+
+    return m_axisOrdering;
+}
+
 void SpatialReference::set(std::string v)
 {
     m_horizontalWkt.clear();
@@ -170,6 +188,12 @@ void SpatialReference::set(std::string v)
         throw pdal_error(oss.str());
     }
 
+    if (m_axisOrdering.size())
+    {
+        // Automatically implies SetAxisMappingStrategy(OAMS_CUSTOM)
+        srs.SetDataAxisToSRSAxisMapping(m_axisOrdering);
+    }
+
     char *poWKT = 0;
     srs.exportToWkt(&poWKT);
     m_wkt = poWKT;
@@ -439,7 +463,7 @@ int SpatialReference::getUTMZone() const
 
 int SpatialReference::computeUTMZone(const BOX3D& cbox) const
 {
-    SrsTransform transform(*this, SpatialReference("EPSG:4326"));
+    SrsTransform transform(*this, SpatialReference("EPSG:4326", {2, 1}));
 
     // We made the argument constant so copy so that we can modify.
     BOX3D box(cbox);

pdal/SpatialReference.hpp

@@ -72,6 +72,15 @@ class PDAL_DLL SpatialReference
     */
     SpatialReference(const std::string& wkt);
 
+    /**
+      Construct a spatial reference from well-known text.
+
+      \param wkt  Well-known text from which to construct SRS.
+      \param ordering  Axis override ordering.
+    */
+    SpatialReference(const std::string& wkt, const std::vector<int>& ordering);
+
+
 
     /**
       Determine if this spatial reference is the same as another.
@@ -152,6 +161,11 @@ class PDAL_DLL SpatialReference
     void dump() const;
     MetadataNode toMetadata() const;
 
+
+    /// Set axis mapping strategy
+    void setAxisOrdering(const std::vector<int>& v);
+    const std::vector<int>& getAxisOrdering() const;
+
     bool isGeographic() const;
     bool isGeocentric() const;
     bool isProjected() const;
@@ -166,6 +180,8 @@ class PDAL_DLL SpatialReference
 private:
     std::string m_wkt;
     mutable std::string m_horizontalWkt;
+    std::vector<int> m_axisOrdering;
+
     friend PDAL_DLL std::ostream& operator<<(std::ostream& ostr,
         const SpatialReference& srs);
     friend PDAL_DLL std::istream& operator>>(std::istream& istr,

pdal/private/SrsTransform.cpp

@@ -51,8 +51,10 @@ SrsTransform::SrsTransform(const SpatialReference& src,
 // discussion for more info.
 //
 #if GDAL_VERSION_MAJOR >= 3
-    srcRef.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
-    dstRef.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
+    if (src.getAxisOrdering().size())
+        srcRef.SetDataAxisToSRSAxisMapping(src.getAxisOrdering());
+    if (dst.getAxisOrdering().size())
+        dstRef.SetDataAxisToSRSAxisMapping(dst.getAxisOrdering());
 #endif
     m_transform.reset(OGRCreateCoordinateTransformation(&srcRef, &dstRef));
 }

test/unit/SpatialReferenceTest.cpp

@@ -393,6 +393,8 @@ TEST(SpatialReferenceTest, merge)
 
     Options o2;
     o2.add("filename", Support::datapath("las/test_epsg_4326.las"));
+    o2.add("in_axis_ordering", "2");
+    o2.add("in_axis_ordering", "1");
     LasReader r2;
     r2.setOptions(o2);
 
@@ -404,6 +406,8 @@ TEST(SpatialReferenceTest, merge)
 
     Options o4;
     o4.add("out_srs", "EPSG:4326");
+    o4.add("out_axis_ordering", "2");
+    o4.add("out_axis_ordering", "1");
     ReprojectionFilter repro;
     repro.setOptions(o4);
     repro.setInput(r1);
@@ -439,7 +443,7 @@ TEST(SpatialReferenceTest, test_bounds)
     std::string wgs84_wkt = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]";
 
 
-    SpatialReference wgs84(wgs84_wkt);
+    SpatialReference wgs84(wgs84_wkt, {2, 1});
     BOX2D box17(289814.15, 4320978.61, 289818.50, 4320980.59);
     pdal::Polygon p(box17);
     p.setSpatialReference(utm17);
@@ -513,6 +517,36 @@ TEST(SpatialReferenceTest, set_srs)
     EXPECT_NE(m.value().find("AUTHORITY[\"EPSG\",\"2029\"]]"),
         std::string::npos);
 }
+
+
+// Test setting EPSG:4326 from User string
+TEST(SpatialReferenceTest, axis_ordering)
+{
+    SpatialReference wgs84("EPSG:4326", {2, 1});
+
+    std::vector<int> axis = wgs84.getAxisOrdering();
+
+    EXPECT_EQ(axis[0], 2);
+    EXPECT_EQ(axis[1], 1);
+
+
+    SpatialReference utm17("EPSG:26917");
+    BOX2D box_utm17(289814.15, 4320978.61, 289818.50, 4320980.59);
+//     BOX2D box_dd(-83.427597f, 39.0126f, -83.427551f, 39.01261f);
+    BOX2D box_dd(39.0126f, -83.427597f, 39.0126f, -83.427551f);
+    pdal::Polygon p(box_utm17);
+    p.setSpatialReference(utm17);
+    p.transform(wgs84);
+
+    BOX3D b2 = p.bounds();
+    EXPECT_FLOAT_EQ(static_cast<float>(b2.minx), -83.427597f);
+    EXPECT_FLOAT_EQ(static_cast<float>(b2.miny), 39.0126f);
+    EXPECT_FLOAT_EQ(static_cast<float>(b2.maxx), -83.427551f);
+    EXPECT_FLOAT_EQ(static_cast<float>(b2.maxy), 39.01261f);
+
+}
+
+
 #endif // PDAL_HAVE_LIBXML2
 
 } // namespace pdal