Dump osm2rdf:length and osm2rdf:area as meter and m² (ad-freiburg#127)

Changes the units of osm2rdf:length to meters, and of osm2rdf:area to square meters.

Meter lengths are computed by summing up the great-circle distances between line anchor points using the haversine formular.

Square meter areas are computed by projecting polygon coordinates using an area-preserving projection (the Lambert cylindrical equal-area projection) and computing the area using Gauss's area formula. This simple projection assumes that the earth is a perfect sphere, but the approximation error should be negligible for typical polygons in OSM. For example, the approximation error of the area of Switzerland is below 0.3%.

Author: patrickbr

include/osm2rdf/osm/Constants.h

@@ -21,7 +21,8 @@
 
 namespace osm2rdf::osm::constants {
 
-static const int AREA_PRECISION = 12;
+static const int AREA_PRECISION = 4;
+static const int LENGTH_PRECISION = 2;
 static const int BASE10_BASE = 10;
 static const double BASE_SIMPLIFICATION_FACTOR = 0.001;
 

src/osm/Area.cpp

@@ -17,10 +17,11 @@
 // You should have received a copy of the GNU General Public License
 // along with osm2rdf.  If not, see <https://www.gnu.org/licenses/>.
 
+#include "osm2rdf/osm/Area.h"
+
 #include <iostream>
 #include <limits>
 
-#include "osm2rdf/osm/Area.h"
 #include "osm2rdf/osm/Constants.h"
 #include "osmium/osm/area.hpp"
 #include "util/geo/Geo.h"
@@ -33,7 +34,35 @@ osm2rdf::osm::Area::Area() {
 
 // ____________________________________________________________________________
 void osm2rdf::osm::Area::finalize() noexcept {
-  _geomArea = ::util::geo::area(_geom);
+  // poly in area preserving projection
+  ::util::geo::DMultiPolygon lambertPoly;
+  lambertPoly.resize(_geom.size());
+
+  double EARTH_RAD = 6371008.7714;  // mean radius
+
+  for (size_t i = 0; i < _geom.size(); i++) {
+    const auto& poly = _geom[i];
+    lambertPoly[i].getOuter().reserve(poly.getOuter().size());
+    for (const auto& p : poly.getOuter()) {
+      lambertPoly[i].getOuter().push_back(
+          ::util::geo::DPoint{EARTH_RAD * (p.getX() * util::geo::RAD),
+                              EARTH_RAD * (sin(p.getY() * util::geo::RAD))});
+    }
+
+    lambertPoly[i].getInners().reserve(poly.getInners().size());
+
+    for (const auto& inner : poly.getInners()) {
+      lambertPoly[i].getInners().push_back({});
+      lambertPoly[i].getInners().back().reserve(inner.size());
+      for (const auto& p : inner) {
+        lambertPoly[i].getInners().back().push_back(
+            ::util::geo::DPoint{EARTH_RAD * (p.getX() * util::geo::RAD),
+                                EARTH_RAD * (sin(p.getY() * util::geo::RAD))});
+      }
+    }
+  }
+
+  _geomArea = ::util::geo::area(lambertPoly);
 }
 
 // ____________________________________________________________________________
@@ -115,8 +144,7 @@ ::util::geo::DPolygon osm2rdf::osm::Area::convexHull() const noexcept {
 }
 
 // ____________________________________________________________________________
-::util::geo::DPolygon osm2rdf::osm::Area::orientedBoundingBox()
-    const noexcept {
+::util::geo::DPolygon osm2rdf::osm::Area::orientedBoundingBox() const noexcept {
   return ::util::geo::convexHull(::util::geo::getOrientedEnvelope(_geom));
 }
 
@@ -129,8 +157,7 @@ const ::util::geo::DPoint osm2rdf::osm::Area::centroid() const noexcept {
 bool osm2rdf::osm::Area::operator==(
     const osm2rdf::osm::Area& other) const noexcept {
   return _id == other._id && _objId == other._objId &&
-         _geomArea == other._geomArea &&
-         _envelope == other._envelope &&
+         _geomArea == other._geomArea && _envelope == other._envelope &&
          _geom == other._geom;
 }
 

src/osm/FactHandler.cpp

@@ -33,6 +33,7 @@
 #include "osm2rdf/ttl/Writer.h"
 
 using osm2rdf::osm::constants::AREA_PRECISION;
+using osm2rdf::osm::constants::LENGTH_PRECISION;
 using osm2rdf::osm::constants::BASE_SIMPLIFICATION_FACTOR;
 using osm2rdf::ttl::constants::CHANGESET_NAMESPACE;
 using osm2rdf::ttl::constants::DATASET_ID;
@@ -441,7 +442,8 @@ void osm2rdf::osm::FactHandler<W>::way(const osm2rdf::osm::Way& way) {
   }
 
   _writer->writeLiteralTripleUnsafe(
-      subj, IRI__OSM2RDF__LENGTH, std::to_string(::util::geo::len(way.geom())),
+      subj, IRI__OSM2RDF__LENGTH,
+      ::util::formatFloat(::util::geo::latLngLen(way.geom()), LENGTH_PRECISION),
       _iriXSDDouble);
 }
 

tests/issues/Issue24.cpp

@@ -75,7 +75,7 @@ TEST(Issue24, areaFromWayHasGeometryAsGeoSPARQL) {
       ".\nosmway:21 "
       "osm2rdfgeom:envelope \"POLYGON((48 7.5,48.1 7.5,48.1 7.6,48 7.6,48 "
       "7.5))\"^^geo:wktLiteral .\nosmway:21 "
-      "osm2rdf:area \"0.01\"^^xsd:double .\n",
+      "osm2rdf:area \"122568687.4654\"^^xsd:double .\n",
       buffer.str());
 
   // Cleanup
@@ -131,7 +131,7 @@ TEST(Issue24, areaFromRelationHasGeometryAsGeoSPARQL) {
       "osmrel:10 geo:hasGeometry osm2rdfgeom:osmrel_10 "
       ".\nosm2rdfgeom:osmrel_10 geo:asWKT \"POLYGON((48 7.5,48 "
       "7.6,48.1 7.6,48.1 7.5,48 7.5))\"^^geo:wktLiteral .\nosmrel:10 "
-      "osm2rdf:area \"0.01\"^^xsd:double .\n",
+      "osm2rdf:area \"122568687.4654\"^^xsd:double .\n",
       buffer.str());
 
   // Cleanup
@@ -348,7 +348,7 @@ TEST(Issue24, wayHasGeometryAsGeoSPARQL) {
       "osm2rdfgeom:envelope \"POLYGON((48 7.5,48.1 7.5,48.1 7.6,48 7.6,48 "
       "7.5))\"^^geo:wktLiteral .\nosmway:42 osm2rdfgeom:obb \"POLYGON((48.1 "
       "7.6,48.1 7.6,48 7.5,48 7.5,48.1 7.6))\"^^geo:wktLiteral "
-      ".\nosmway:42 osm2rdf:length \"0.141421\"^^xsd:double .\n",
+      ".\nosmway:42 osm2rdf:length \"15674.68\"^^xsd:double .\n",
       buffer.str());
 
   // Cleanup

tests/osm/FactHandler.cpp

@@ -98,7 +98,8 @@ TEST(OSM_FactHandler, areaFromWay) {
       "7.5))\"^^geo:wktLiteral .\nosmway:21 osm2rdfgeom:envelope \"POLYGON((48 "
       "7.5,48.1 7.5,48.1 7.6,48 7.6,48 7.5))\"^^geo:wktLiteral .\nosmway:21 "
       "osm2rdfgeom:obb \"POLYGON((48.1 7.5,48.1 7.6,48 7.6,48 7.5,48.1 "
-      "7.5))\"^^geo:wktLiteral .\nosmway:21 osm2rdf:area \"0.01\"^^xsd:double "
+      "7.5))\"^^geo:wktLiteral .\nosmway:21 osm2rdf:area "
+      "\"122568687.4654\"^^xsd:double "
       ".\n",
       buffer.str());
 
@@ -161,7 +162,8 @@ TEST(OSM_FactHandler, areaFromRelation) {
       "7.5))\"^^geo:wktLiteral .\nosmrel:10 osm2rdfgeom:envelope \"POLYGON((48 "
       "7.5,48.1 7.5,48.1 7.6,48 7.6,48 7.5))\"^^geo:wktLiteral .\nosmrel:10 "
       "osm2rdfgeom:obb \"POLYGON((48.1 7.5,48.1 7.6,48 7.6,48 7.5,48.1 "
-      "7.5))\"^^geo:wktLiteral .\nosmrel:10 osm2rdf:area \"0.01\"^^xsd:double "
+      "7.5))\"^^geo:wktLiteral .\nosmrel:10 osm2rdf:area "
+      "\"122568687.4654\"^^xsd:double "
       ".\n",
       buffer.str());
 
@@ -537,7 +539,7 @@ TEST(OSM_FactHandler, way) {
       "7.6,48 7.5))\"^^geo:wktLiteral .\n"
       "osmway:42 osm2rdfgeom:obb \"POLYGON((48.1 7.6,48.1 7.6,48 7.5,48 "
       "7.5,48.1 7.6))\"^^geo:wktLiteral .\n"
-      "osmway:42 osm2rdf:length \"0.141421\"^^xsd:double .\n",
+      "osmway:42 osm2rdf:length \"15674.68\"^^xsd:double .\n",
       buffer.str());
 
   // Cleanup
@@ -604,7 +606,7 @@ TEST(OSM_FactHandler, wayAddWayNodeOrder) {
       "osm2rdfgeom:envelope \"POLYGON((48 7.5,48.1 7.5,48.1 7.6,48 7.6,48 "
       "7.5))\"^^geo:wktLiteral .\nosmway:42 osm2rdfgeom:obb \"POLYGON((48.1 "
       "7.6,48.1 7.6,48 7.5,48 7.5,48.1 7.6))\"^^geo:wktLiteral .\nosmway:42 "
-      "osm2rdf:length \"0.141421\"^^xsd:double .\n",
+      "osm2rdf:length \"15674.68\"^^xsd:double .\n",
       buffer.str());
 
   // Cleanup
@@ -667,14 +669,15 @@ TEST(OSM_FactHandler, wayAddWayNodeSpatialMetadataShortWay) {
       ".\n_:0_1 osmway:member_id osmnode:2 .\n_:0_1 osmway:member_pos "
       "\"1\"^^xsd:integer .\n_:0_0 osmway:next_node osmnode:2 .\n_:0_0 "
       "osmway:next_node_distance \"15657.137001\"^^xsd:decimal .\nosmway:42 "
-      "geo:hasGeometry osm2rdfgeom:osmway_42 .\nosm2rdfgeom:osmway_42 geo:asWKT "
+      "geo:hasGeometry osm2rdfgeom:osmway_42 .\nosm2rdfgeom:osmway_42 "
+      "geo:asWKT "
       "\"LINESTRING(48 7.5,48.1 7.6)\"^^geo:wktLiteral .\nosmway:42 "
       "osm2rdfgeom:convex_hull \"POLYGON((48 7.5,48.1 7.6,48 "
       "7.5))\"^^geo:wktLiteral .\nosmway:42 osm2rdfgeom:envelope \"POLYGON((48 "
       "7.5,48.1 7.5,48.1 7.6,48 7.6,48 7.5))\"^^geo:wktLiteral .\nosmway:42 "
       "osm2rdfgeom:obb \"POLYGON((48.1 7.6,48.1 7.6,48 7.5,48 7.5,48.1 "
       "7.6))\"^^geo:wktLiteral .\nosmway:42 osm2rdf:length "
-      "\"0.141421\"^^xsd:double .\n",
+      "\"15674.68\"^^xsd:double .\n",
       buffer.str());
 
   // Cleanup
@@ -744,10 +747,11 @@ TEST(OSM_FactHandler, wayAddWayNodeSpatialMetadataLongerWay) {
       "\"3\"^^xsd:integer "
       ".\n_:0_2 osmway:next_node osmnode:3 .\n_:0_2 osmway:next_node_distance "
       "\"11024.108103\"^^xsd:decimal .\nosmway:42 geo:hasGeometry "
-      "osm2rdfgeom:osmway_42 .\nosm2rdfgeom:osmway_42 geo:asWKT \"LINESTRING(48 "
+      "osm2rdfgeom:osmway_42 .\nosm2rdfgeom:osmway_42 geo:asWKT "
+      "\"LINESTRING(48 "
       "7.5,48.1 "
       "7.6,48.1 7.5,48 7.5)\"^^geo:wktLiteral .\nosmway:42 osm2rdf:length "
-      "\"0.341421\"^^xsd:double .\n",
+      "\"37843.09\"^^xsd:double .\n",
       buffer.str());
 
   // Cleanup
@@ -817,7 +821,7 @@ TEST(OSM_FactHandler, wayAddWayMetaData) {
       "osmway:42 osmway:is_closed \"false\"^^xsd:boolean .\n"
       "osmway:42 osmway:nodeCount \"2\"^^xsd:integer .\n"
       "osmway:42 osmway:uniqueNodeCount \"2\"^^xsd:integer .\n"
-      "osmway:42 osm2rdf:length \"0.141421\"^^xsd:double .\n",
+      "osmway:42 osm2rdf:length \"15674.68\"^^xsd:double .\n",
       buffer.str());
 
   // Cleanup