From f81f116581cc864a6f2b5718641c1518f3d7783e Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Mon, 4 Mar 2024 23:00:56 +0100
Subject: [PATCH] OGR VRT: add support for coordinate precision

---
 autotest/ogr/ogr_vrt.py             | 28 +++++++++++++++++++++++++++
 doc/source/drivers/vector/vrt.rst   |  8 ++++++++
 ogr/ogrsf_frmts/vrt/data/ogrvrt.xsd |  3 +++
 ogr/ogrsf_frmts/vrt/ogr_vrt.h       |  2 ++
 ogr/ogrsf_frmts/vrt/ogrvrtlayer.cpp | 30 +++++++++++++++++++++++++++++
 5 files changed, 71 insertions(+)

diff --git a/autotest/ogr/ogr_vrt.py b/autotest/ogr/ogr_vrt.py
index f6b1d8659c86..ddf5fcdce121 100755
--- a/autotest/ogr/ogr_vrt.py
+++ b/autotest/ogr/ogr_vrt.py
@@ -3330,3 +3330,31 @@ def test_ogr_vrt_field_names_same_case():
     f = lyr.GetNextFeature()
     assert f["id"] == "foo"
     assert f["id_from_uc"] == "bar"
+
+
+###############################################################################
+# Test geometry coordinate precision support
+
+
+def test_ogr_vrt_geom_coordinate_preicsion():
+
+    ds = ogr.Open(
+        """<OGRVRTDataSource>
+  <OGRVRTLayer name="poly">
+    <SrcDataSource>data/poly.shp</SrcDataSource>
+    <GeometryField>
+        <GeometryType>wkbPolygon</GeometryType>
+        <XYResolution>1e-5</XYResolution>
+        <ZResolution>1e-3</ZResolution>
+        <MResolution>1e-2</MResolution>
+    </GeometryField>
+  </OGRVRTLayer>
+</OGRVRTDataSource>
+"""
+    )
+    lyr = ds.GetLayer(0)
+    geom_fld = lyr.GetLayerDefn().GetGeomFieldDefn(0)
+    prec = geom_fld.GetCoordinatePrecision()
+    assert prec.GetXYResolution() == 1e-5
+    assert prec.GetZResolution() == 1e-3
+    assert prec.GetMResolution() == 1e-2
diff --git a/doc/source/drivers/vector/vrt.rst b/doc/source/drivers/vector/vrt.rst
index b1de74ac1050..82c32fa8ae21 100644
--- a/doc/source/drivers/vector/vrt.rst
+++ b/doc/source/drivers/vector/vrt.rst
@@ -210,6 +210,14 @@ layer name, and may have the following subelements:
      **SrcRegion**
   *  **ExtentXMin**, **ExtentYMin**, **ExtentXMax** and **ExtentXMax**
      (optional) : same syntax as OGRVRTLayer-level elements of same name
+  *  **XYResolution** (optional, GDAL >= 3.9):
+     Resolution for the coordinate precision of the X and Y coordinates.
+     Expressed in the units of the X and Y axis of the SRS
+  *  **ZResolution** (optional, GDAL >= 3.9):
+     Resolution for the coordinate precision of the Z coordinates.
+     Expressed in the units of the Z axis of the SRS
+  *  **MResolution** (optional, GDAL >= 3.9):
+     Resolution for the coordinate precision of the M coordinates.
 
   If no **GeometryField** element is specified, all the geometry fields of
   the source layer will be exposed by the VRT layer. In order not to
diff --git a/ogr/ogrsf_frmts/vrt/data/ogrvrt.xsd b/ogr/ogrsf_frmts/vrt/data/ogrvrt.xsd
index 6adc7911cedb..d7458319a691 100644
--- a/ogr/ogrsf_frmts/vrt/data/ogrvrt.xsd
+++ b/ogr/ogrsf_frmts/vrt/data/ogrvrt.xsd
@@ -393,6 +393,9 @@
                 <xs:element name="SrcRegion" type="SrcRegionType"/>
                 <xs:element name="SRS" type="nonEmptyStringType"/>
                 <xs:group ref="ExtentType"/>
