refactor: epsg codes enum

Author: amfage

packages/topo-imagery-common/src/topo_imagery_common/epsg.py

@@ -0,0 +1,12 @@
+from enum import Enum
+
+
+class EpsgNumber(int, Enum):
+    NZTM_2000 = 2193
+    """New Zealand Transverse Mercator 2000"""
+    WGS_1984 = 4326
+    """World Geodetic System 1984"""
+    CITM_2000 = 3793
+    """Chatham Islands Transverse Mercator 2000"""
+    NZVD_2016 = 7839
+    """New Zealand Vertical Datum 2016"""

packages/topo-imagery-gdal/src/topo_imagery_gdal/gdal/gdal_commands.py

@@ -3,8 +3,8 @@
 from decimal import Decimal
 
 from linz_logger import get_log
+from topo_imagery_common.epsg import EpsgNumber
 from topo_imagery_gdal.gdal.gdal_bands import get_gdal_band_offset
-from topo_imagery_gdal.gdal.gdal_helper import EpsgNumber
 from topo_imagery_gdal.gdal.gdal_presets import (
     BASE_COG,
     COMPRESS_LZW,
@@ -101,7 +101,7 @@ def get_cutline_command(cutline: str | None) -> list[str]:
 def get_build_vrt_command(
     files: list[str],
     output: str = "output.vrt",
-    epsg: int = 2193,
+    epsg: int = EpsgNumber.NZTM_2000,
     add_alpha: bool = False,
     resolution: list[Decimal] | None = None,
 ) -> list[str]:

packages/topo-imagery-gdal/src/topo_imagery_gdal/gdal/gdal_helper.py

@@ -1,7 +1,6 @@
 import json
 import os
 import subprocess
-from enum import Enum
 from shutil import rmtree
 from tempfile import mkdtemp
 from typing import cast
@@ -10,6 +9,7 @@
 from topo_imagery_common.aws.aws_helper import is_s3
 from topo_imagery_common.files.files_helper import get_file_name_from_path
 from topo_imagery_common.files.fs import copy
+from topo_imagery_common.epsg import EpsgNumber
 from topo_imagery_common.log.time_helper import time_in_ms
 from topo_imagery_gdal.gdal.gdalinfo import GdalInfo
 
@@ -18,13 +18,6 @@ class GDALExecutionException(Exception):
     pass
 
 
-class EpsgNumber(int, Enum):
-    NZTM_2000 = 2193
-    """New Zealand Transverse Mercator 2000"""
-    WGS_1984 = 4326
-    """World Geodetic System 1984"""
-
-
 def get_vfs_path(path: str) -> str:
     """Make the path as a GDAL Virtual File Systems path.
 

packages/topo-imagery-gdal/src/topo_imagery_gdal/standardising.py

@@ -11,6 +11,7 @@
 from tifffile import TiffFile
 from topo_imagery_common.aws.aws_helper import is_s3
 from topo_imagery_common.cli.cli_helper import TileFiles
+from topo_imagery_common.epsg import EpsgNumber
 from topo_imagery_common.files.files_helper import ContentType, is_tiff
 from topo_imagery_common.files.fs import exists, read, write, write_all, write_sidecars
 from topo_imagery_common.log.time_helper import time_in_ms
@@ -211,7 +212,7 @@ def standardising(
 def create_vrt(
     source_tiffs: list[str],
     target_path: str,
-    epsg: int = 2193,
+    epsg: int = EpsgNumber.NZTM_2000,
     add_alpha: bool = False,
     resolution: list[Decimal] | None = None,
 ) -> str:

packages/topo-imagery-gdal/src/topo_imagery_gdal/tile/tile_index.py

@@ -1,6 +1,7 @@
 import re
 from typing import NamedTuple
 
+from topo_imagery_common.epsg import EpsgNumber
 from topo_imagery_gdal.tile.util import charcodeat
 
 SHEET_WIDTH = 24_000
@@ -16,9 +17,9 @@
 CHAR_A = charcodeat("A", 0)
 CHAR_S = charcodeat("S", 0)
 
-MAINLAND_EPSG = 2193
+MAINLAND_EPSG = EpsgNumber.NZTM_2000
 """ EPSG code of the mainland NZTM50 mapsheet grid (NZGD2000 / NZTM2000) """
-CHATHAM_EPSG = 3793
+CHATHAM_EPSG = EpsgNumber.CITM_2000
 """ EPSG code of the Chatham Islands mapsheet grid (NZGD2000 / CITM2000) """
 
 
@@ -63,41 +64,21 @@ class Bounds(NamedTuple):
     "CI06": Point(x=3_506_000, y=5_104_000),  # Pitt Island (Rangiauria)
 }
 
-
-def get_chatham_mapsheet_offset(sheet_code: str) -> Point:
-    """Look up the origin point for a Chatham Islands mapsheet code.
-
-    Args:
-        sheet_code: Chatham Islands topo 50 map sheet code eg "CI06"
-
-    Returns:
-        Point: The top left point of the mapsheet, in EPSG:3793
-
-    Example:
-        >>> get_chatham_mapsheet_offset("CI06")
-        Point(x=3506000, y=5104000)
-    """
-    origin = CHATHAM_SHEET_ORIGINS.get(sheet_code[:4])
-    if origin is None:
-        raise ValueError(f"Unknown Chatham Islands map sheet: {sheet_code}; known sheets: {sorted(CHATHAM_SHEET_ORIGINS)}")
-    return origin
-
-
 def get_bounds_from_name(tile_name: str, target_epsg: int = MAINLAND_EPSG) -> Bounds:
     """Get the origin coordinates and size of the tile from its name.
 
     Args:
         tile_name: the tile name as `sheetCode_gridSize_tileId`
-        target_epsg: EPSG code of the mapsheet grid the tile is named against; only
+        target_epsg: EPSG code of the mapsheet grid for the tile.
             `MAINLAND_EPSG` (2193) and `CHATHAM_EPSG` (3793) are supported
 
     Returns:
         a `Bounds` object
     """
-    if target_epsg == CHATHAM_EPSG:
-        get_offset = get_chatham_mapsheet_offset
-    elif target_epsg == MAINLAND_EPSG:
+    if target_epsg == MAINLAND_EPSG:
         get_offset = get_mapsheet_offset
+    elif target_epsg == CHATHAM_EPSG:
+        get_offset = get_chatham_mapsheet_offset
     else:
         raise ValueError(
             f"Unsupported target EPSG:{target_epsg} for mapsheet lookup; supported: {MAINLAND_EPSG}, {CHATHAM_EPSG}"
@@ -169,6 +150,26 @@ def get_mapsheet_offset(sheet_code: str) -> Point:
     return Point(x=SHEET_WIDTH * x + SHEET_ORIGIN_LEFT, y=SHEET_ORIGIN_TOP - SHEET_HEIGHT * y)
 
 
+def get_chatham_mapsheet_offset(sheet_code: str) -> Point:
+    """Look up the origin point for a Chatham Islands mapsheet code.
+
+    Args:
+        sheet_code: Chatham Islands topo 50 map sheet code eg "CI06"
+
+    Returns:
+        Point: The top left point of the mapsheet, in EPSG:3793
+
+    Example:
+        >>> get_chatham_mapsheet_offset("CI06")
+        Point(x=3506000, y=5104000)
+    """
+    origin = CHATHAM_SHEET_ORIGINS.get(sheet_code[:4])
+    if origin is None:
+        raise ValueError(f"Unknown Chatham Islands map sheet: {sheet_code}. Known sheets: {sorted(CHATHAM_SHEET_ORIGINS)}")
+    return origin
+
+
+
 def get_tile_offset(grid_size: int, x: int, y: int) -> Bounds:
     """Get the tile offset from its coordinate and the grid size
 

packages/topo-imagery-gdal/test/gdal/gdal_commands_test.py

@@ -2,13 +2,13 @@
 
 import pytest
 from pytest_subtests import SubTests
+from topo_imagery_common.epsg import EpsgNumber
 from topo_imagery_gdal.gdal.gdal_commands import (
     get_buffer_distance,
     get_cutline_command,
     get_footprint_command,
     get_gdal_command,
 )
-from topo_imagery_gdal.gdal.gdal_helper import EpsgNumber
 from topo_imagery_gdal.gdal.gdal_presets import CompressionPreset, HillshadePreset
 
 

packages/topo-imagery-gdal/test/tile/tile_index_data.py

@@ -1,6 +1,6 @@
 from typing import Any
 
-# Values transcribed from https://data.linz.govt.nz/layer/50089-nz-chatham-island-linz-map-sheets-topo-150k/
+# Values source: https://data.linz.govt.nz/layer/50089-nz-chatham-island-linz-map-sheets-topo-150k/
 CHATHAM_SHEET_DATA: list[dict[Any, Any]] = [
     {"origin": {"x": 3458000, "y": 5176000}, "code": "CI01"},
     {"origin": {"x": 3482000, "y": 5176000}, "code": "CI02"},

scripts/generate_hillshade.py

@@ -9,6 +9,7 @@
 from topo_imagery_common.cli.cli_helper import InputParameterError, TileFiles, load_input_files
 from topo_imagery_common.cli.common_args import CommonArgumentParser
 from topo_imagery_common.datetimes import RFC_3339_DATETIME_FORMAT
+from topo_imagery_common.epsg import EpsgNumber
 from topo_imagery_common.files.files_helper import SUFFIX_JSON, ContentType, is_tiff
 from topo_imagery_common.files.fs import exists, read, write, write_all
 from topo_imagery_common.log.time_helper import time_in_ms
@@ -112,7 +113,7 @@ def create_hillshade(
 
         # COGify the hillshade output, using ZSTD compression
         run_gdal(
-            get_gdal_command(CompressionPreset.DEM_ZSTD.value, 2193),
+            get_gdal_command(CompressionPreset.DEM_ZSTD.value, EpsgNumber.NZTM_2000),
             input_file=hillshade_working_path,
             output_file=hillshade_cog_working_path,
         )

scripts/pdal/pdal_commands.py

@@ -5,6 +5,7 @@
 
 from linz_logger import get_log
 from topo_imagery_common.aws.aws_helper import is_s3
+from topo_imagery_common.epsg import EpsgNumber
 from topo_imagery_common.files.fs import copy
 from topo_imagery_common.log.time_helper import time_in_ms
 
@@ -53,7 +54,7 @@ def get_pdal_command(command: str, options: list[str]) -> list[str]:
 pdal_translate_add_proj_command = get_pdal_command(
     "translate",
     [
-        "--readers.las.spatialreference=EPSG:2193+7839",
+        f"--readers.las.spatialreference=EPSG:{EpsgNumber.NZTM_2000.value}+{EpsgNumber.NZVD_2016.value}",
         "--writers.las.filesource_id=0",
         "--writers.las.forward=all",
     ],