"""How to represent points clouds, specifying specific data signatures (e.g. choice of features).

Country-specific definitions inherit from abstract class LidarDataLogic, which implements general logics
to split each point cloud into subtiles, format these subtiles, and save a learning-ready dataset splitted into
train, val. A test set is also created, which is simply a copy of the selected test point clouds. This is so
test phase conditions are similar to "in-the-wild" prediction conditions.

In particular, subclasses implement a "load_las" method that is used by the datamodule at test and inference time.

"""

from abc import ABC, abstractmethod
import argparse
import os
import glob
import os.path as osp
from shutil import copyfile
import shutil
from tqdm import tqdm
import laspy
import numpy as np
import pandas as pd
import torch
from torch_geometric.data import Data

SPLIT_PHASES = ["val", "train", "test"]

# Files for the creation of a toy dataset.
LAS_SUBSET_FOR_TOY_DATASET = (
    "tests/data/toy_dataset_src/870000_6618000.subset.50mx100m.las"
)
SPLIT_CSV_FOR_TOY_DATASET = "tests/data/toy_dataset_src/toy_dataset_split.csv"


class LidarDataLogic(ABC):
    """Abstract class to load, chunk, and save a point cloud dataset according to a train/val/test split.
    load_las and its needed parameters ares specified in child classes.

    """

    return_number_normalization_max_value = 7

    def __init__(
        self,
        input_data_dir: str = None,
        prepared_data_dir: str = None,
        split_csv: str = None,
        input_tile_width_meters: int = 1000,
        subtile_width_meters: int = 50,
        **kwargs,
    ):
        self.input_data_dir = input_data_dir
        self.prepared_data_dir = prepared_data_dir
        self.split_csv = split_csv
        self.input_tile_width_meters = input_tile_width_meters
        self.subtile_width_meters = subtile_width_meters

    @abstractmethod
    def load_las(self, las_filepath: str) -> Data:
        """Load a point cloud in LAS format to memory and turn it into torch-geometric Data object.

        Args:
            las_filepath (str): path to the LAS file.

        Returns:
            Data: The point cloud formatted for later deep learning training.

        """
        raise NotImplementedError

    def prepare(self):
        """Prepare a dataset for model training and model evaluation.

        Iterates through LAS files listed in a csv metadata file.
        A `split` column specifies the train/val/test split of the dataset to be created.

        """
        split_df = pd.read_csv(self.split_csv)
        for phase in tqdm(SPLIT_PHASES, desc="Phases"):
            basenames = split_df[split_df.split == phase].basename.tolist()
            print(f"Subset: {phase}")
            print("  -  ".join(basenames))
            for file_basename in tqdm(basenames, desc="Files"):
                filepath = self._find_file_in_dir(self.input_data_dir, file_basename)
                output_subdir_path = osp.join(self.prepared_data_dir, phase)
                if phase == "test":
                    # Simply copy the LAS to the new test folder.
                    os.makedirs(output_subdir_path, exist_ok=True)
                    target_file = osp.join(output_subdir_path, file_basename)
                    copyfile(filepath, target_file)
                elif phase in ["train", "val"]:
                    # Load LAS into memory as a Data object with selected features,
                    # then iteratively extract 50m*50m subtiles by filtering along x
                    # then y axis. Serialize the resulting Data object using torch.save.
                    output_subdir_path = osp.join(
                        output_subdir_path, osp.basename(filepath)
                    )
                    os.makedirs(output_subdir_path, exist_ok=True)
                    self.split_and_save(filepath, output_subdir_path)
                else:
                    raise KeyError("Phase should be one of train/val/test.")

    def split_and_save(self, filepath: str, output_subdir_path: str) -> None:
        """Parse a LAS, extract and save each subtile as a Data object.

        Args:
            filepath (str): input LAS file
            output_subdir_path (str): output directory to save splitted `.data` objects.

        """
        range_by_axis = np.arange(
            self.input_tile_width_meters // self.subtile_width_meters + 1
        )

        data = self.load_las(filepath)
        idx = 0
        for _ in tqdm(range_by_axis):
            if not len(data.pos):
                # Ignore if empty
                break
            data_x_band = self._extract_by_x(data)
            for _ in range_by_axis:
                if not len(data_x_band.pos):
                    # Ignore if empty
                    break
                subtile_data = self._extract_by_y(data_x_band)
                self._save(subtile_data, output_subdir_path, idx)
                idx += 1

    def _find_file_in_dir(self, input_data_dir: str, basename: str) -> str:
        """Query files matching a basename in input_data_dir and its subdirectories.

        Args:
            input_data_dir (str): data directory

        Returns:
            [str]: first file path matching the query.

        """
        query = f"{input_data_dir}/**/{basename}"
        files = glob.glob(query, recursive=True)
        return files[0]

    def _extract_by_axis(self, data: Data, axis=0) -> Data:
        """Filter a data object on a chosen axis, using a relative position .
        Modifies the original data object so that extracted future filters are faster.

        Args:
            data (Data): a pyg Data object with pos, x, and y attributes.
            relative_pos (int): where the data to extract start on chosen axis (typically in range 0-1000)
            axis (int, optional): 0 for x and 1 for y axis. Defaults to 0.

        Returns:
            Data: the data that is at most subtile_width_meters above relative_pos on the chosen axis.
        """
        sub_tile_data = data.clone()
        pos_axis = sub_tile_data.pos[:, axis]
        absolute_low = pos_axis.min(0)
        absolute_high = absolute_low + self.subtile_width_meters
        mask = (absolute_low <= pos_axis) & (pos_axis <= absolute_high)

        # select
        sub_tile_data.pos = sub_tile_data.pos[mask]
        sub_tile_data.x = sub_tile_data.x[mask]
        sub_tile_data.y = sub_tile_data.y[mask]

        data.pos = data.pos[~mask]
        data.x = data.x[~mask]
        data.y = data.y[~mask]
        return sub_tile_data

    def _extract_by_x(self, data: Data) -> Data:
        """extract_by_axis applied on first axis x"""
        return self._extract_by_axis(data, axis=0)

    def _extract_by_y(self, data: Data) -> Data:
        """extract_by_axis applied on second axis y"""
        return self._extract_by_axis(data, axis=1)

    def _save(self, subtile_data: Data, output_subdir_path: str, idx: int) -> None:
        """Save the subtile data object with torch.

        Args:
            subtile_data (Data): the object to save.
            output_subdir_path (str): the subfolder to save it.
            idx (int): an arbitrary but unique subtile identifier.
        """
        subtile_save_path = osp.join(output_subdir_path, f"{str(idx).zfill(4)}.data")
        torch.save(subtile_data, subtile_save_path)


