Merge pull request OSGeo#9144 from dshean/gdal_edit_epoch

gdal_edit: add -a_coord_epoch option and update doc

autotest/pyscripts/test_gdal_edit.py

@@ -444,3 +444,26 @@ def test_gdal_edit_py_unsetrpc(script_path):
         ds = None
     finally:
         gdal.GetDriverByName("GTiff").Delete(filename)
+
+
+###############################################################################
+# Test -a_coord_epoch
+
+
+def test_gdal_edit_py_epoch(script_path):
+
+    filename = "tmp/test_gdal_edit_py.tif"
+    shutil.copy(test_py_scripts.get_data_path("gcore") + "byte.tif", filename)
+
+    try:
+        test_py_scripts.run_py_script(
+            script_path, "gdal_edit", filename + " -a_coord_epoch 2021.3"
+        )
+
+        ds = gdal.Open(filename)
+        srs = ds.GetSpatialRef()
+        assert srs.GetCoordinateEpoch() == 2021.3
+        ds = None
+
+    finally:
+        gdal.GetDriverByName("GTiff").Delete(filename)

doc/source/programs/gdal_edit.rst

@@ -15,9 +15,10 @@ Synopsis
 
 .. code-block::
 
-    gdal_edit [--help] [--help-general] [-ro] [-a_srs <srs_def>]
+    gdal_edit [--help] [--help-general] [-ro] [-a_srs <srs_def>] 
             [-a_ullr <ulx> <uly> <lrx> <lry>] [-a_ulurll <ulx> <uly> <urx> <ury> <llx> <lly>]
             [-tr <xres> <yres>] [-unsetgt] [-unsetrpc] [-a_nodata <value>] [-unsetnodata]
+            [-a_coord_epoch <epoch>] [-unsetepoch]
             [-unsetstats] [-stats] [-approx_stats]
             [-setstats <min> <max> <mean> <stddev>]
             [-scale <value>] [-offset <value>] [-units <value>]
@@ -58,6 +59,18 @@ It works only with raster formats that support update access to existing dataset
     coordinate system will be removed (for TIFF/GeoTIFF, might not be well
     supported besides that).
 
+.. option:: -a_coord_epoch <epoch>
+
+    Assign/override the coordinate epoch of the dataset, as a decimal year (e.g., 2021.3).
+
+    .. versionadded:: 3.9
+
+.. option:: -unsetepoch
+
+    Remove the coordinate epoch of the dataset.
+
+    .. versionadded:: 3.9
+
 .. option:: -a_ullr  <ulx> <uly> <lrx> <lry>
 
     Assign/override the georeferenced bounds of the dataset.

doc/source/user/coordinate_epoch.rst

@@ -29,21 +29,20 @@ a CRS is a dynamic one.
 The :cpp:func:`OGRSpatialReference::SetCoordinateEpoch` and
 :cpp:func:`OGRSpatialReference::GetCoordinateEpoch` methods can be used to
 set/retrieve a coordinate epoch associated with a CRS. The coordinate epoch is
-expressed as a decimal year (e.g. 2021.3).
+expressed as a decimal year (e.g. 2021.3 for April 21, 2021).
 
 Formally, the coordinate epoch of an observation belongs to the
-observation.  However, almost all formats do not allow for storing
+observation. However, almost all formats do not allow for storing
 per-observation epoch, and typical usage is a set of observations with
-the same epoch.  Therefore we store the epoch as property of the CRS,
-and the meaning of this is as if that epoch value were part of every
-observation.  This choice eases processing, storage and format
-complexity for most usage.  For now, this means that a dataset where
-points have different epochs cannot be handled.
+the same epoch. Therefore we store the epoch as property of the CRS,
+and assume that it is valid for every observation. This choice eases processing,
+storage and format complexity for most usage. For now, this means that a dataset 
+containing observations or points with different epochs cannot be handled.
 
 For vector formats, per-geometry coordinate epoch could also make sense, but as
 most formats only support a per-layer CRS, we also for now limit support of
