diff --git a/avaframe/com1DFA/com1DFA.py b/avaframe/com1DFA/com1DFA.py index be07a7b10..55b9740a6 100644 --- a/avaframe/com1DFA/com1DFA.py +++ b/avaframe/com1DFA/com1DFA.py @@ -1184,7 +1184,14 @@ def initializeSimulation(cfg, outDir, demOri, inputSimLines, logName): if cfg["EXPORTS"].getboolean("exportRasters"): outDir = pathlib.Path(cfgGen["avalancheDir"], "Outputs", "internalRasters") fU.makeADir(outDir) - IOf.writeResultToRaster(dem["originalHeader"], relRaster, (outDir / "releaseRaster"), flip=True) + useCompression = cfg["EXPORTS"].getboolean("useCompression") + IOf.writeResultToRaster( + dem["originalHeader"], + relRaster, + (outDir / "releaseRaster"), + flip=True, + useCompression=useCompression, + ) log.info( "Release area raster derived from %s saved to %s" % (releaseLine["initializedFrom"], str(outDir / "releaseRaster")) @@ -1268,11 +1275,13 @@ def initializeSimulation(cfg, outDir, demOri, inputSimLines, logName): # export secondary release raster used for computations (after cutting potential overlap with release) if cfg["EXPORTS"].getboolean("exportRasters"): outDir = pathlib.Path(cfg["GENERAL"]["avalancheDir"], "Outputs", "internalRasters") + useCompression = cfg["EXPORTS"].getboolean("useCompression") IOf.writeResultToRaster( dem["originalHeader"], secRelRaster, (outDir / ("secondaryReleaseRaster_%d" % secIndex)), flip=True, + useCompression=useCompression, ) log.info( "SecondaryRelease area raster derived from %s saved to %s" @@ -1284,11 +1293,13 @@ def initializeSimulation(cfg, outDir, demOri, inputSimLines, logName): # export entrainment raster used for computations (after cutting potential overlap with release or secondary release) if cfg["EXPORTS"].getboolean("exportRasters"): outDir = pathlib.Path(cfg["GENERAL"]["avalancheDir"], "Outputs", "internalRasters") + useCompression = cfg["EXPORTS"].getboolean("useCompression") IOf.writeResultToRaster( dem["originalHeader"], entrMassRaster / cfg["GENERAL"].getfloat("rhoEnt"), (outDir / "entrainmentRaster"), flip=True, + useCompression=useCompression, ) log.info( "Entrainment area raster derived from %s saved to %s" @@ -2017,8 +2028,7 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si # make sure to save all desiered resuts for first and last time step for # the report resTypesReport = fU.splitIniValueToArraySteps(cfg["REPORT"]["plotFields"]) - # always add particles to first and last time step - resTypesLast = list(set(resTypes + resTypesReport + ["particles"])) + resTypesLast = list(set(resTypes + resTypesReport)) # derive friction type # turn friction model into integer frictModelsList = [ @@ -2059,7 +2069,7 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si # setup a result fields info data frame to save max values of fields and avalanche front resultsDF = setupresultsDF(resTypesLast, cfg["VISUALISATION"].getboolean("createRangeTimeDiagram")) - # TODO: add here different time stepping options + # Add different time stepping options here log.debug("Use standard time stepping") # Initialize time and counters nSave = 1 @@ -2074,6 +2084,12 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si # export initial time step if cfg["EXPORTS"].getboolean("exportData"): exportFields(cfg, t, fields, dem, outDir, cuSimName, TSave="initial") + + if "particles" in resTypes: + outDirData = outDir / "particles" + fU.makeADir(outDirData) + savePartToPickle(particles, outDirData, cuSimName) + # export particles properties for visulation if cfg["VISUALISATION"].getboolean("writePartToCSV"): particleTools.savePartToCsv( @@ -2085,10 +2101,6 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si countParticleCsv = countParticleCsv + 1 # export particles dictionaries of saving time steps - # (if particles is not in resType, only first and last time step are saved) - outDirData = outDir / "particles" - fU.makeADir(outDirData) - savePartToPickle(particles, outDirData, cuSimName) zPartArray0 = copy.deepcopy(particles["z"]) @@ -2187,7 +2199,8 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si exportFields(cfg, t, fields, dem, outDir, cuSimName, TSave="intermediate") # export particles dictionaries of saving time steps - savePartToPickle(particles, outDirData, cuSimName) + if "particles" in resTypes: + savePartToPickle(particles, outDirData, cuSimName) # export particles properties for visulation if cfg["VISUALISATION"].getboolean("writePartToCSV"): @@ -2315,7 +2328,8 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si exportFields(cfg, t, fields, dem, outDir, cuSimName, TSave="final") # export particles dictionaries of saving time steps - savePartToPickle(particles, outDirData, cuSimName) + if "particles" in resTypes: + savePartToPickle(particles, outDirData, cuSimName) else: # fetch contourline info contourDictXY = outCom1DFA.fetchContCoors( @@ -2371,7 +2385,7 @@ def setupresultsDF(resTypes, cfgRangeTime): resultsDF: dataframe data frame with on line for iniital time step and max and mean values of fields """ - + # TODO catch empty resTypes resDict = {"timeStep": [0.0]} for resT in resTypes: if resT != "particles" and resT != "FTDet": @@ -2566,7 +2580,7 @@ def computeEulerTimeStep( cfg["ResistanceModel"].lower() == "default" and cfg.getboolean("detrainment") ): # update resistance area fields using thresholds - fields = com1DFATools.updateResCoeffFields(fields, cfg, float(particles["t"]), dem) + fields = com1DFATools.updateResCoeffFields(fields, cfg) if debugPlot: outCom1DFA.plotResFields(fields, cfg, particles["tPlot"], dem, particles["mTot"]) @@ -2603,7 +2617,6 @@ def computeEulerTimeStep( particles = DFAfunC.updatePositionC(cfg, particles, dem, force, fields, typeStop=0) tCPUPos = time.time() - startTime tCPU["timePos"] = tCPU["timePos"] + tCPUPos - # Split particles if cfg.getint("splitOption") == 0: # split particles with too much mass @@ -2968,12 +2981,21 @@ def exportFields( # convert from J/cell to kJ/m² # (by dividing the peak kinetic energy per cell by the real area of the cell) resField = resField * 0.001 / dem["areaRaster"] + dataName = cuSimName + "_" + resType + "_" + "t%.2f" % (timeStep) # create directory outDirPeak = outDir / "peakFiles" / "timeSteps" fU.makeADir(outDirPeak) outFile = outDirPeak / dataName - IOf.writeResultToRaster(dem["originalHeader"], resField, outFile, flip=True) + useCompression = cfg["EXPORTS"].getboolean("useCompression") + IOf.writeResultToRaster( + dem["originalHeader"], resField, outFile, flip=True, useCompression=useCompression + ) + log.debug( + "Results parameter: %s has been exported to Outputs/peakFiles for time step: %.2f " + % (resType, timeStep) + ) + if TSave == "final": log.debug( "Results parameter: %s exported to Outputs/peakFiles for time step: %.2f - FINAL time step " @@ -2984,11 +3006,9 @@ def exportFields( outDirPeakAll = outDir / "peakFiles" fU.makeADir(outDirPeakAll) outFile = outDirPeakAll / dataName - IOf.writeResultToRaster(dem["originalHeader"], resField, outFile, flip=True) - else: - log.debug( - "Results parameter: %s has been exported to Outputs/peakFiles for time step: %.2f " - % (resType, timeStep) + useCompression = cfg["EXPORTS"].getboolean("useCompression") + IOf.writeResultToRaster( + dem["originalHeader"], resField, outFile, flip=True, useCompression=useCompression ) diff --git a/avaframe/com1DFA/com1DFACfg.ini b/avaframe/com1DFA/com1DFACfg.ini index 707b34aa6..e00963570 100644 --- a/avaframe/com1DFA/com1DFACfg.ini +++ b/avaframe/com1DFA/com1DFACfg.ini @@ -545,4 +545,6 @@ unitpfv = ms-1 exportData = True # export release and optional entrainment raster files derived from shp files saved to Outputs/com1DFA/internalRasters exportRasters = False +# use LZW compression when writing TIFF raster files +useCompression = True diff --git a/avaframe/com1DFA/com1DFATools.py b/avaframe/com1DFA/com1DFATools.py index 10f813c73..5fde3a4a7 100644 --- a/avaframe/com1DFA/com1DFATools.py +++ b/avaframe/com1DFA/com1DFATools.py @@ -399,7 +399,7 @@ def checkCfgInfoType(cfgInfo): return typeCfgInfo -def updateResCoeffFields(fields, cfg, t, dem): +def updateResCoeffFields(fields, cfg): """update fields of cRes and detK, coefficients of resistance parameter and detrainment parameter according to the thresholds of FV and FT if FV OR FT below min thresholds -> apply detrainment in that area @@ -412,10 +412,6 @@ def updateResCoeffFields(fields, cfg, t, dem): dictionary with cResRasterOrig, cResRaster and detRasterOrig, detRaster fields cfg: configparser object configuration of com1DFA, thresholds - t: float - time step - dem: dict - dictionary with info on DEM Returns ------------ @@ -453,10 +449,6 @@ def updateResCoeffFields(fields, cfg, t, dem): "Resistance area removed %d cells because FV or FT exceeded %.2f ms-1, %.2f m" % (lTh, vMax, thMax) ) - outFileRes = pathlib.Path( - cfg["avalancheDir"], "Outputs", "com1DFA", "reports", ("resArea_t%.2f" % t) - ) - IOf.writeResultToRaster(dem["originalHeader"], fields["cResRasterOrig"], outFileRes, flip=True) # update fields dictionary fields["cResRaster"] = cResRaster diff --git a/avaframe/com1DFA/deriveParameterSet.py b/avaframe/com1DFA/deriveParameterSet.py index c5c8e8469..ae633f395 100644 --- a/avaframe/com1DFA/deriveParameterSet.py +++ b/avaframe/com1DFA/deriveParameterSet.py @@ -1017,7 +1017,10 @@ def checkExtentAndCellSize(cfg, inputFile, dem, fileType): raise FileExistsError(message) # write raster to file - outFile = IOf.writeResultToRaster(dem["header"], inputField["rasterData"], outFile, flip=True) + useCompression = cfg["EXPORTS"].getboolean("useCompression") + outFile = IOf.writeResultToRaster( + dem["header"], inputField["rasterData"], outFile, flip=True, useCompression=useCompression + ) log.info("Saved remeshed raster to %s" % outFile) returnStr = str(pathlib.Path("remeshedRasters", outFile.name)) remeshedFlag = "Yes" diff --git a/avaframe/in2Trans/rasterUtils.py b/avaframe/in2Trans/rasterUtils.py index 5e64a6ed2..a66d1b780 100644 --- a/avaframe/in2Trans/rasterUtils.py +++ b/avaframe/in2Trans/rasterUtils.py @@ -1,5 +1,5 @@ """ - Raster (ascii and tif) file reader and handler +Raster (ascii and tif) file reader and handler """ @@ -150,7 +150,7 @@ def isEqualASCheader(headerA, headerB): ) -def writeResultToRaster(header, resultArray, outFileName, flip=False): +def writeResultToRaster(header, resultArray, outFileName, flip=False, useCompression=True): """Write 2D array to a raster file with header and save to location of outFileName Parameters @@ -165,6 +165,9 @@ class with methods that give cellsize, nrows, ncols, xllcenter flip: boolean if True, flip the rows of the resultArray when writing. AF considers the first line in a data array to be the southernmost one. Some formats (e.g. tif) have the northernmost line first + useCompression: boolean + True if compression should be used on writing tiff files (lzw) + Returns ------- @@ -172,32 +175,32 @@ class with methods that give cellsize, nrows, ncols, xllcenter to file being written """ - if header["driver"] == "AAIGrid": - outFile = outFileName.parent / (outFileName.name + ".asc") - elif header["driver"] == "GTiff": - outFile = outFileName.parent / (outFileName.name + ".tif") - - # try: - rasterOut = rasterio.open( - outFile, - "w", - driver=header["driver"], - crs=header["crs"], - nodata=header["nodata_value"], - transform=header["transform"], - height=resultArray.shape[0], - width=resultArray.shape[1], - count=1, - dtype=resultArray.dtype, - # decimal_precision=3, - ) - if flip: - rasterOut.write(np.flipud(resultArray), 1) - else: - rasterOut.write(resultArray, 1) - rasterOut.close() - # except: - # log.error("could not write {} to {}".format(resultArray, outFileName)) + driver = header["driver"] + if driver not in ("AAIGrid", "GTiff"): + raise ValueError(f"Unsupported driver: {driver}") + + extMap = {"AAIGrid": ".asc", "GTiff": ".tif"} + outFile = outFileName.parent / (outFileName.name + extMap[driver]) + + commonKwargs = { + "driver": driver, + "crs": header["crs"], + "nodata": header["nodata_value"], + "transform": header["transform"], + "height": resultArray.shape[0], + "width": resultArray.shape[1], + "count": 1, + "dtype": resultArray.dtype, + } + + extraKwargs = {} + if useCompression and driver == "GTiff": + extraKwargs = {"compress": "lzw"} + + with rasterio.open(outFile, "w", **commonKwargs, **extraKwargs) as rasterOut: + data = np.flipud(resultArray) if flip else resultArray + rasterOut.write(data, 1) + return outFile diff --git a/avaframe/in3Utils/fileHandlerUtils.py b/avaframe/in3Utils/fileHandlerUtils.py index 3139cfcfb..1f217df5e 100644 --- a/avaframe/in3Utils/fileHandlerUtils.py +++ b/avaframe/in3Utils/fileHandlerUtils.py @@ -368,7 +368,7 @@ def splitTimeValueToArrayInterval(cfgValues, endTime): Returns -------- items : 1D numpy array - time step values as 1D numpy array + sorted time step values as 1D numpy array """ if ":" in cfgValues: diff --git a/avaframe/tests/test_com1DFA.py b/avaframe/tests/test_com1DFA.py index 6b1c46734..6ac6f78cc 100644 --- a/avaframe/tests/test_com1DFA.py +++ b/avaframe/tests/test_com1DFA.py @@ -1746,6 +1746,7 @@ def test_exportFields(tmp_path): cfg = configparser.ConfigParser() cfg["GENERAL"] = {"resType": "ppr|pft|FT"} cfg["REPORT"] = {"plotFields": "ppr|pft|pfv|pke"} + cfg["EXPORTS"] = {"useCompression": "True"} Tsave = [0, 10, 15, 25, 40] demHeader = {} demHeader["cellsize"] = 1 diff --git a/avaframe/tests/test_com1DFATools.py b/avaframe/tests/test_com1DFATools.py index 108535b79..d34df7867 100644 --- a/avaframe/tests/test_com1DFATools.py +++ b/avaframe/tests/test_com1DFATools.py @@ -119,7 +119,7 @@ def test_updateResCoeffFields(tmp_path): dem["originalHeader"]["transform"] = transform dem["originalHeader"]["crs"] = rasterio.crs.CRS() - fields = com1DFATools.updateResCoeffFields(fields, cfg["GENERAL"], 0.0, dem) + fields = com1DFATools.updateResCoeffFields(fields, cfg["GENERAL"]) assert fields["cResRaster"].any() == False @@ -135,7 +135,7 @@ def test_updateResCoeffFields(tmp_path): "FT": FT, } - fields = com1DFATools.updateResCoeffFields(fields, cfg["GENERAL"], 0.0, dem) + fields = com1DFATools.updateResCoeffFields(fields, cfg["GENERAL"]) assert fields["cResRaster"].any() == True assert fields["cResRaster"][0, 4] == 1 diff --git a/avaframe/tests/test_deriveParameterSet.py b/avaframe/tests/test_deriveParameterSet.py index fdd63e3b9..03f89f2db 100644 --- a/avaframe/tests/test_deriveParameterSet.py +++ b/avaframe/tests/test_deriveParameterSet.py @@ -608,6 +608,7 @@ def test_checkExtentAndCellSize(tmp_path): cfg = configparser.ConfigParser() cfg["GENERAL"] = {"resizeThreshold": 3.0, "meshCellSize": 1.0, "meshCellSizeThreshold": 0.0001} cfg["GENERAL"]["avalancheDir"] = str(testDir) + cfg["EXPORTS"] = {"useCompression": "True"} inDir = testDir / "Inputs" inDirR = inDir / "RASTERS" fU.makeADir(inDirR)