class FrenchLidarDataLogic(LidarDataLogic):

    x_features_names = [
        "intensity",
        "return_number",
        "number_of_returns",
        "red",
        "green",
        "blue",
        "nir",
        "rgb_avg",
        "ndvi",
    ]
    colors_normalization_max_value = 255 * 256

    # Exclusion of 65: artefacts, and 66: virtual points.
    excluded_classes = [65, 66]

    @classmethod
    def load_las(self, las_filepath: str):
        f"""Loads a point cloud in LAS format to memory and turns it into torch-geometric Data object.

        Builds a composite (average) color channel on the fly.

        Calculate NDVI on the fly.

        x_features_names are {' - '.join(self.x_features_names)}

        Args:
            las_filepath (str): path to the LAS file.

        Returns:
            Data: the point cloud formatted for later deep learning training.

        """
        las = laspy.read(las_filepath)

        # Filter out unused classes
        las.points = las.points[~np.isin(las.classification, self.excluded_classes)]
        # Positions and base features
        pos = np.asarray([las.x, las.y, las.z], dtype=np.float32).transpose()
        x = np.asarray(
            [
                las[x_name]
                for x_name in [
                    "intensity",
                    "return_number",
                    "number_of_returns",
                    "red",
                    "green",
                    "blue",
                    "nir",
                ]
            ],
            dtype=np.float32,
        ).transpose()

        # normalization
        return_number_idx = self.x_features_names.index("return_number")
        occluded_points = x[:, return_number_idx] > 1

        x[:, return_number_idx] = (x[:, return_number_idx]) / (
            self.return_number_normalization_max_value
        )
        num_return_idx = self.x_features_names.index("number_of_returns")
        x[:, num_return_idx] = (x[:, num_return_idx]) / (
            self.return_number_normalization_max_value
        )

        for idx, c in enumerate(self.x_features_names):
            if c in ["red", "green", "blue", "nir"]:
                print(
                    x[:, idx].max()
                )  # DEBUG: just to be sure that it is the same as before
                assert x[:, idx].max() <= self.colors_normalization_max_value
                x[:, idx] = x[:, idx] / self.colors_normalization_max_value
                x[occluded_points, idx] = 0

        red = x[:, self.x_features_names.index("red")]
        green = x[:, self.x_features_names.index("green")]
        blue = x[:, self.x_features_names.index("blue")]

        # Additional features :
        # Average color, that will be normalized on the fly based on single-sample
        rgb_avg = np.asarray([red, green, blue], dtype=np.float32).mean(axis=0)

        # NDVI
        nir = x[:, self.x_features_names.index("nir")]
        ndvi = (nir - red) / (nir + red + 10**-6)
        x = np.concatenate([x, rgb_avg[:, None], ndvi[:, None]], axis=1)

        try:
            # for LAS format V1.2
            y = las.classification.array.astype(int)
        except Exception:
            # for  LAS format V1.4
            y = las.classification.astype(int)

        return Data(
            pos=pos,
            x=x,
            y=y,
            las_filepath=las_filepath,
            x_features_names=self.x_features_names,
        )


