From a8df62fe09ce8c39753c420988158c99b6433939 Mon Sep 17 00:00:00 2001 From: David Shean Date: Thu, 25 Jan 2024 21:59:46 -0800 Subject: [PATCH] gdal_edit: add -a_coord_epoch option with test, update doc --- autotest/pyscripts/test_gdal_edit.py | 23 ++++++++++++ doc/source/programs/gdal_edit.rst | 15 +++++++- doc/source/user/coordinate_epoch.rst | 19 +++++----- .../gdal-utils/osgeo_utils/gdal_edit.py | 35 +++++++++++++++++++ 4 files changed, 82 insertions(+), 10 deletions(-) mode change 100644 => 100755 swig/python/gdal-utils/osgeo_utils/gdal_edit.py diff --git a/autotest/pyscripts/test_gdal_edit.py b/autotest/pyscripts/test_gdal_edit.py index f56a28f9cf79..6d10bf9aaafc 100755 --- a/autotest/pyscripts/test_gdal_edit.py +++ b/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) diff --git a/doc/source/programs/gdal_edit.rst b/doc/source/programs/gdal_edit.rst index 6bc0bbe07deb..053787c85576 100644 --- a/doc/source/programs/gdal_edit.rst +++ b/doc/source/programs/gdal_edit.rst @@ -15,9 +15,10 @@ Synopsis .. code-block:: - gdal_edit [--help] [--help-general] [-ro] [-a_srs ] + gdal_edit [--help] [--help-general] [-ro] [-a_srs ] [-a_ullr ] [-a_ulurll ] [-tr ] [-unsetgt] [-unsetrpc] [-a_nodata ] [-unsetnodata] + [-a_coord_epoch ] [-unsetepoch] [-unsetstats] [-stats] [-approx_stats] [-setstats ] [-scale ] [-offset ] [-units ] @@ -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 + + 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 Assign/override the georeferenced bounds of the dataset. diff --git a/doc/source/user/coordinate_epoch.rst b/doc/source/user/coordinate_epoch.rst index 9ce7e5ac242c..4e9e852377ae 100644 --- a/doc/source/user/coordinate_epoch.rst +++ b/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. diff --git a/swig/python/gdal-utils/osgeo_utils/gdal_edit.py b/swig/python/gdal-utils/osgeo_utils/gdal_edit.py old mode 100644 new mode 100755 index 341c89e8dc3a..a8864e6a13c4 --- a/swig/python/gdal-utils/osgeo_utils/gdal_edit.py +++ b/swig/python/gdal-utils/osgeo_utils/gdal_edit.py @@ -52,6 +52,7 @@ def Usage(isError): " [-colorinterp_ {red|green|blue|alpha|gray|undefined]]...", file=f, ) + print(" [-a_coord_epoch ] [-unsetepoch]", file=f) print(" [-unsetstats] [-stats] [-approx_stats]", file=f) print(" [-setstats ]", 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,