diff --git a/avaframe/in2Trans/rasterUtils.py b/avaframe/in2Trans/rasterUtils.py index a66d1b780..c3537c1cb 100644 --- a/avaframe/in2Trans/rasterUtils.py +++ b/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) + ) diff --git a/avaframe/tests/data/testShpConv/testLineGood.shx b/avaframe/tests/data/testShpConv/testLineGood.shx index ad7ed884f..27a1c2a1c 100644 Binary files a/avaframe/tests/data/testShpConv/testLineGood.shx and b/avaframe/tests/data/testShpConv/testLineGood.shx differ diff --git a/avaframe/tests/test_rasterUtils.py b/avaframe/tests/test_rasterUtils.py new file mode 100644 index 000000000..77991eb85 --- /dev/null +++ b/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()