class SwissTopoLidarDataLogic(LidarDataLogic):
    x_features_names = [
        "intensity",
        "return_number",
        "number_of_returns",
        "red",
        "green",
        "blue",
        "rgb_avg",
    ]
    colors_normalization_max_value = 256

    @classmethod
    def load_las(self, las_filepath: str) -> Data:
        """Loads a point cloud in LAS format to memory and turns it into torch-geometric Data object.

        Builds a composite (average) color channel on the fly.

        Args:
            las_filepath (str): path to the LAS file.

        Returns:
            Data: the point cloud formatted for later deep learning training.

        """
        las = laspy.read(las_filepath)
        pos = np.asarray([las.x, las.y, las.z], dtype=np.float32).transpose()

        x = np.asarray(
            [
                las[x_name]
                for x_name in [
                    "intensity",
                    "return_number",
                    "number_of_returns",
                    "red",
                    "green",
                    "blue",
                ]
            ],
            dtype=np.float32,
        ).transpose()

        return_number_idx = self.x_features_names.index("return_number")
        occluded_points = x[:, return_number_idx] > 1

        x[:, return_number_idx] = (x[:, return_number_idx]) / (
            self.return_number_normalization_max_value
        )
        num_return_idx = self.x_features_names.index("number_of_returns")
        x[:, num_return_idx] = (x[:, num_return_idx]) / (
            self.return_number_normalization_max_value
        )

        for idx, c in enumerate(self.x_features_names):
            if c in ["red", "green", "blue"]:
                assert x[:, idx].max() <= self.colors_normalization_max_value
                x[:, idx] = x[:, idx] / self.colors_normalization_max_value
                x[occluded_points, idx] = 0

        rgb_avg = (
            np.asarray(
                [las[x_name] for x_name in ["red", "green", "blue"]], dtype=np.float32
            )
            .transpose()
            .mean(axis=1, keepdims=True)
        )

        x = np.concatenate([x, rgb_avg], axis=1)

        try:
            # for LAS format V1.2
            y = las.classification.array.astype(int)
        except Exception:
            # for  LAS format V1.4
            y = las.classification.astype(int)

        return Data(
            pos=pos,
            x=x,
            y=y,
            las_filepath=las_filepath,
            x_features_names=self.x_features_names,
        )


def make_toy_dataset_from_test_file(
    prepared_data_dir: str,
    src_las_path: str = LAS_SUBSET_FOR_TOY_DATASET,
    split_csv: str = SPLIT_CSV_FOR_TOY_DATASET,
):
    """Prepare a toy dataset from a single, small LAS file.

    The file is first duplicated to get 2 LAS in each split (train/val/test),
    and then each file is splitted into .data files, resulting in a training-ready
    dataset loacted in td_prepared

    Args:
        src_las_path (str): input, small LAS file to generate toy dataset from
        split_csv (str): Path to csv with a `basename` (e.g. '123_456.las') and
        a `split` (train/val/test) columns specifying the dataset split.
        prepared_data_dir (str): where to copy files (`raw` subfolder) and to prepare
        dataset files (`prepared` subfolder)

    Returns:
        str: path to directory containing prepared dataset.

    """
    # Copy input file for full test isolation
    td_raw = osp.join(prepared_data_dir, "raw")
    td_prepared = osp.join(prepared_data_dir, "prepared")
    os.makedirs(td_raw)
    os.makedirs(td_prepared)
    # Make a "raw", unporcessed dataset with six files.
    basename = osp.basename(src_las_path)
    for s in ["train1", "train2", "val1", "val2", "test1", "test2"]:
        copy_path = osp.join(td_raw, basename.replace(".las", f".{s}.las"))
        shutil.copy(src_las_path, copy_path)

    # Prepare a Deep-Learning-ready dataset, using the split defined in the csv.
    data_prepper = FrenchLidarDataLogic(
        input_data_dir=td_raw,
        prepared_data_dir=td_prepared,
        split_csv=split_csv,
        input_tile_width_meters=110,
        subtile_width_meters=50,
    )
    data_prepper.prepare()
    return td_prepared


