fire2a.cell2fire

  1#!python
  2# fmt: off
  3"""
  4Classes and methods auxiliary to using Cell2FireW and its QGIS integration
  5
  6Currently:
  7* Writes C2FW's firebreak specification (a.csv file) from a QGIS raster layer
  8
  9Processes C2FW plain text outputs 
 10* Fire Scars (Grids directories) into rasters and polygons
 11* Statistics (Intensity, HitRos, ... directories) into rasters
 12
 13Auxiliary:
 14* get_scars_files (all together but a bit slower than the next two methods)
 15* get_scars_indexed (part 1/2)
 16* group_scars (part 2/2)
 17
 18Sample Usage:
 19```bash
 20python -m fire2a.cell2fire -vvv --base-raster ../fuels.asc --authid EPSG:25831 --scar-sample Grids/Grids2/F
 21orestGrid01.csv --scar-raster scar_raster.tif --scar-poly scar_poly.shp --burn-prob burn_prob.tif --stat-sample RateOfSpread/ROSFile11.asc --stat-raster ros.tif --stat-summary ros_stats.tif
 22```
 23"""
 24# fmt: on
 25__author__ = "Fernando Badilla"
 26__revision__ = "$Format:%H$"
 27
 28import logging
 29import sys
 30from pathlib import Path
 31
 32from qgis.core import QgsRasterLayer
 33
 34from fire2a import setup_file
 35from fire2a.utils import count_header_lines, fprint, loadtxt_nodata
 36
 37logger = logging.getLogger(__name__)
 38"""captures the local logger"""
 39
 40NAME, FILEPATH = setup_file(name="cell2fire", filepath=Path("/home/fdo/source/fire2a-lib/src/fire2a"))
 41"""setups the file name and path for logger (set path manually when pasting into ipython or jupyter session)"""
 42
 43
 44def raster_layer_to_firebreak_csv(layer: QgsRasterLayer, firebreak_val: int = 1, output_file="firebreaks.csv") -> None:
 45    """Write a (Cell2) Fire Simulator (C2FW) firebreak csv file from a QGIS raster layer  
 46    Usage as cli argument `Cell2Fire --FirebreakCells firebreaks.csv ...`
 47
 48    Args:
 49    -   layer (QgsRasterLayer): A QGIS raster layer, default is the active layer
 50    -   firebreak_val (int): The value used to identify the firebreaks, default is 666
 51    -   output_file (str or Path): The path to the output csv file, default is firebreaks.csv
 52
 53    QGIS Desktop Example: Choose method A or B, draw with serval (2)
 54
 55    1.A. Use the 'Create constant raster layer' tool to create one with the same extent extent, crs and pixel size than the fuels raster layer. Recommended: 
 56    - constant value = 0
 57    - (Advanced Parameters) output raster data type = Byte
 58
 59    1.B. Use the 'raster calculator' with a base layer and 0 in the formula
 60
 61    2. Use 'Serval' plugin-tool to draw with the mouse the firebreaks on the new raster layer (values =1). Reload or save as to see changes.
 62
 63    QGIS Python Console Example:  
 64    ```
 65    layer = iface.activeLayer()  
 66    from fire2a.firebreaks import raster_layer_to_firebreak_csv
 67    raster_layer_to_firebreak_csv(layer)
 68    import processing
 69    processing.run("fire2a:cell2firesimulator", ...
 70    ```
 71    See also: https://fire2a.github.io/docs/docs/qgis-toolbox/c2f_firebreaks.html
 72    """  # fmt: skip
 73    from numpy import array as np_array
 74    from numpy import where as np_where
 75
 76    from fire2a.raster import get_rlayer_data, xy2id
 77
 78    width = layer.width()
 79    data = get_rlayer_data(layer)
 80
 81    # numpy is hh,ww indexing
 82    yy, xx = np_where(data == firebreak_val)
 83    ids = np_array([xy2id(x, y, width) for x, y in zip(xx, yy)])
 84    ids += 1
 85
 86    with open(output_file, "w") as f:
 87        f.write("Year,Ncell\n")
 88        f.write(f"1,{','.join(map(str,ids))}\n")
 89
 90
 91def get_scars_files(sample_file: Path):
 92    """Get sorted lists of (non-empty) files matching the pattern `root/parent(+any digit)/children(+any digit).(any extension)`
 93
 94    Normally used to read Cell2FireW scars `results/Grids/Grids*/ForestGrid*.csv`
 95
 96    Args:
 97    - sample_file (Path): A sample file to extract the name and extension, parent and root directories
 98
 99    Returns a tuple:
100    - bool: True if successful, False otherwise.  
101    - str: Error message, if any.  
102    - Path: `root` - all other paths are relative to this and must be used as root/parent or root/child.  
103    - list[Path]: `parents` - sorted list of parent directories.  
104    - list[int]: `parents_ids` - corresponding simulation ids of these parents.  
105    - list[list[Path]]: `children` - sorted list of lists of children files (grouped by simulation)
106    - list[list[int]]: `children_ids` - list of lists of corresponding period ids of each simulation
107
108    Sample Usage:
109    ```python
110    ret_val, msg, root, parent_dirs, parent_ids, children, children_ids = get_scars_files(Path(sample_file))
111    if not ret_val:
112        logger.error(msg)
113    ```
114    """  # fmt: skip
115    from re import search as re_search
116
117    ext = sample_file.suffix
118    if match := re_search(r"(\d+)$", sample_file.stem):
119        num = match.group()
120    else:
121        msg = f"sample_file: {sample_file} does not contain a number at the end"
122        return False, msg, None, None, None, None, None
123    file_name_wo_num = sample_file.stem[: -len(num)]
124    parent = sample_file.absolute().parent
125    root = sample_file.absolute().parent.parent
126    parent = parent.relative_to(root)
127    if match := re_search(r"(\d+)$", parent.name):
128        num = match.group()
129    else:
130        msg = f"sample_file:{sample_file} parent:{parent} does not contain a number at the end"
131        return False, msg, root, None, None, None, None
132
133    parent_wo_num = parent.name[: -len(num)]
134    parent_dirs = []
135    parent_ids = []
136    for par in root.glob(parent_wo_num + "[0-9]*"):
137        if par.is_dir():
138            par = par.relative_to(root)
139            parent_ids += [int(re_search(r"(\d+)$", par.name).group(0))]
140            parent_dirs += [par]
141    adict = dict(zip(parent_dirs, parent_ids))
142    parent_dirs.sort(key=lambda x: adict[x])
143    parent_ids.sort()
144
145    children = []
146    children_ids = []
147    for par in parent_dirs:
148        chl_files = []
149        chl_ids = []
150        for afile in (root / par).glob(file_name_wo_num + "[0-9]*" + ext):
151            if afile.is_file() and afile.stat().st_size > 0:
152                afile = afile.relative_to(root)
153                chl_ids += [int(re_search(r"(\d+)$", afile.stem).group(0))]
154                chl_files += [afile]
155        adict = dict(zip(chl_files, chl_ids))
156        chl_files.sort(key=lambda x: adict[x])
157        chl_ids.sort()
158        children += [chl_files]
159        children_ids += [chl_ids]
160
161    # msg = f"Got {len(parent_dirs)} parent directories with {sum([len(chl) for chl in children])} children files"
162    msg = ""
163    return True, msg, root, parent_dirs, parent_ids, children, children_ids
164
165
166def get_scars_indexed(sample_file: Path):
167    """Get a sorted list of files with the pattern `root/parent(+any digit)/children(+any digit).(any extension)`
168
169    Args:
170    - sample_file (Path): A sample file to extract the extension, children name (wo ending number),  parent (wo ending number) and root directory
171
172    Returns a tuple:
173    - return_value (bool): True if successful, False otherwise
174    - return_message (str): Debug/Error message if any
175    - root (Path): all paths are relative to this and must be used as root/file
176    - parent_wo_num (str): parent name without the ending number
177    - child_wo_num (str): children name without the ending number
178    - extension (str): file extension
179    - files (list[Path]): sorted list of (relative paths) files
180    - indexes (list[Tuple[int,int]]]): list of tuples of simulation and period ids
181
182    Sample Usage
183    ```python
184    retval, retmsg, root, parent_wo_num, child_wo_num, ext, files, indexes = get_scars_indexed(sample_file)
185    if not ret_val:
186        logger.error(msg)
187    ```
188    """
189    from os import sep
190    from re import findall as re_findall
191    from re import search as re_search
192
193    from numpy import array as np_array
194    from numpy import fromiter as np_fromiter
195
196    ext = sample_file.suffix
197    if match := re_search(r"(\d+)$", sample_file.stem):
198        num = match.group()
199    else:
200        msg = f"sample_file: {sample_file} does not contain a number at the end"
201        return False, msg, None, None, None, ext, None, None
202    child_wo_num = sample_file.stem[: -len(num)]
203    parent = sample_file.absolute().parent
204    root = sample_file.absolute().parent.parent
205    parent = parent.relative_to(root)
206    if match := re_search(r"(\d+)$", parent.name):
207        num = match.group()
208    else:
209        msg = f"sample_file:{sample_file} parent:{parent} does not contain a number at the end"
210        return False, msg, root, None, child_wo_num, None, None
211
212    parent_wo_num = parent.name[: -len(num)]
213
214    files = np_array(
215        [
216            afile.relative_to(root)
217            for afile in root.glob(parent_wo_num + "[0-9]*" + sep + child_wo_num + "[0-9]*" + ext)
218            if afile.is_file() and afile.stat().st_size > 0
219        ]
220    )
221
222    if sep == "\\":
223        sep = "\\\\"
224    indexes = np_fromiter(
225        # re_findall(parent_wo_num + r"(\d+)" + sep + child_wo_num + r"(\d+)" + ext, " ".join(map(str, files))),
226        re_findall(r"(\d+)" + sep + child_wo_num + r"(\d+)", " ".join(map(str, files))),
227        dtype=[("sim", int), ("per", int)],
228        count=len(files),
229    )
230
231    files = files[indexes.argsort(order=("sim", "per"))]
232    indexes.sort()
233
234    msg = ""
235    # msg = f"Got {len(files)} files\n"
236    # msg += f"{len(np_unique(indexes['sim']))} simulations\n"
237    # msg += f"{len(np_unique(indexes['per']))} different periods"
238
239    return True, msg, root, parent_wo_num, child_wo_num, ext, files, indexes
240
241
242def group_scars(root, parent_wo_num, child_wo_num, ext, files, indexes):
243    """Group scars files by simulation and period
244
245    Args:
246    - root (Path): root directory
247    - parent_wo_num (str): parent name without the ending number
248    - child_wo_num (str): children name without the ending number
249    - ext (str): file extension
250    - files (list[Path]): list of files
251    - indexes (list[Tuple[int,int]]): list of tuples of simulation and period ids
252
253    Returns:
254    - parent_ids (list[int]): list of simulation ids
255    - parent_dirs (list[Path]): list of parent directories
256    - children_ids (list[Tuple[int,int]]): list of tuples of simulation and period ids
257    - children_files (list[Path]): list of children files
258    - final_scars_ids (list[Tuple[int,int]]): list of tuples of simulation and period ids
259    - final_scars_files (list[Path]): list of final scars files
260
261    Sample:
262    ```python
263    retval, retmsg, root, parent_wo_num, child_wo_num, ext, files, indexes = get_scars_indexed(sample_file)
264    if not retval:
265        logger.error(retmsg)
266        sys.exit(1)
267    parent_ids, parent_dirs, children_ids, children_files, final_scars_ids, final_scars_files = group_scars(
268        root, parent_wo_num, child_wo_num, ext, files, indexes
269    )
270    ```
271    """
272    from numpy import unique as np_unique
273    from numpy import where as np_where
274
275    parent_ids = [sim_id for sim_id in np_unique(indexes["sim"])]
276    children_ids = [[(sim_id, per_id) for per_id in indexes["per"][indexes["sim"] == sim_id]] for sim_id in parent_ids]
277    children_files = [[afile for afile in files[np_where(indexes["sim"] == pid)[0]]] for pid in parent_ids]
278
279    final_idx = [np_where(indexes["sim"] == pid)[0][-1] for pid in parent_ids]
280
281    parent_dirs = [afile.parent for afile in files[final_idx]]
282
283    final_scars_files = files[final_idx]
284    final_scars_ids = indexes[final_idx]
285
286    return parent_ids, parent_dirs, children_ids, children_files, final_scars_ids, final_scars_files
287
288
289def build_scars(
290    scar_raster: str,
291    scar_poly: str,
292    burn_prob: str,
293    sample_file: Path,
294    W: int,
295    H: int,
296    geotransform: tuple,
297    authid: str,
298    callback=None,
299    feedback=None,
300):
301    """Builds the final scars raster, evolution scars polygons and/or burn probability raster files
302
303    Args:
304    - scar_raster (str): The output file name for the final scars raster
305    - scar_poly (str): The output file name for the evolution scars polygons
306    - burn_prob (str): The output file name for the burn probability raster
307    - sample_file (Path): Any scar sample file usually `results/Grids/Grids1/ForestGrid1.csv`
308    - W (int): Width of the raster
309    - H (int): Height of the raster
310    - geotransform (tuple): The geotransform of the raster
311    - authid (str): The projection of the raster
312    - callback (function): A function to call with the progress percentage
313    - feedback (QgsFeedback): A feedback object
314
315    Returns:
316    - int: 0 if successful, 1 otherwise
317    """
318    from numpy import any as np_any
319    from numpy import float32 as np_float32
320    from numpy import int8 as np_int8
321    from numpy import loadtxt as np_loadtxt
322    from numpy import zeros as np_zeros
323    from osgeo import gdal, ogr, osr
324
325    from fire2a.processing_utils import get_output_raster_format, get_vector_driver_from_filename
326
327    gdal.UseExceptions()
328
329    retval, retmsg, root, parent_wo_num, child_wo_num, ext, files, indexes = get_scars_indexed(sample_file)
330    if not retval:
331        fprint(retmsg, level="error", feedback=feedback, logger=logger)
332        return 1
333    parent_ids, parent_dirs, children_ids, children_files, final_scars_ids, final_scars_files = group_scars(
334        root, parent_wo_num, child_wo_num, ext, files, indexes
335    )
336
337    if burn_prob:
338        burn_prob_arr = np_zeros((H, W), dtype=np_float32)
339    else:
340        burn_prob_arr = None
341
342    if scar_raster:
343        driver_name = get_output_raster_format(scar_raster, feedback=feedback)
344        scar_raster_ds = gdal.GetDriverByName(driver_name).Create(
345            scar_raster, W, H, len(final_scars_ids), gdal.GDT_Byte
346        )
347        scar_raster_ds.SetGeoTransform(geotransform)
348        scar_raster_ds.SetProjection(authid)
349    else:
350        scar_raster_ds = None
351
352    def final_scar_step(i, data, afile, scar_raster, scar_raster_ds, burn_prob, burn_prob_arr, feedback=None):
353        if scar_raster:
354            band = scar_raster_ds.GetRasterBand(i)
355            # band.SetUnitType("burned")
356            if 0 != band.SetNoDataValue(0):
357                fprint(
358                    f"Set NoData failed for Final Scar {i}: {afile}", level="warning", feedback=feedback, logger=logger
359                )
360            if 0 != band.WriteArray(data):
361                fprint(
362                    f"WriteArray failed for Final Scar {i}: {afile}", level="warning", feedback=feedback, logger=logger
363                )
364            if i % 100 == 0:
365                scar_raster_ds.FlushCache()
366        if burn_prob:
367            if np_any(data == -1):
368                mask = data != -1
369                burn_prob_arr[ mask ] += data[ mask ]
370            else:
371                burn_prob_arr += data
372
373    if scar_poly:
374        # raster for each grid
375        src_ds = gdal.GetDriverByName("MEM").Create("", W, H, 1, gdal.GDT_Byte)
376        src_ds.SetGeoTransform(geotransform)
377        src_ds.SetProjection(authid)  # export coords to file
378
379        # datasource for shadow geometry vector layer (polygonize output)
380        ogr_ds = ogr.GetDriverByName("Memory").CreateDataSource("")
381
382        # spatial reference
383        sp_ref = osr.SpatialReference()
384        sp_ref.SetFromUserInput(authid)
385
386        if scar_poly.startswith("memory:"):
387            driver_name = "GPKG"
388            fprint(f"Memory layer, using {driver_name} driver", level="info", feedback=feedback, logger=logger)
389        else:
390            driver_name = get_vector_driver_from_filename(scar_poly)
391            fprint(f"NOT Mem lyr, using {driver_name} driver", level="info", feedback=feedback, logger=logger)
392
393        drv = ogr.GetDriverByName(driver_name)
394        if 0 != drv.DeleteDataSource(scar_poly):
395            fprint(f"Failed to delete {scar_poly}", level="error", feedback=feedback, logger=logger)
396        otrods = drv.CreateDataSource(scar_poly)
397        otrolyr = otrods.CreateLayer("propagation_scars", srs=sp_ref, geom_type=ogr.wkbPolygon)
398        otrolyr.CreateField(ogr.FieldDefn("simulation", ogr.OFTInteger))
399        otrolyr.CreateField(ogr.FieldDefn("time", ogr.OFTInteger))
400        otrolyr.CreateField(ogr.FieldDefn("area", ogr.OFTInteger))
401        otrolyr.CreateField(ogr.FieldDefn("perimeter", ogr.OFTInteger))
402
403        count_evo = 0
404        count_fin = 0
405        total_files = sum(len(files) for files in children_files)
406
407        for sim_id, ids, files in zip(parent_ids, children_ids, children_files):
408            count_fin += 1
409            for (_, per_id), afile in zip(ids, files):
410                count_evo += 1
411
412                try:
413                    data = np_loadtxt(root / afile, delimiter=",", dtype=np_int8)
414                except:
415                    fprint(
416                        f"Error reading {afile}, retrying with nodata = 0",
417                        level="error",
418                        feedback=feedback,
419                        logger=logger,
420                    )
421                    data = loadtxt_nodata(root / afile, delimiter=",", dtype=np_int8, no_data=0)
422
423                if not np_any(data == 1):
424                    fprint(
425                        f"no fire in {afile}, skipping propagation polygon",
426                        level="warning",
427                        feedback=feedback,
428                        logger=logger,
429                    )
430                else:
431                    src_band = src_ds.GetRasterBand(1)
432                    src_band.SetNoDataValue(0)
433                    src_band.WriteArray(data)
434
435                    ogr_layer = ogr_ds.CreateLayer("", srs=sp_ref)
436                    gdal.Polygonize(src_band, src_band, ogr_layer, -1, ["8CONNECTED=8"])
437
438                    feat = ogr_layer.GetNextFeature()
439                    geom = feat.GetGeometryRef()
440                    featureDefn = otrolyr.GetLayerDefn()
441                    feature = ogr.Feature(featureDefn)
442                    feature.SetGeometry(geom)
443                    feature.SetField("simulation", int(sim_id))
444                    feature.SetField("time", int(per_id))
445                    feature.SetField("area", int(geom.GetArea()))
446                    feature.SetField("perimeter", int(geom.Boundary().Length()))
447                    otrolyr.CreateFeature(feature)
448                    if count_evo % 100 == 0:
449                        otrods.FlushCache()
450
451                if callback:
452                    callback(count_evo / len(files) * 100, f"Processed Propagation-Scar {count_evo}/{len(indexes)}")
453                else:
454                    fprint(
455                        f"Processed Propagation-Scar {count_evo}/{len(indexes)}",
456                        level="info",
457                        feedback=feedback,
458                        logger=logger,
459                    )
460
461            if scar_raster or burn_prob:
462                final_scar_step(
463                    count_fin, data, afile, scar_raster, scar_raster_ds, burn_prob, burn_prob_arr, feedback=feedback
464                )
465                if callback:
466                    callback(None, f"Processed +Final-Scar {count_evo}/{len(indexes)}")
467                else:
468                    fprint(
469                        f"Processed Final-Scar {count_fin}/{len(files)}", level="info", feedback=feedback, logger=logger
470                    )
471        # clean up
472        if scar_raster:
473            scar_raster_ds.FlushCache()
474            scar_raster_ds = None
475        otrods.FlushCache()
476        otrods = None
477        # otrolyr.SyncToDisk() CRASHES QGIS
478        # otrolyr.FlushCache()
479        otrolyr = None
480        src_ds.FlushCache()
481        src_ds = None
482    else:
483        # final scar loop
484        count_fin = 0
485        for (sim_id, per_id), afile in zip(final_scars_ids, final_scars_files):
486            count_fin += 1
487            try:
488                data = np_loadtxt(root / afile, delimiter=",", dtype=np_int8)
489            except:
490                fprint(
491                    f"Error reading {afile}, retrying with nodata = 0", level="error", feedback=feedback, logger=logger
492                )
493                data = loadtxt_nodata(root / afile, delimiter=",", dtype=np_int8, no_data=0)
494            if not np_any(data == 1):
495                fprint(f"no fire in Final-Scar {afile}", level="warning", feedback=feedback, logger=logger)
496            final_scar_step(
497                count_fin, data, afile, scar_raster, scar_raster_ds, burn_prob, burn_prob_arr, feedback=feedback
498            )
499            if callback:
500                callback(count_fin / len(final_scars_files) * 100, f"Processed Final-Scar {count_fin}/{len(files)}")
501            else:
502                fprint(f"Processed Final-Scar {count_fin}/{len(files)}", level="info", feedback=feedback, logger=logger)
503        if scar_raster:
504            scar_raster_ds.FlushCache()
505            scar_raster_ds = None
506
507    if burn_prob:
508        driver_name = get_output_raster_format(burn_prob, feedback=feedback)
509        burn_prob_ds = gdal.GetDriverByName(driver_name).Create(burn_prob, W, H, 1, gdal.GDT_Float32)
510        burn_prob_ds.SetGeoTransform(geotransform)
511        burn_prob_ds.SetProjection(authid)
512        band = burn_prob_ds.GetRasterBand(1)
513        # band.SetUnitType("probability")
514        if 0 != band.SetNoDataValue(0):
515            fprint(
516                f"Set NoData failed for Burn Probability {burn_prob}", level="warning", feedback=feedback, logger=logger
517            )
518        if 0 != band.WriteArray(burn_prob_arr / len(final_scars_files)):  # type: ignore
519            fprint(
520                f"WriteArray failed for Burn Probability {burn_prob}", level="warning", feedback=feedback, logger=logger
521            )
522        burn_prob_ds.FlushCache()
523        burn_prob_ds = None
524
525    return 0
526
527
528def glob_numbered_files(sample_file: Path) -> tuple[list[Path], Path, str, str]:
529    """Get a list of files with the same name (+ any digit) and extension and the directory and name of the sample file
530
531    Args:
532    - sample_file (Path): A sample file to extract the extension, name and directory
533      - e.g., results/Intensity/Intensity1.asc, RateOfSpread/ROSFile2.asc, FlameLength/FL.asc, CrownFractionBurn/Cfb.asc, etc.
534
535    Returns a tuple:
536    - files (list[Path]): sorted list of files
537    - adir (Path): directory of the sample file
538    - aname (str): name of the sample file
539    - ext (str): extension of the sample file
540    """
541    from re import search
542
543    ext = sample_file.suffix
544    if match := search(r"(\d+)$", sample_file.stem):
545        num = match.group()
546    else:
547        raise ValueError(f"sample_file: {sample_file} does not contain a number at the end")
548    aname = sample_file.stem[: -len(num)]
549    adir = sample_file.absolute().parent
550    files = []
551    for afile in sorted(adir.glob(aname + "[0-9]*" + ext)):
552        if afile.is_file() and afile.stat().st_size > 0:
553            files += [afile]
554    # QgsMessageLog.logMessage(f"files: {files}, adir: {adir}, aname: {aname}, ext: {ext}", "fire2a", Qgis.Info)
555    return files, adir, aname, ext
556
557
558def build_stats(
559    stat_raster: str,
560    stat_summary: str,
561    sample_file: Path,
562    W: int,
563    H: int,
564    geotransform: tuple,
565    authid: str,
566    callback=None,
567    feedback=None,
568):
569    """Builds final statistics raster (1 band per-simulation) and summary raster (2 bands: mean against pixel burn count and stdev against total number of simulations) files
570    """
571    from numpy import float32 as np_float32
572    from numpy import loadtxt as np_loadtxt
573    from numpy import sqrt as np_sqrt
574    from numpy import zeros as np_zeros
575    from osgeo import gdal
576
577    from fire2a.processing_utils import get_output_raster_format
578
579    gdal.UseExceptions()
580
581    files, root, aname, ext = glob_numbered_files(sample_file)
582    num_headers = count_header_lines(files[0], sep=" ", feedback=feedback)
583
584    if stat_raster:
585        driver_name = get_output_raster_format(stat_raster, feedback=feedback)
586        stat_raster_ds = gdal.GetDriverByName(driver_name).Create(stat_raster, W, H, len(files), gdal.GDT_Float32)
587        stat_raster_ds.SetGeoTransform(geotransform)
588        stat_raster_ds.SetProjection(authid)
589    else:
590        stat_raster_ds = None
591
592    if stat_summary:
593        summed = np_zeros((H, W), dtype=np_float32)
594        sumsquared = np_zeros((H, W), dtype=np_float32)
595        burncount = np_zeros((H, W), dtype=np_float32)
596    else:
597        summed = None
598        sumsquared = None
599        burncount = None
600
601    count = 0
602    for afile in files:
603        count += 1
604        try:
605            data = np_loadtxt(root / afile, delimiter=" ", dtype=np_float32, skiprows=num_headers)
606        except Exception as e:
607            fprint(
608                f"Error reading {afile}, retrying with nodata = -9999: {e}",
609                level="error",
610                feedback=feedback,
611                logger=logger,
612            )
613            data = loadtxt_nodata(root / afile, delimiter=" ", dtype=np_float32, skiprows=num_headers, no_data=-9999)
614
615        if stat_raster_ds:
616            band = stat_raster_ds.GetRasterBand(count)
617            if 0 != band.SetNoDataValue(0):
618                fprint(
619                    f"Set NoData failed for Statistic {count}: {afile}",
620                    level="warning",
621                    feedback=feedback,
622                    logger=logger,
623                )
624            if 0 != band.WriteArray(data):
625                fprint(
626                    f"WriteArray failed for Statistic {count}: {afile}",
627                    level="warning",
628                    feedback=feedback,
629                    logger=logger,
630                )
631            if count % 100 == 0:
632                stat_raster_ds.FlushCache()
633
634        if stat_summary:
635            mask = data != -9999
636            tmp = data[mask]
637            summed[mask] += tmp
638            sumsquared[mask] += tmp ** 2
639            burncount[mask & (data != 0)] += 1
640
641        if callback:
642            callback(count / len(files) * 100, f"Processed Statistic {count}/{len(files)}")
643        else:
644            fprint(f"Processed Statistic {count}/{len(files)}", level="info", feedback=feedback, logger=logger)
645
646    if stat_raster_ds:
647        stat_raster_ds.FlushCache()
648        stat_raster_ds = None
649
650    if stat_summary:
651        driver_name = get_output_raster_format(stat_summary, feedback=feedback)
652        stat_summary_ds = gdal.GetDriverByName(driver_name).Create(stat_summary, W, H, 2, gdal.GDT_Float32)
653        stat_summary_ds.SetGeoTransform(geotransform)
654        stat_summary_ds.SetProjection(authid)
655        # mean
656        mask = burncount != 0
657        mean = np_zeros((H, W), dtype=np_float32) - 9999
658        mean[mask] = summed[mask] / burncount[mask]
659        band = stat_summary_ds.GetRasterBand(1)
660        if 0 != band.SetNoDataValue(-9999):
661            fprint(
662                f"Set NoData failed for Statistic {count}: {afile}", level="warning", feedback=feedback, logger=logger
663            )
664        if 0 != band.WriteArray(mean):
665            fprint(
666                f"WriteArray failed for Statistic {count}: {afile}", level="warning", feedback=feedback, logger=logger
667            )
668        # std
669        # from all simulations
670        N = len(files)
671        stddev = np_sqrt(sumsquared / N - (summed/N)**2)
672        # beacause this is always zero:
673        # stddev = np_zeros((H, W), dtype=np_float32) - 9999
674        # zero_mask = burncount == 0
675        # burncount[zero_mask] = 1
676        # stddev = np_sqrt(sumsquared / burncount - mean ** 2)
677        # stddev[~mask] = -9999
678        band = stat_summary_ds.GetRasterBand(2)
679        if 0 != band.SetNoDataValue(-9999):
680            fprint(
681                f"Set NoData failed for Statistic {count}: {afile}", level="warning", feedback=feedback, logger=logger
682            )
683        if 0 != band.WriteArray(stddev):
684            fprint(
685                f"WriteArray failed for Statistic {count}: {afile}", level="warning", feedback=feedback, logger=logger
686            )
687        stat_summary_ds.FlushCache()
688        stat_summary_ds = None
689
690    return 0
691
692
693def arg_parser(argv):
694    """Arguments parses to coordinate: scar and stats processing, verbosity, logfile
695
696    Args:
697    - argv (list): list of strings
698      - e.g., interactive: `args = arg_parser(['-vvv'])`
699
700    Returns:
701    - argparse.Namespace: parsed arguments
702    """
703    import argparse
704
705    parser = argparse.ArgumentParser(
706        description=r"Fire2a-Cell2FireW Related algorithms CLI. Implemented here: Cell2FireW plain scars and statistic outputs into rasters and polygons.",
707        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
708    )
709
710    scars = parser.add_argument_group(
711        "Scars: Transform C2FW scar outputs into rasters and polygons, i.e.: results/Grids/Grids*/ForestGrid*.csv"
712    )
713    scars.add_argument(
714        "--scar-sample",
715        help="Matching the pattern 'root/parent(+any digit)/children(+any digit).(any csv extension)' ",
716    )
717    scars.add_argument(
718        "--scar-raster", default="", help="Output file name for the final scars raster, each band a simulation"
719    )
720    scars.add_argument(
721        "--scar-poly",
722        default="",
723        help="Output file name for the evolution scars polygons, multiple features with simulation and period attributes",
724    )
725    scars.add_argument("--burn-prob", default="", help="Output file name for the burn probability raster, one band")
726
727    stats = parser.add_argument_group(
728        "Stats: Transform C2FW outputs into rasters, .e.g.: results/Intensity/Intensity1.asc, RateOfSpread/ROSFile2.asc, FlameLength/FL.asc, CrownFractionBurn/Cfb.asc, etc."
729    )
730    stats.add_argument(
731        "--stat-sample",
732        default="",
733        help="Matching the pattern 'statistic_name(+any digit).asc' ",
734    )
735    stats.add_argument(
736        "--stat-raster", default="", help="Output file name for the statistic raster, each band a simulation"
737    )
738    stats.add_argument(
739        "--stat-summary", default="", help="Output file name for the raster, two bands: mean and std-deviation"
740    )
741
742    parser.add_argument("--verbose", "-v", action="count", default=0, help="WARNING:1, INFO:2, DEBUG:3")
743    parser.add_argument(
744        "--logfile",
745        "-l",
746        action="store_true",
747        help="enable 5 log files named " + NAME + ".log (verbose must be enabled)",
748        default=None,
749    )
750    parser.add_argument(
751        "--base-raster", required=True, help="Raster to base the geotransform, width, heigth and authid"
752    )
753    parser.add_argument(
754        "--authid", required=False, help="Auth id to override (or missing on the base raster (e.g., using .asc raster)"
755    )
756    args = parser.parse_args(argv)
757    if args.logfile:
758        args.logfile = NAME + ".log"
759    return args
760
761
762def main(argv=None):
763    """
764
765    args = arg_parser(['-vvv'])
766    """
767
768    if argv is sys.argv:
769        argv = sys.argv[1:]
770    args = arg_parser(argv)
771
772    if args.verbose != 0:
773        global logger
774        from fire2a import setup_logger
775
776        logger = setup_logger(verbosity=args.verbose, logfile=args.logfile)
777        # set other modules logging level
778        logging.getLogger("asyncio").setLevel(logging.INFO)
779
780    logger.info("args %s", args)
781
782    if not args.base_raster:
783        logger.error("No base raster file name provided")
784        return 1
785
786    if not Path(args.base_raster).is_file():
787        logger.error("Base raster %s is not a file", args.base_raster)
788        return 1
789
790    from fire2a.raster import read_raster
791
792    _, raster_props = read_raster(args.base_raster, data=False, info=True)
793
794    if not (authid := raster_props["Projection"]):
795        if not (authid := args.authid):
796            logger.error("No authid found on the base raster or provided")
797            return 1
798    logger.info("Read base raster, using authid: %s", authid)
799
800    retval = {}
801    if args.scar_sample and (args.scar_poly or args.scar_raster or args.burn_prob):
802        retval["scar"] = build_scars(
803            args.scar_raster,
804            args.scar_poly,
805            args.burn_prob,
806            Path(args.scar_sample),
807            raster_props["RasterXSize"],
808            raster_props["RasterYSize"],
809            raster_props["Transform"],
810            authid,
811        )
812        logger.info("built scars return value: %s (0 means sucess)", retval["scar"])
813
814    if args.stat_sample and args.stat_raster:
815        retval["stat"] = build_stats(
816            args.stat_raster,
817            args.stat_summary,
818            Path(args.stat_sample),
819            raster_props["RasterXSize"],
820            raster_props["RasterYSize"],
821            raster_props["Transform"],
822            authid,
823        )
824        logger.info("built statistic return value: %s (0 means success)", retval["stat"])
825
826    print(retval)
827    return sum(retval.values())
828
829
830if __name__ == "__main__":
831    sys.exit(main(sys.argv))