[Liblas-commits] r1107 - in trunk: include/liblas/detail python/tests src/detail

liblas-commits at liblas.org liblas-commits at liblas.org
Tue Mar 17 12:58:25 EDT 2009


Author: hobu
Date: Tue Mar 17 12:58:22 2009
New Revision: 1107
URL: http://liblas.org/changeset/1107

Log:
Fix #108, LASWriter::SetSRS and LASReader::SetSRS don't work

Modified:
   trunk/include/liblas/detail/reader.hpp
   trunk/include/liblas/detail/writer.hpp
   trunk/python/tests/SRS.txt
   trunk/src/detail/reader.cpp
   trunk/src/detail/writer.cpp

Modified: trunk/include/liblas/detail/reader.hpp
==============================================================================
--- trunk/include/liblas/detail/reader.hpp	(original)
+++ trunk/include/liblas/detail/reader.hpp	Tue Mar 17 12:58:22 2009
@@ -89,9 +89,9 @@
     LASSRS m_in_srs;
     
 
+    OGRCoordinateTransformationH m_transform;
     OGRSpatialReferenceH m_in_ref;
     OGRSpatialReferenceH m_out_ref;
-    OGRCoordinateTransformationH m_transform;
 
 
 private:

Modified: trunk/include/liblas/detail/writer.hpp
==============================================================================
--- trunk/include/liblas/detail/writer.hpp	(original)
+++ trunk/include/liblas/detail/writer.hpp	Tue Mar 17 12:58:22 2009
@@ -81,6 +81,8 @@
     LASSRS m_in_srs;
     
     OGRCoordinateTransformationH m_transform;
+    OGRSpatialReferenceH m_in_ref;
+    OGRSpatialReferenceH m_out_ref;
     
 private:
 

Modified: trunk/python/tests/SRS.txt
==============================================================================
--- trunk/python/tests/SRS.txt	(original)
+++ trunk/python/tests/SRS.txt	Tue Mar 17 12:58:22 2009
@@ -28,16 +28,22 @@
   >>> test_srs()
   True
 