def _get_data_preparation_parser():
    """Gets a parser with parameters for dataset preparation.

    Returns:
        argparse.ArgumentParser: the parser.
    """
    parser = argparse.ArgumentParser(
        description="Prepare a Lidar dataset for deep learning.",
        epilog="If you need a toy dataset, you only need to"
        " set --origin=FR_TOY and specify a value for --prepared_data_dir.",
    )
    parser.add_argument(
        "--input_data_dir",
        type=str,
        default="./data/raw/",
        help="Path to directory with las files with unique basename and `.las` suffix.",
    )
    parser.add_argument(
        "--split_csv",
        type=str,
        default="./split.csv",
        help="Path to csv with a basename (e.g. '123_456.las') and split (train/val/test) columns specifying the dataset split.",
    )
    parser.add_argument(
        "--prepared_data_dir",
        type=str,
        default="./data/prepared/",
        help="Path to folder to save Data object train/val/test subfolders.",
    )
    parser.add_argument(
        "--origin", type=str, default="FR", choices=["FR", "CH", "FR_TOY"]
    )

    return parser


def main():
    """Code to execute to prepare a new set of Lidar tiles for training 3D segmentaiton neural networks.

    From a data directory containing point cloud in LAS format, and a scv specifying the dataset
    train/val/test split for each file (columns: split, basename, example: "val","123_456.las"),
    split the dataset, chunk the point cloud tiles into smaller subtiles, and prepare each sample
    as a pytorch geometric Data object.

    Echo numbers and colors are scaled to be in 0-1 range. Intensity and average color
    are not scaled b/c they are expected to be standardized later.

    To show help relative to data preparation, run:

        cd myria3d/datamodules/

        conda activate myria3d

        python loading.py -h


    """

    parser = _get_data_preparation_parser()
    args = parser.parse_args()
    if args.origin == "FR_TOY":
        make_toy_dataset_from_test_file(args.prepared_data_dir)
    elif args.origin == "FR":
        data_prepper = FrenchLidarDataLogic(**args.__dict__)
        data_prepper.prepare()
    elif args.origin == "CH":
        data_prepper = SwissTopoLidarDataLogic(**args.__dict__)
        data_prepper.prepare()
    else:
        raise KeyError(f"Data origin is invalid (currently: {args.origin})")


if __name__ == "__main__":
    main()

"""How to represent points clouds, specifying specific data signatures (e.g. choice of features).

Country-specific definitions inherit from abstract class LidarDataLogic, which implements general logics
to split each point cloud into subtiles, format these subtiles, and save a learning-ready dataset splitted into
train, val. A test set is also created, which is simply a copy of the selected test point clouds. This is so
test phase conditions are similar to "in-the-wild" prediction conditions.

In particular, subclasses implement a "load_las" method that is used by the datamodule at test and inference time.

"""

from abc import ABC, abstractmethod
import argparse
import os
import glob
import os.path as osp
from shutil import copyfile
import shutil
from tqdm import tqdm
import laspy
import numpy as np
import pandas as pd
import torch
from torch_geometric.data import Data

SPLIT_PHASES = ["val", "train", "test"]

# Files for the creation of a toy dataset.
LAS_SUBSET_FOR_TOY_DATASET = (
    "tests/data/toy_dataset_src/870000_6618000.subset.50mx100m.las"
)
SPLIT_CSV_FOR_TOY_DATASET = "tests/data/toy_dataset_src/toy_dataset_split.csv"