-coordinate epoch at the layer level. The coordinate transformation mechanics
-itself can support per-vertex coordinate epoch.
+coordinate epoch at the layer level. The underlying coordinate transformation mechanics
+can support per-vertex coordinate epoch.
 
 Support in raster and vector formats
 ------------------------------------
@@ -162,6 +161,8 @@ Support in utilities
 :program:`gdalinfo` and :program:`ogrinfo` report the coordinate epoch, when
 attached to a dataset/layer SRS.
 
+:program:`gdal_edit.py` has a ``-a_coord_epoch`` option to define the epoch of a dataset in place. 
+
 :program:`gdal_translate` and :program:`ogr2ogr` have a ``-a_coord_epoch`` option to be used
 together with ``-a_srs``, and otherwise preserve the coordinate epoch in the output SRS
 from the source SRS when no SRS related options are specified.

swig/python/gdal-utils/osgeo_utils/gdal_edit.py

@@ -52,6 +52,7 @@ def Usage(isError):
         "                 [-colorinterp_<X> {red|green|blue|alpha|gray|undefined]]...",
         file=f,
     )
+    print("                 [-a_coord_epoch <epoch>] [-unsetepoch]", file=f)
     print("                 [-unsetstats] [-stats] [-approx_stats]", file=f)
     print("                 [-setstats <min> <max> <mean> <stddev>]", file=f)
     print(
@@ -107,6 +108,8 @@ def gdal_edit(argv):
     xres = None
     yres = None
     unsetgt = False
+    epoch = None
+    unsetepoch = False
     unsetstats = False
     stats = False
     setstats = False
@@ -131,6 +134,9 @@ def gdal_edit(argv):
         elif argv[i] == "-a_srs" and i < len(argv) - 1:
             srs = argv[i + 1]
             i = i + 1
+        elif argv[i] == "-a_coord_epoch" and i < len(argv) - 1:
+            epoch = float(argv[i + 1])
+            i = i + 1
         elif argv[i] == "-a_ullr" and i < len(argv) - 4:
             ulx = float(argv[i + 1])
             i = i + 1
@@ -199,6 +205,8 @@ def gdal_edit(argv):
             unsetgt = True
         elif argv[i] == "-unsetrpc":
             unsetrpc = True
+        elif argv[i] == "-unsetepoch":
+            unsetepoch = True
         elif argv[i] == "-unsetstats":
             unsetstats = True
         elif argv[i] == "-approx_stats":
@@ -279,9 +287,11 @@ def gdal_edit(argv):
 
     if (
         srs is None
+        and epoch is None
         and lry is None
         and yres is None
         and not unsetgt
+        and not unsetepoch
         and not unsetstats
         and not stats
         and not setstats
@@ -330,6 +340,11 @@ def gdal_edit(argv):
         print("", file=sys.stderr)
         return Usage(isError=True)
 
+    if unsetepoch and epoch:
+        print("-unsetepoch and -a_coord_epoch options are exclusive.", file=sys.stderr)
+        print("", file=sys.stderr)
+        return Usage(isError=True)
+
     if open_options is not None:
         if ro:
             ds = gdal.OpenEx(datasetname, gdal.OF_RASTER, open_options=open_options)
@@ -379,6 +394,26 @@ def gdal_edit(argv):
         if not gcp_list:
             ds.SetProjection(wkt)
 
+    if epoch is not None:
+        sr = ds.GetSpatialRef()
+        if sr is None:
+            print(
+                "Dataset SRS is undefined, cannot set epoch. See the -a_srs option.",
+                file=sys.stderr,
+            )
+            return -1
+        sr.SetCoordinateEpoch(epoch)
+        ds.SetSpatialRef(sr)
+
+    if unsetepoch:
+        sr = ds.GetSpatialRef()
+        if sr is None:
+            print("Dataset SRS is undefined, with no epoch specified.", file=sys.stderr)
+            return -1
+        # Set to 0.0, which is what GetCoordinateEpoch() returns for SRS with no epoch defined
+        sr.SetCoordinateEpoch(0.0)
+        ds.SetSpatialRef(sr)
+
     if lry is not None:
         gt = [
             ulx,