+                <xs:element name="XYResolution" type="xs:double" minOccurs="0" maxOccurs="1"/>
+                <xs:element name="ZResolution" type="xs:double" minOccurs="0" maxOccurs="1"/>
+                <xs:element name="MResolution" type="xs:double" minOccurs="0" maxOccurs="1"/>
             </xs:choice>
         </xs:sequence>
         <xs:attributeGroup ref="GeometryFieldTypeAttrGroupWithoutSrc"/>
diff --git a/ogr/ogrsf_frmts/vrt/ogr_vrt.h b/ogr/ogrsf_frmts/vrt/ogr_vrt.h
index a4587c955b0f..8ed870f7b016 100644
--- a/ogr/ogrsf_frmts/vrt/ogr_vrt.h
+++ b/ogr/ogrsf_frmts/vrt/ogr_vrt.h
@@ -83,6 +83,8 @@ class OGRVRTGeomFieldProps
 
     OGREnvelope sStaticEnvelope;
 
+    OGRGeomCoordinatePrecision sCoordinatePrecision{};
+
     OGRVRTGeomFieldProps();
     ~OGRVRTGeomFieldProps();
 };
diff --git a/ogr/ogrsf_frmts/vrt/ogrvrtlayer.cpp b/ogr/ogrsf_frmts/vrt/ogrvrtlayer.cpp
index 12fb32820457..4decb5cf644c 100644
--- a/ogr/ogrsf_frmts/vrt/ogrvrtlayer.cpp
+++ b/ogr/ogrsf_frmts/vrt/ogrvrtlayer.cpp
@@ -513,6 +513,34 @@ bool OGRVRTLayer::ParseGeometryField(CPLXMLNode *psNode,
     poProps->bNullable =
         CPLTestBool(CPLGetXMLValue(psNode, "nullable", "TRUE"));
 
+    if (GetSrcLayerDefn()->GetGeomFieldCount() == 1)
+    {
+        poProps->sCoordinatePrecision =
+            GetSrcLayerDefn()->GetGeomFieldDefn(0)->GetCoordinatePrecision();
+    }
+    else if (poProps->eGeometryStyle == VGS_Direct && poProps->iGeomField >= 0)
+    {
+        poProps->sCoordinatePrecision =
+            GetSrcLayerDefn()
+                ->GetGeomFieldDefn(poProps->iGeomField)
+                ->GetCoordinatePrecision();
+    }
+    if (const char *pszXYResolution =
+            CPLGetXMLValue(psNode, "XYResolution", nullptr))
+    {
+        poProps->sCoordinatePrecision.dfXYResolution = CPLAtof(pszXYResolution);
+    }
+    if (const char *pszZResolution =
+            CPLGetXMLValue(psNode, "ZResolution", nullptr))
+    {
+        poProps->sCoordinatePrecision.dfZResolution = CPLAtof(pszZResolution);
+    }
+    if (const char *pszMResolution =
+            CPLGetXMLValue(psNode, "MResolution", nullptr))
+    {
+        poProps->sCoordinatePrecision.dfMResolution = CPLAtof(pszMResolution);
+    }
+
     return true;
 }
 
@@ -811,6 +839,8 @@ bool OGRVRTLayer::FullInitialize()
                                     apoGeomFieldProps[i]->eGeomType);
         oFieldDefn.SetSpatialRef(apoGeomFieldProps[i]->poSRS);
         oFieldDefn.SetNullable(apoGeomFieldProps[i]->bNullable);
+        oFieldDefn.SetCoordinatePrecision(
+            apoGeomFieldProps[i]->sCoordinatePrecision);
         poFeatureDefn->AddGeomFieldDefn(&oFieldDefn);
     }