class LidarDataLogic(ABC):
    """Abstract class to load, chunk, and save a point cloud dataset according to a train/val/test split.
    load_las and its needed parameters ares specified in child classes.

    """

    return_number_normalization_max_value = 7

    def __init__(
        self,
        input_data_dir: str = None,
        prepared_data_dir: str = None,
        split_csv: str = None,
        input_tile_width_meters: int = 1000,
        subtile_width_meters: int = 50,
        **kwargs,
    ):
        self.input_data_dir = input_data_dir
        self.prepared_data_dir = prepared_data_dir
        self.split_csv = split_csv
        self.input_tile_width_meters = input_tile_width_meters
        self.subtile_width_meters = subtile_width_meters

    @abstractmethod
    def load_las(self, las_filepath: str) -> Data:
        """Load a point cloud in LAS format to memory and turn it into torch-geometric Data object.

        Args:
            las_filepath (str): path to the LAS file.

        Returns:
            Data: The point cloud formatted for later deep learning training.

        """
        raise NotImplementedError

    def prepare(self):
        """Prepare a dataset for model training and model evaluation.

        Iterates through LAS files listed in a csv metadata file.
        A `split` column specifies the train/val/test split of the dataset to be created.

        """
        split_df = pd.read_csv(self.split_csv)
        for phase in tqdm(SPLIT_PHASES, desc="Phases"):
            basenames = split_df[split_df.split == phase].basename.tolist()
            print(f"Subset: {phase}")
            print("  -  ".join(basenames))
            for file_basename in tqdm(basenames, desc="Files"):
                filepath = self._find_file_in_dir(self.input_data_dir, file_basename)
                output_subdir_path = osp.join(self.prepared_data_dir, phase)
                if phase == "test":
                    # Simply copy the LAS to the new test folder.
                    os.makedirs(output_subdir_path, exist_ok=True)
                    target_file = osp.join(output_subdir_path, file_basename)
                    copyfile(filepath, target_file)
                elif phase in ["train", "val"]:
                    # Load LAS into memory as a Data object with selected features,
                    # then iteratively extract 50m*50m subtiles by filtering along x
                    # then y axis. Serialize the resulting Data object using torch.save.
                    output_subdir_path = osp.join(
                        output_subdir_path, osp.basename(filepath)
                    )
                    os.makedirs(output_subdir_path, exist_ok=True)
                    self.split_and_save(filepath, output_subdir_path)
                else:
                    raise KeyError("Phase should be one of train/val/test.")

    def split_and_save(self, filepath: str, output_subdir_path: str) -> None:
        """Parse a LAS, extract and save each subtile as a Data object.

        Args:
            filepath (str): input LAS file
            output_subdir_path (str): output directory to save splitted `.data` objects.

        """
        range_by_axis = np.arange(
            self.input_tile_width_meters // self.subtile_width_meters + 1
        )

        data = self.load_las(filepath)
        idx = 0
        for _ in tqdm(range_by_axis):
            if not len(data.pos):
                # Ignore if empty
                break
            data_x_band = self._extract_by_x(data)
            for _ in range_by_axis:
                if not len(data_x_band.pos):
                    # Ignore if empty
                    break
                subtile_data = self._extract_by_y(data_x_band)
                self._save(subtile_data, output_subdir_path, idx)
                idx += 1

    def _find_file_in_dir(self, input_data_dir: str, basename: str) -> str:
        """Query files matching a basename in input_data_dir and its subdirectories.

        Args:
            input_data_dir (str): data directory

        Returns:
            [str]: first file path matching the query.

        """
        query = f"{input_data_dir}/**/{basename}"
        files = glob.glob(query, recursive=True)
        return files[0]

    def _extract_by_axis(self, data: Data, axis=0) -> Data:
        """Filter a data object on a chosen axis, using a relative position .
        Modifies the original data object so that extracted future filters are faster.

        Args:
            data (Data): a pyg Data object with pos, x, and y attributes.
            relative_pos (int): where the data to extract start on chosen axis (typically in range 0-1000)
            axis (int, optional): 0 for x and 1 for y axis. Defaults to 0.

        Returns:
            Data: the data that is at most subtile_width_meters above relative_pos on the chosen axis.
        """
        sub_tile_data = data.clone()
        pos_axis = sub_tile_data.pos[:, axis]
        absolute_low = pos_axis.min(0)
        absolute_high = absolute_low + self.subtile_width_meters
        mask = (absolute_low <= pos_axis) & (pos_axis <= absolute_high)

        # select
        sub_tile_data.pos = sub_tile_data.pos[mask]
        sub_tile_data.x = sub_tile_data.x[mask]
        sub_tile_data.y = sub_tile_data.y[mask]

        data.pos = data.pos[~mask]
        data.x = data.x[~mask]
        data.y = data.y[~mask]
        return sub_tile_data

    def _extract_by_x(self, data: Data) -> Data:
        """extract_by_axis applied on first axis x"""
        return self._extract_by_axis(data, axis=0)

    def _extract_by_y(self, data: Data) -> Data:
        """extract_by_axis applied on second axis y"""
        return self._extract_by_axis(data, axis=1)

    def _save(self, subtile_data: Data, output_subdir_path: str, idx: int) -> None:
        """Save the subtile data object with torch.

        Args:
            subtile_data (Data): the object to save.
            output_subdir_path (str): the subfolder to save it.
            idx (int): an arbitrary but unique subtile identifier.
        """
        subtile_save_path = osp.join(output_subdir_path, f"{str(idx).zfill(4)}.data")
        torch.save(subtile_data, subtile_save_path)