+#  ...         s2.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.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]"""
+#  ...         s2.wkt = """PROJCS["NAD83 / UTM zone 16N",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-87],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],AUTHORITY["EPSG","26916"],AXIS["Easting",EAST],AXIS["Northing",NORTH]]"""
+
   >>> def test_srs():
   ...     if not liblas.HAVE_GDAL: 
   ...         return True
   ...     if liblas.HAVE_GDAL:
   ...         s2 = srs.SRS()
   ...         s2.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.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]"""
+  ...         p = f.read(0)
+  ...         p.descale(f.header)
+  ...         if int(p.x) != 47069244: return False
   ...         f.set_srs(s2)
-  ...         p = f.get(0)
-  ...         print p.x, p.y
-  ...         return s.wkt == """PROJCS["WGS 84 / UTM zone 15N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-93],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","32615"]]"""
+  ...         p = f.read(0)
+  ...         p.descale(f.header)
+  ...         return int(round(p.x)) == -93  and int(round(p.y)) == 90
   ...     return False  
 
   >>> test_srs()

Modified: trunk/src/detail/reader.cpp
==============================================================================
--- trunk/src/detail/reader.cpp	(original)
+++ trunk/src/detail/reader.cpp	Tue Mar 17 12:58:22 2009
@@ -77,6 +77,13 @@
     if (m_transform) {
         OCTDestroyCoordinateTransformation(m_transform);
     }
+    if (m_in_ref) {
+        OSRDestroySpatialReference(m_in_ref);
+    }
+    if (m_out_ref) {
+        OSRDestroySpatialReference(m_out_ref);
+    }
+
 #endif
 }
 
@@ -92,6 +99,7 @@
     point.SetZ(record.z);
     
     if (m_transform) Project(point);
+
     point.SetIntensity(record.intensity);
     point.SetScanFlags(record.flags);
     point.SetClassification(record.classification);
@@ -165,16 +173,17 @@
     m_in_ref = OSRNewSpatialReference(0);
     m_out_ref = OSRNewSpatialReference(0);
 
-    const char* in_wkt = m_in_srs.GetWKT().c_str();
-    if (OSRImportFromWkt(m_in_ref, (char**) &in_wkt) != OGRERR_NONE) 
+    int result = OSRSetFromUserInput(m_in_ref, m_in_srs.GetWKT().c_str());
+    if (result != OGRERR_NONE) 
     {
-        throw std::runtime_error("Could not import input spatial reference for Reader::");
+        std::ostringstream msg; 
+        msg << "Could not import input spatial reference for Reader::" << CPLGetLastErrorMsg() << result;
+        std::string message(msg.str());
+        throw std::runtime_error(message);
     }
     
-    const char* out_wkt = m_out_srs.GetWKT().c_str();
-    printf("outwkt: %s", out_wkt);
-    int result = OSRImportFromWkt(m_out_ref, (char**) &out_wkt);
-    if (result!= OGRERR_NONE) 
+    result = OSRSetFromUserInput(m_out_ref, m_out_srs.GetWKT().c_str());
+    if (result != OGRERR_NONE) 
     {
         std::ostringstream msg; 
         msg << "Could not import output spatial reference for Reader::" << CPLGetLastErrorMsg() << result;
@@ -198,10 +207,14 @@
     
     ret = OCTTransform(m_transform, 1, &x, &y, &z);
     
-    if (ret != OGRERR_NONE) {
-        throw std::runtime_error("could not project point!");
+    if (!ret) {
+        std::ostringstream msg; 
+        msg << "Could not project point for Reader::" << CPLGetLastErrorMsg() << ret;
+        std::string message(msg.str());
+        throw std::runtime_error(message);
     }
     
+
     point.SetX(x);
     point.SetY(y);
     point.SetZ(z);

Modified: trunk/src/detail/writer.cpp
==============================================================================
--- trunk/src/detail/writer.cpp	(original)
+++ trunk/src/detail/writer.cpp	Tue Mar 17 12:58:22 2009
@@ -60,7 +60,7 @@
 
 namespace liblas { namespace detail {
 
-Writer::Writer(std::ostream& ofs) : m_ofs(ofs), m_transform(0)
+Writer::Writer(std::ostream& ofs) : m_ofs(ofs), m_transform(0), m_in_ref(0), m_out_ref(0)
 {
 }
 
@@ -69,7 +69,13 @@
 #ifdef HAVE_GDAL
     if (m_transform) {
         OCTDestroyCoordinateTransformation(m_transform);
-    }    
+    }
+    if (m_in_ref) {
+        OSRDestroySpatialReference(m_in_ref);
+    }
+    if (m_out_ref) {
+        OSRDestroySpatialReference(m_out_ref);
+    }
 #endif
 }
 
@@ -114,6 +120,65 @@
     }
 }
 
+
+void Writer::SetSRS(const LASSRS& srs)
+{
+    m_out_srs = srs;
+#ifdef HAVE_GDAL
+    m_in_ref = OSRNewSpatialReference(0);
+    m_out_ref = OSRNewSpatialReference(0);
+
+    int result = OSRSetFromUserInput(m_in_ref, m_in_srs.GetWKT().c_str());
+    if (result != OGRERR_NONE) 
+    {
+        std::ostringstream msg; 
+        msg << "Could not import input spatial reference for Reader::" << CPLGetLastErrorMsg() << result;
+        std::string message(msg.str());
+        throw std::runtime_error(message);
+    }
+    
+    result = OSRSetFromUserInput(m_out_ref, m_out_srs.GetWKT().c_str());
+    if (result != OGRERR_NONE) 
+    {
+        std::ostringstream msg; 
+        msg << "Could not import output spatial reference for Reader::" << CPLGetLastErrorMsg() << result;
+        std::string message(msg.str());
+        throw std::runtime_error(message);
+    }
+
+    m_transform = OCTNewCoordinateTransformation( m_in_ref, m_out_ref);
+    
+#endif
+}
+
+void Writer::Project(LASPoint& point)
+{
+#ifdef HAVE_GDAL
+    
+    // FIXME: This needs to ensure it is projecting scaled points
+    int ret = 0;
+    double x = point.GetX();
+    double y = point.GetY();
+    double z = point.GetZ();
+    
+    ret = OCTTransform(m_transform, 1, &x, &y, &z);
+    
+    if (!ret) {
+        std::ostringstream msg; 
+        msg << "Could not project point for Writer::" << CPLGetLastErrorMsg() << ret;
+        std::string message(msg.str());
+        throw std::runtime_error(message);
+    }
+    
+
+    point.SetX(x);
+    point.SetY(y);
+    point.SetZ(z);
+#else
+    UNREFERENCED_PARAMETER(point);
+#endif
+}
+
 void Writer::SetSRS(const LASSRS& srs)
 {
     m_out_srs = srs;


More information about the Liblas-commits mailing list