[in2]: add function to convert raster to nc file

avaframe/in2Trans/rasterUtils.py

@@ -6,6 +6,9 @@
 import logging
 import rasterio
 import numpy as np
+import re
+import subprocess
+import pathlib
 
 # create local logger
 log = logging.getLogger(__name__)
@@ -227,3 +230,90 @@ def getRasterFileTypeFromHeader(header):
         raise AssertionError("Unsupported driver for DEM %s" % header["driver"])
 
     return fileType
+
+
+def _extractTime(filename):
+    match = re.search(r"_t([\d.]+)\.", filename.name)
+    return float(match.group(1)) if match else float("inf")
+
+
+def convertRasterToNcFile(inputDir, infileExt, infileSuffix, outputDir, outFileName="", crs=None):
+    """convert .asc or .tif raster file to netCDF file,
+    if the rasters have timesteps in the name, they are sorted and numerated that they work with ParaView
+
+    Parameters
+    ----------
+    inputDir : pathlib object or str
+        path to folder containing the raster files
+    infileExt: str, list
+        file extension e.g. asc, tif
+    infileSuffix: str
+        file name part
+    outputDir : pathlib object or str
+        path to folder where netCDF files will be stored
+    outFileName: str
+        name of output netCDF file (if "", the same name of input raster file is used)
+    crs: str
+        overwrite crs
+    """
+    inputDir = pathlib.Path(inputDir)
+    outputDir = pathlib.Path(outputDir)
+
+    outputDir.mkdir(exist_ok=True)
+
+    if infileSuffix != "":
+        inputFiles = [file for file in inputDir.glob(f"*{infileSuffix}*.{infileExt}")]
+    else:
+        inputFiles = [file for file in inputDir.glob(f"*.{infileExt}")]
+
+    inputFiles.sort(
+        key=_extractTime,
+    )
+
+    numberConverted = 0
+    for i, filename in enumerate(inputFiles):
+        if outFileName == "":
+            cleanName = re.sub(r"_t[\d.]+", "", filename.stem)
+            outputFileName = f"{cleanName}_{i:04d}.nc"
+        else:
+            outputFileName = f"{outFileName}_{i:04d}.nc"
+
+        try:
+            if crs is None:
+                subprocess.run(
+                    [
+                        "gdal_translate",
+                        "-of",
+                        "netCDF",
+                        inputDir / filename.name,
+                        outputDir / outputFileName,
+                    ],
+                    check=True,
+                    capture_output=True,
+                    text=True,
+                )
+            else:
+                subprocess.run(
+                    [
+                        "gdal_translate",
+                        "-of",
+                        "netCDF",
+                        "-a_srs",
+                        crs,
+                        inputDir / filename.name,
+                        outputDir / outputFileName,
+                    ],
+                    check=True,
+                    capture_output=True,
+                    text=True,
+                )
+            log.info("%s is converted into %s and saved in %s" % (filename, outputFileName, outputDir))
+            numberConverted += 1
+        except subprocess.CalledProcessError as e:
+            message = "%s could not be converted: %s." % (filename, e.stderr)
+            log.info(message)
+
+    log.info(
+        "%s raster files (from %s available files) are converted and saved in %s"
+        % (numberConverted, len(inputFiles), outputDir)
+    )

avaframe/tests/test_rasterUtils.py

@@ -0,0 +1,51 @@
+"""Tests for module rasterUtils"""
+
+import pathlib
+import numpy as np
+import rasterio
+
+import avaframe.in2Trans.rasterUtils as rasterUtils
+
+
+def test_convertRasterToNcFile(tmp_path):
+    inputDir = pathlib.Path(tmp_path, "testDir")
+    inputDir.mkdir(exist_ok=True)
+
+    # save test rasters as tif
+    inputRasterNames = ["testRaster_t0.00", "testRaster_t0.01", "testRaster_t20.00"]
+
+    testRaster = np.zeros((10, 10))
+    cellsize = 10
+    nrows, ncols = testRaster.shape
+
+    header = {
+        "cellsize": cellsize,
+        "nrows": nrows,
+        "ncols": ncols,
+        "xllcenter": 0,
+        "yllcenter": 0,
+        "nodata_value": -9999,
+        "driver": "GTiff",
+        "crs": "EPSG:4326",
+    }
+    # convert lower-left center to upper-left corner
+    x_ul = header["xllcenter"] - cellsize / 2
+    y_ul = header["yllcenter"] + nrows * cellsize - cellsize / 2
+
+    transform = rasterio.transform.from_origin(x_ul, y_ul, cellsize, cellsize)
+    header["transform"] = transform
+    for name in inputRasterNames:
+        # write flipped raster, the read raster function does also flip the raster
+        rasterUtils.writeResultToRaster(header, testRaster, inputDir / name, useCompression=False, flip=True)
+    del testRaster
+
+    inFileExt = "tif"
+    inFileSuf = "test"
+    outputDir = inputDir / "ncFiles"
+
+    rasterUtils.convertRasterToNcFile(inputDir, inFileExt, inFileSuf, outputDir)
+
+    ncFileNames = ["testRaster_0000", "testRaster_0001", "testRaster_0002"]
+    for ncFileName in ncFileNames:
+        ncPath = outputDir / f"{ncFileName}.nc"
+        assert ncPath.is_file()