class FrenchLidarDataLogic(LidarDataLogic):

    x_features_names = [
        "intensity",
        "return_number",
        "number_of_returns",
        "red",
        "green",
        "blue",
        "nir",
        "rgb_avg",
        "ndvi",
    ]
    colors_normalization_max_value = 255 * 256

    # Exclusion of 65: artefacts, and 66: virtual points.
    excluded_classes = [65, 66]

    @classmethod
    def load_las(self, las_filepath: str):
        f"""Loads a point cloud in LAS format to memory and turns it into torch-geometric Data object.

        Builds a composite (average) color channel on the fly.

        Calculate NDVI on the fly.

        x_features_names are {' - '.join(self.x_features_names)}

        Args:
            las_filepath (str): path to the LAS file.

        Returns:
            Data: the point cloud formatted for later deep learning training.

        """
        las = laspy.read(las_filepath)

        # Filter out unused classes
        las.points = las.points[~np.isin(las.classification, self.excluded_classes)]
        # Positions and base features
        pos = np.asarray([las.x, las.y, las.z], dtype=np.float32).transpose()
        x = np.asarray(
            [
                las[x_name]
                for x_name in [
                    "intensity",
                    "return_number",
                    "number_of_returns",
                    "red",
                    "green",
                    "blue",
                    "nir",
                ]
            ],
            dtype=np.float32,
        ).transpose()

        # normalization
        return_number_idx = self.x_features_names.index("return_number")
        occluded_points = x[:, return_number_idx] > 1

        x[:, return_number_idx] = (x[:, return_number_idx]) / (
            self.return_number_normalization_max_value
        )
        num_return_idx = self.x_features_names.index("number_of_returns")
        x[:, num_return_idx] = (x[:, num_return_idx]) / (
            self.return_number_normalization_max_value
        )

        for idx, c in enumerate(self.x_features_names):
            if c in ["red", "green", "blue", "nir"]:
                print(
                    x[:, idx].max()
                )  # DEBUG: just to be sure that it is the same as before
                assert x[:, idx].max() <= self.colors_normalization_max_value
                x[:, idx] = x[:, idx] / self.colors_normalization_max_value
                x[occluded_points, idx] = 0

        red = x[:, self.x_features_names.index("red")]
        green = x[:, self.x_features_names.index("green")]
        blue = x[:, self.x_features_names.index("blue")]

        # Additional features :
        # Average color, that will be normalized on the fly based on single-sample
        rgb_avg = np.asarray([red, green, blue], dtype=np.float32).mean(axis=0)

        # NDVI
        nir = x[:, self.x_features_names.index("nir")]
        ndvi = (nir - red) / (nir + red + 10**-6)
        x = np.concatenate([x, rgb_avg[:, None], ndvi[:, None]], axis=1)

        try:
            # for LAS format V1.2
            y = las.classification.array.astype(int)
        except Exception:
            # for  LAS format V1.4
            y = las.classification.astype(int)

        return Data(
            pos=pos,
            x=x,
            y=y,
            las_filepath=las_filepath,
            x_features_names=self.x_features_names,
        )


class SwissTopoLidarDataLogic(LidarDataLogic):
    x_features_names = [
        "intensity",
        "return_number",
        "number_of_returns",
        "red",
        "green",
        "blue",
        "rgb_avg",
    ]
    colors_normalization_max_value = 256

    @classmethod
    def load_las(self, las_filepath: str) -> Data:
        """Loads a point cloud in LAS format to memory and turns it into torch-geometric Data object.

        Builds a composite (average) color channel on the fly.

        Args:
            las_filepath (str): path to the LAS file.

        Returns:
            Data: the point cloud formatted for later deep learning training.

        """
        las = laspy.read(las_filepath)
        pos = np.asarray([las.x, las.y, las.z], dtype=np.float32).transpose()

        x = np.asarray(
            [
                las[x_name]
                for x_name in [
                    "intensity",
                    "return_number",
                    "number_of_returns",
                    "red",
                    "green",
                    "blue",
                ]
            ],
            dtype=np.float32,
        ).transpose()

        return_number_idx = self.x_features_names.index("return_number")
        occluded_points = x[:, return_number_idx] > 1

        x[:, return_number_idx] = (x[:, return_number_idx]) / (
            self.return_number_normalization_max_value
        )
        num_return_idx = self.x_features_names.index("number_of_returns")
        x[:, num_return_idx] = (x[:, num_return_idx]) / (
            self.return_number_normalization_max_value
        )

        for idx, c in enumerate(self.x_features_names):
            if c in ["red", "green", "blue"]:
                assert x[:, idx].max() <= self.colors_normalization_max_value
                x[:, idx] = x[:, idx] / self.colors_normalization_max_value
                x[occluded_points, idx] = 0

        rgb_avg = (
            np.asarray(
                [las[x_name] for x_name in ["red", "green", "blue"]], dtype=np.float32
            )
            .transpose()
            .mean(axis=1, keepdims=True)
        )

        x = np.concatenate([x, rgb_avg], axis=1)

        try:
            # for LAS format V1.2
            y = las.classification.array.astype(int)
        except Exception:
            # for  LAS format V1.4
            y = las.classification.astype(int)

        return Data(
            pos=pos,
            x=x,
            y=y,
            las_filepath=las_filepath,
            x_features_names=self.x_features_names,
        )


def make_toy_dataset_from_test_file(
    prepared_data_dir: str,
    src_las_path: str = LAS_SUBSET_FOR_TOY_DATASET,
    split_csv: str = SPLIT_CSV_FOR_TOY_DATASET,
):
    """Prepare a toy dataset from a single, small LAS file.

    The file is first duplicated to get 2 LAS in each split (train/val/test),
    and then each file is splitted into .data files, resulting in a training-ready
    dataset loacted in td_prepared

    Args:
        src_las_path (str): input, small LAS file to generate toy dataset from
        split_csv (str): Path to csv with a `basename` (e.g. '123_456.las') and
        a `split` (train/val/test) columns specifying the dataset split.
        prepared_data_dir (str): where to copy files (`raw` subfolder) and to prepare
        dataset files (`prepared` subfolder)

    Returns:
        str: path to directory containing prepared dataset.

    """
    # Copy input file for full test isolation
    td_raw = osp.join(prepared_data_dir, "raw")
    td_prepared = osp.join(prepared_data_dir, "prepared")
    os.makedirs(td_raw)
    os.makedirs(td_prepared)
    # Make a "raw", unporcessed dataset with six files.
    basename = osp.basename(src_las_path)
    for s in ["train1", "train2", "val1", "val2", "test1", "test2"]:
        copy_path = osp.join(td_raw, basename.replace(".las", f".{s}.las"))
        shutil.copy(src_las_path, copy_path)

    # Prepare a Deep-Learning-ready dataset, using the split defined in the csv.
    data_prepper = FrenchLidarDataLogic(
        input_data_dir=td_raw,
        prepared_data_dir=td_prepared,
        split_csv=split_csv,
        input_tile_width_meters=110,
        subtile_width_meters=50,
    )
    data_prepper.prepare()
    return td_prepared


def _get_data_preparation_parser():
    """Gets a parser with parameters for dataset preparation.

    Returns:
        argparse.ArgumentParser: the parser.
    """
    parser = argparse.ArgumentParser(
        description="Prepare a Lidar dataset for deep learning.",
        epilog="If you need a toy dataset, you only need to"
        " set --origin=FR_TOY and specify a value for --prepared_data_dir.",
    )
    parser.add_argument(
        "--input_data_dir",
        type=str,
        default="./data/raw/",
        help="Path to directory with las files with unique basename and `.las` suffix.",
    )
    parser.add_argument(
        "--split_csv",
        type=str,
        default="./split.csv",
        help="Path to csv with a basename (e.g. '123_456.las') and split (train/val/test) columns specifying the dataset split.",
    )
    parser.add_argument(
        "--prepared_data_dir",
        type=str,
        default="./data/prepared/",
        help="Path to folder to save Data object train/val/test subfolders.",
    )
    parser.add_argument(
        "--origin", type=str, default="FR", choices=["FR", "CH", "FR_TOY"]
    )

    return parser


def main():
    """Code to execute to prepare a new set of Lidar tiles for training 3D segmentaiton neural networks.

    From a data directory containing point cloud in LAS format, and a scv specifying the dataset
    train/val/test split for each file (columns: split, basename, example: "val","123_456.las"),
    split the dataset, chunk the point cloud tiles into smaller subtiles, and prepare each sample
    as a pytorch geometric Data object.

    Echo numbers and colors are scaled to be in 0-1 range. Intensity and average color
    are not scaled b/c they are expected to be standardized later.

    To show help relative to data preparation, run:

        cd myria3d/datamodules/

        conda activate myria3d

        python loading.py -h


    """

    parser = _get_data_preparation_parser()
    args = parser.parse_args()
    if args.origin == "FR_TOY":
        make_toy_dataset_from_test_file(args.prepared_data_dir)
    elif args.origin == "FR":
        data_prepper = FrenchLidarDataLogic(**args.__dict__)
        data_prepper.prepare()
    elif args.origin == "CH":
        data_prepper = SwissTopoLidarDataLogic(**args.__dict__)
        data_prepper.prepare()
    else:
        raise KeyError(f"Data origin is invalid (currently: {args.origin})")


if __name__ == "__main__":
    main()

How to represent points clouds, specifying specific data signatures (e.g. choice of features).  Country-specific definitions inherit from abstract class LidarDataLogic, which implements general logics to split each point cloud into subtiles, format these subtiles, and save a learning-ready dataset splitted into train, val. A test set is also created, which is simply a copy of the selected test point clouds. This is so test phase conditions are similar to "in-the-wild" prediction conditions.  In particular, subclasses implement a "load_las" method that is used by the datamodule at test and inference time. # Files for the creation of a toy dataset. Abstract class to load, chunk, and save a point cloud dataset according to a train/val/test split.     load_las and its needed parameters ares specified in child classes. Load a point cloud in LAS format to memory and turn it into torch-geometric Data object.          Args:             las_filepath (str): path to the LAS file.          Returns:             Data: The point cloud formatted for later deep learning training. Prepare a dataset for model training and model evaluation.          Iterates through LAS files listed in a csv metadata file.         A `split` column specifies the train/val/test split of the dataset to be created. # Simply copy the LAS to the new test folder. # Load LAS into memory as a Data object with selected features, # then iteratively extract 50m*50m subtiles by filtering along x # then y axis. Serialize the resulting Data object using torch.save. Parse a LAS, extract and save each subtile as a Data object.          Args:             filepath (str): input LAS file             output_subdir_path (str): output directory to save splitted `.data` objects. # Ignore if empty # Ignore if empty Query files matching a basename in input_data_dir and its subdirectories.          Args:             input_data_dir (str): data directory          Returns:             [str]: first file path matching the query. Filter a data object on a chosen axis, using a relative position .         Modifies the original data object so that extracted future filters are faster.          Args:             data (Data): a pyg Data object with pos, x, and y attributes.             relative_pos (int): where the data to extract start on chosen axis (typically in range 0-1000)             axis (int, optional): 0 for x and 1 for y axis. Defaults to 0.          Returns:             Data: the data that is at most subtile_width_meters above relative_pos on the chosen axis. # select extract_by_axis applied on first axis x extract_by_axis applied on second axis y Save the subtile data object with torch.          Args:             subtile_data (Data): the object to save.             output_subdir_path (str): the subfolder to save it.             idx (int): an arbitrary but unique subtile identifier. # Exclusion of 65: artefacts, and 66: virtual points. Loads a point cloud in LAS format to memory and turns it into torch-geometric Data object.          Builds a composite (average) color channel on the fly.          Calculate NDVI on the fly.          x_features_names are {' - '.join(self.x_features_names)}          Args:             las_filepath (str): path to the LAS file.          Returns:             Data: the point cloud formatted for later deep learning training. # Filter out unused classes # Positions and base features # normalization # DEBUG: just to be sure that it is the same as before # Additional features : # Average color, that will be normalized on the fly based on single-sample # NDVI # for LAS format V1.2 # for  LAS format V1.4 Loads a point cloud in LAS format to memory and turns it into torch-geometric Data object.          Builds a composite (average) color channel on the fly.          Args:             las_filepath (str): path to the LAS file.          Returns:             Data: the point cloud formatted for later deep learning training. # for LAS format V1.2 # for  LAS format V1.4 Prepare a toy dataset from a single, small LAS file.      The file is first duplicated to get 2 LAS in each split (train/val/test),     and then each file is splitted into .data files, resulting in a training-ready     dataset loacted in td_prepared      Args:         src_las_path (str): input, small LAS file to generate toy dataset from         split_csv (str): Path to csv with a `basename` (e.g. '123_456.las') and         a `split` (train/val/test) columns specifying the dataset split.         prepared_data_dir (str): where to copy files (`raw` subfolder) and to prepare         dataset files (`prepared` subfolder)      Returns:         str: path to directory containing prepared dataset. # Copy input file for full test isolation # Make a "raw", unporcessed dataset with six files. # Prepare a Deep-Learning-ready dataset, using the split defined in the csv. Gets a parser with parameters for dataset preparation.      Returns:         argparse.ArgumentParser: the parser. Code to execute to prepare a new set of Lidar tiles for training 3D segmentaiton neural networks.      From a data directory containing point cloud in LAS format, and a scv specifying the dataset     train/val/test split for each file (columns: split, basename, example: "val","123_456.las"),     split the dataset, chunk the point cloud tiles into smaller subtiles, and prepare each sample     as a pytorch geometric Data object.      Echo numbers and colors are scaled to be in 0-1 range. Intensity and average color     are not scaled b/c they are expected to be standardized later.      To show help relative to data preparation, run:          cd myria3d/datamodules/          conda activate myria3d          python loading.py -h