## Atelier HPC SAGEO 2025

### Session 24/02/2025
#### Besoin
Cas d'application géomatique/raster

Chaine de traitement vue sur 3 outils différents
- OpenMOLE
- Jupyter avec Python
- Jupyter avec R
- {targets} + clustermq + SLURM
#### Hébergement

- IN2P3 : pas possible
- Relancer CRIIAN
- Onyxia
- https://datalab.sspcloud.fr/catalog/ide?search=para
- DASK
- [Argo workflows](https://argo-workflows.readthedocs.io/en/latest/)
-

#### Supports

Supports:
- Notebook Jupyter (Python/R)
- Git action de recherche MAGIS (OpenMole)

Humanum forge AR HPC :
https://gitlab.huma-num.fr/ar-magis-hpc

#### Presentation Chaine de traitement

Classification supervisé image sentinel 2, peupleraies

Image fourni par tuile (100 * 100km) => 90 tuiles
Modèle Random Forest échantillon entrainement classif supervisé => peuplier qui nous intéresse
A deployer sur les 90 tuiles, utilisation du Cluster CNES
Carte classification peuplier non peuplier, avec carte de confiance pour chaque pixel.

Pas parfait, parfois peuplier la ou il n'y en a pas. On avait pas de masque pour classer foret/nonforet donc il fallait faire du post traitement.
Trouver les règles expertes pour ne pas avoir besoin de ce masque hors forêt.
Classif france entière, 2017 - 2022, 90 tuiles par année, 90 cartes confiance pour ces 90tuiles => post traitement

Code dans un notebook Jupyter
- ? pour R / Rasterio python
- [rgdal2](https://github.com/thk686/rgdal2) pour R / equivalent python
- https://cran.r-project.org/web/packages/gdalraster/index.html
- https://cran.r-project.org/web/packages/vapour/index.html
- https://github.com/gearslaboratory/gdalUtils (obsolete ?)
- [stars](https://r-spatial.github.io/stars/articles/stars1.html) ? [terra] (https://rspatial.github.io/terra/index.html) ? pour R / Scipy pour python

##### 1ere etape

Binarisé les cartes de 1 à 8 (une par essence forestière), recup pixel avec peuplier

Intéressant ou pas côté pédagogie ?

##### 2eme etape

Conserver les valeurs de confiance pour les peuplier : lire l'image de classif et garder que la confiance associé au peuplier => valeur égale a 1

Traite par paire 1 carte classif, 1 carte confiance

##### 1 er post traitement : masque de pente

Seuil de pente, peuplier pousse pas au dela de 30%, donc on masque ce qui est au dessus.
MNT IGN => pente à l'échelle nationale
Pour chaque tuile, garde pixel si pente < 30%
Masque avec GDALCalque

##### 2 eme post traitement : masque de hauteur

Parfois confiance elevé, donc regarder aussi en fonction d'un seuil de haute.
Si <3m c'est pas de la forêt pour enlever les cultures
Carte nationale hauteur arbre, echelle france
Pour chaque tuile, correspondance avec <3m, et on filtre.

##### 3 eme post traitement : médian, lissage

Lissage pixel

##### 4 eme post traitement : masque de pente et hauteur

Si malgré tout çà on a une confiance < 50% c'est que c'est pas du peuplier
On a testé 30 et 50%

##### 5 eme post traitement : tamisage groupe de pixel

Tamisage pour > 0.5 hectares (50 pixels)

**Perspective :**

Comparer les classifications en multi-temporels, comparaison entre années pour chercher la cohérence entre année sur ce qui est ou pas peuplier.
Coupe c'est facile, mais quand c'est de l'aléatoire d'une année sur l'autre c'est problématique.

#### Scripts (durée, complexité)

Pur Python, pas d'optimisation, pas de paralléisation sur lecture gdal/raster
=> parallélisable :
- interne au niveau de la tuile
- entre tuile
- entre années

1 semaine de traitement pour france entière, Dell 4 coeurs, 8go ram

- 1 tuile : 2 Go
- 46Go Raster pour france entière pour 6 année ? (à vérifier)

#### Progression pédagogique

Matin :

- Introduction HPC => intro charabia HPC, mesocentre, etc.
- Introduction Chaine de Traitement : raster, etc.

- OpenMOLE :
- introduction openmole
Verrous : Quid des données un peu grosses => stockage in situ mésocentre, ou stockage LAL ?
- raster chaine traitement peuplier
- intro vecteur si il y a le temps => vector matching

Après midi :

**Introduction au raster et aux opérations **
- lecture, écriture avec GDAL
- manipulation simple de raster
=> application à la binarisation

**HPC**

Paralélisation entre noeuds, paralléisation sur plusieurs CPU d'un même noeud (gdal xarray numpy?)

*TODO LIST*

- Forge humanum :
- créer dépot Atelier
- créer site web associé

- A test (exemple minimum), deux notebooks différents (SEB pour le 28)
- quarto + jupyter python/R => un peu lourd (lance un serveur Jupyter)
- jupyter python/R => myst

- Relance mail + voir ensuite Jocelyn / Onyxia (Seb)
- Multilingue Python/R => regarder les équivalents Python/R (Nicolas)

- script OpenMOLE pour python : commencer (Juste)

- Preparation données, sous jeu de données uploadé => (Yousra) => Sharedocs ou autre (en fonction de la réponse humanum) ?

- Voir Dask et targets(R) => (tous/toustes en fonction des affinités et du temps dispo)

- Mail DSI Avignon avec file les infos des ip publiques ExModelo(SEB)

### Session 28/02/2025
- Zone étude : Nouvelle Aquitaine 27 tuiles
- tuiles de confiance : quelques Mo environ
- tuiles de classification : 500 Mo environ
- tuiles de pente MTN 10m: 1 Mo par tuile, France entière 22800 dalles (22,8 Go)
- Carte de hauteur : 22 Go France entière
- Problématique résolue pour le mésocentre calcul, çà sera le CRIANN à Rouen sur leur espace Jupyter , notre contact c'est Benoist Gaston : https://services.criann.fr/services/hpc/cluster-austral/guide/jupyter/

- Installation de R via conda/mamba
- https://www.happykhan.com/posts/conda-install-r/

- Conversion base de code vers R
- vapour est proche de GDAL mais ne fait que lire des fichiers donc petite déception
- Il est possible de faire des appels vers GDAL avec system() et link2GI1 semble pourvoir aider à s'assurer que le bouzin fonctionne
- terra est devenu l'outil de référence pour la manipulation de raster mais il faut retrouver les fonctions qui vont bien. De plus, avec [geotargets](https://geotargets.njtierney.com/), ça va bien dans le sens où on veut aller
- stars est orienté cube de données mais propose des fonctions intéressantes (ici une comparaison avec le paquet raster qui a été déprécié)
- [gdalraster](https://usdaforestservice.github.io/gdalraster/) binding directe à GDAL mais pas compatible avec geotargets encore (à creuser ? en profiter pour contribuer ?)
- Les deux méthodes qui m'inquiètent le plus sont celles qui proviennent de scipy.ndimages : `median_filter` et `generic_filter`. Pas sûr de pouvoir trouver un équivalent R (je n'ai pas approfondi plus que ça).
- `generic_filter` n'est pas appelé
- Appels Numpy à traduire
- `np.nan_to_num()`
- application de la valeur médiane aux pixels qui ont une valeur nulle
- Remplacement dans terra : https://rspatial.github.io/terra/reference/subst.html
- https://rspatial.github.io/terra/reference/focal.html
- `np.nanmedian()`
- calcul de la valeur médian des pixels valides (sans les NA)

- Test jupyter Lab / Jupyter Server avec Jupyterbook v2, çà marche pas trop mal, mais on n'échappe pas à l'execution d'un kernel jupyter dans tous les cas.

- Toutefois c'est moins lourd que ce ne que je pensais en terme d'install, pas plus que R en tout cas.

- Je vais tester ensuite on peut faire un notebook R aussi avec Myst / JupyterLab

### TODO :

- Forge avec fichier myst Python/R => Seb

- Fichier yaml conda (pour le 13 ):
- Windows et Linux => yousra (nico est sous windows)
- Environnement R => Nicolas, seb
- Environnement Python /=> Yousra, seb

- Penser la décomposition / optimisation
- Soit une grosse fonction avec toutes les étapes
- Fonction `smooth_all_rasters()` qui opère sur toutes les tuiles avec une boucle, paf on la met en fonction et on la daskifie/targetifie
- Test sur Workstation du Notebook => temps / durée du calcul en multi coeur pour info à Benoist Gaston => Sébastien, Yousra
- Soit plein de petite fonctions
- Création tout de suite de l'arborescence de fichiers pour optimiser

- Pédagogie :
- Penser à l'explication du passage de code sous forme de boucle à un code executable HPC / Dask et Targets
- Peut etre ajouter des exercices, en se basant sur la chaine existante : tamisage, etc. => créer une nouvelle targets, ou un nouveau truc dask

### Session 05/03/2025

- Yousra a déposé le jeu de données lundi sur Humanum, possibilité d'y accéder via webdav
- Yousra a commencé à reprendre le code pour diminuer les dépendances (Rasterio -> GDAL direct)
- regénération des données intermédiaires, blocage sur la pente

- COG proposé comme format alternatif pour stocker les tiff, avec du streaming,
- Voir [Optimiser vos rasters et générer des mosaïques au format COG (Cloud Optimized GeoTIFF) avec GDAL](https://geotribu.fr/articles/2025/2025-02-11_bonnes-pratiques-generation-raster-cog-avec-gdal/)
- Utilisation de COG avec OpenMole, reste la question de comment on charge les données ?

- Création des comptes CRIANN
- dépôt de dossier pour débloquer un compte pour Juste
- Demande de contact DSI en cours pour éviter les problèmes de connexion

- Jupyter sur le CRIANN
- 1 répertoire pour chaque langage
- 1 document identique
- à pousser sur la forge
- tutoriel à rédiger à propos de Myst

- Premiers essais avec targets/geotargets/terra (Nicolas)

### Session 13/03/2025

Dask :
- Optimisation inter images
- Optimisation intra image, découpage en chunk, mais pas forcément possible avec les images qui ont besoin de calculer des voisinnages ...

Ressources :

- Chaine => GDAL--> rasterio --> xarray / RioXArray --> Dask

- Rasterio XArray: https://docs.xarray.dev/en/stable/user-guide/io.html#rasterio & https://corteva.github.io/rioxarray/stable/readme.html
- Rasterio Rio XArray : https://corteva.github.io/rioxarray/stable/getting_started/getting_started.html
- GeoTIFF : https://knowledge.dea.ga.gov.au/notebooks/How_to_guides/Opening_GeoTIFFs_NetCDFs/
- XArray et Dask : https://docs.xarray.dev/en/stable/user-guide/dask.html#parallel-computing-with-dask
- Cog + Geotiff + Xarray : https://gallery.pangeo.io/repos/pangeo-data/landsat-8-tutorial-gallery/landsat8.html#Dask-Chunks-and-Cloud-Optimized-Geotiffs
- Fenetre mobile avec function XArray : https://tutorial.xarray.dev/intermediate/01-high-level-computation-patterns.html#sliding-windows-of-fixed-length-rolling et https://tutorial.xarray.dev/intermediate/01-high-level-computation-patterns.html#apply-an-existing-numpy-only-function-with-reduce

Visu : https://examples.holoviz.org/gallery/geospatial.html & https://tutorial.xarray.dev/intermediate/hvplot.html

#### Discussions Benoist Gaston (Criann)

- Jupyter Hub fonctionnel
- des demandes pour des kernel R donc plutôt motivé pour mettre en place

Pour le Dasbboard il y a un proxy à renseigner (clic sur bouton marche pas)

- Compte eleves générique déjà existant, on va passer par çà, soit avec des comptes liés à un projet spécifique, mais en tout cas il faut que çà soit homogene.

Soit projet associé à la formation, soit création d'un projet spécifique, mais il faudra la liste des eleves/apprenants => Volumétrie 20 ok, avec une première date de fin préinscription au 15 avril.

Interface StartMyServer => choix d'un profil de serveur (défini par Benoist en fonction des besoins)
Volumétrie à définir, réservation des temps de calcul peut se faire plus tardivement (pas la veille evidemment, 2 semaines au préalable ?)

Espace de stockage : Quota par défaut : 50 Go par utilisateur, quota par projet
si 100 Go : 2 années, sortie qq mégas par image

AUSTRAL: principalement job Python mais communauté IA demande Jupyter pour calibrer les modèles
Kernel R à positionner

Gestion de l'environnement conda et pip

Arborescence des fichiers : pas possible de le définir par avance mais possibilité d'enregister la dernière position

Dask peut lancer SLURM
Affichage du dashboard (Cluster Map et Task Stream)

Worker : nombre à définir 4 min -> 24 max
Jupyterlab dans des job slurm, avec var d'environnement affecté
Donc configuration d'un job/slurmcluster qui serait pas compatible, et il faudrait mettre en commentaire. => pas tout à fait le meme environnement qu'en local.

Point d'attention : serveur de 192 coeurs de calcul => plusieurs utilisateurs sur une machine, donc plusieurs ports ouverts par DASK pour afficher le dashboard

Session Jupyter exécutée sur un serveur, dans le cadre de la session, Job SLURM pas associé au job Jupyter mais variables d'environnement SLURM pour Jupyter peuvent être transmises.

Montrer passage à l'échelle entre LocalCluster (noeuds dédiés Jupyter, ressources moindres) et avec SLURMCluster.

On peut associer fichier de profil pour le job qui sera lancé par SLURM dans le jupyter. On peut associer une resa, certain nombre d'item déjà pris en compte. Mieux de faire une résa CPU, il faudra le prévoir un peu en amont de notre coté. C'est mieux que cela soit explicite plutot que d'avoir un profil préconfiguré.

Sur les noeud Austral on a REDHAT. Installation à la main sur de la RED HAT, dans des rep spécifiques. Outil module pour charger et décharger les variables d'environnements. Conda pour installer, mais module pour charger ces environnements là.

Pour les YML : export

### Session 20/03/2025

#### image Docker
- choix debian Bookworm/unstable:
- bookworm : Python 3.10 (ok)mais R 4.2 (en retard, pas tous les packages dispo)
- unstable : Python 3.13 (trop avancé) mais R 4.4
- décision partir de bookworm puis installer R 4.4

### Code R

`terra::focal()` devrait pouvoir remplacer le filtre median.

- Taille de la matrice 3x3 pixels
- Plusieurs options sont disponibles et doivent être évaluées:
- option par défaut `na.rm=FALSE`
`focal(r, 3, f)`
- option `na.rm=TRUE`
`focal(r, 3, f, na.rm=TRUE)`
- option `na.policy="only"` : change uniquement les cellules qui sont des NA
`focal(r, 3, f, na.policy="only", na.rm=TRUE)`
- option `na.policy="omit"` : ne change pas les cellules NA
`focal(r, 3, f, na.policy="omit", na.rm=TRUE)`

Les options `na.rm=FALSE` et `na.policy="omit"` peuvent être écartées car le code python remplace les NA par une valeur médiane (à l'échelle de l'image ? mais pas de la fenêtre).
Si `na.rm=TRUE`, alors le pixel considéré prend la valeur de la médiane des pixels adjacents hors NA.
Mais pour `na.policy="only"`, qu'est-ce qu'il se passe ?

Rappel que le but de l'atelier est de mettre en oeuvre la parallélisation si les résultats différentent en python et R ce n'est pas très important. La comparaison peut être utile si on va jusqu'à une publication.

#### Code de python

- problème compatibilité gdal/rasterio/xarray
- Yousra a trouvé une solution

- on créé deux md pour documenter les deux types de distributions

#### Distribution

- 1 année , distribution sur 27 tuiles
- point où appeler la parralélisation à identifier
- remplacer boucles for qui parcours les tuiles
- chaque tuile est traitée de manière indépendante, pas de gestion de la frontière
- Parraléliser dans le `main()` mais pas dans les fonctions
- Pour une tuile faire toute la chaine
- ou pour chaque fonction, passer toutes les tuiles
- temps de mise en route,
- 27 jobs vs (27 tuiles x 6 étapes) jobs
- gestion des entrées/sorties pour chaque étape
- si 1 job, les étapes intermédiaires peuvent être conservées en mémoire mais pas écrites dans un ficier
- dépends de comment sont configuré SLURM
- 30 s en amont et en sortie du job

### Session 26/03/2025

- Utiliser GDAL sans utiliser xarray
- essai Yousra: fonctionne pour les 3 premières étapes
- fonctionne sur les threads mais pas les process
- problème avec GDAL ? liaison GDAL/Dask pas documenté
- objet GDAL pas sériables => blocage au niveau du format
- Conversion rasterio => xarray faisable
- plutôt simple
- plus pythonique qu'avec GDAL
- env python (librairies geospatiales)
- proj=9.6.0=h0054346_0
- libgdal-core=3.10.2=hae73b24_2
- gdal=3.10.2=py310h79127d3_2
- geos=3.13.1=h97f6797_0
- geotiff=1.7.4=h239500f_2

https://rasterio.readthedocs.io/en/stable/topics/switch.html

https://pangeo-data.github.io/pangeo-openeo-BiDS-2023/part3/data_exploitability_pangeo.html

### Session 30/04/2025

- Mail aux participants (avec OpenSSH ? A voir avec la réponse), venir avec un Navigateur + Internet *Todo =>* Sébastien
- Verifier présence d'un poste de démo dans la salle sinon voir avec la DSI pour au moins un accès filaire + numéro téléphone DSI / Benoist "aucasou""
- Problème de seuil sur la fonction Python, mais Yousra va corriger
- Comparaison entre le python et le R en terme de résultats
- Définition du `pattern` dans `tar_target` avec un map
- [Exemple iceBreaker](https://github.com/njtierney/icebreaker/blob/5e330d6bc941e006880e2430d9c7d24016a21a76/_targets.R#L38C1-L42C5)
- Définition du pattern dans [targets](https://github.com/njtierney/icebreaker/blob/5e330d6bc941e006880e2430d9c7d24016a21a76/_targets.R#L38C1-L42C5)
- Programme Matin :
- Chapitre 0 : Présentation atelier (10 - 15 min)
- Motivations
- Objectifs
- Remerciements
- Chapitre 1 : Présentation de la chaine de traitement (45 min)
- Les données (et les formats) ()
- Pourquoi pas gdal directement ?
- La chaine de traitement (30 min)
- Diagramme d'exécution des traitements (extrait du rapport)
- Structure de données utilisés, equivalence (Python / R)
- Explication chaine blocs (10 min)
- Chaine de bloc initialement codé : la boucle
- Type de Parallélisation envisagés
- Niveau de l'année
- Niveau de la tuile => (tous les traitements sur 1 tuile) * nb tuiles
- Niveau de la tuile => (1 traitement pour l'ensemble des tuiles) * nb traitements
- Chapitre 2 : Chaine de traitement
- Présentation librairies utilisés
- Présentation Format Données (Xarray, SpatDataRaster, etc.)
- Execution de la chaine initiale (loop) sur une 1 tuile
- Calcul pour voir que c'est impossible
- Mise en pratique de la parallelisation de cette chaine
- Chapitre 2.1 Fonction 1
- Bloc Python / Bloc R / Comparaison
- Chapitre 2.3 Fonction 2
- Chapitre 3 : Parallélisation (début après midi)
- 3.1 avec R
- présentation targets => [reprendre élément doc Nicolas](https://gt-notebook.gitpages.huma-num.fr/workshop/je_notebook_2024/notebooks-r-plus-reproductibles/support/) (programmation fonctionnelle, Make-like)
- présentation géotargets => à faire
- présentation cluster / crew
- présentation SLURM vs LOCAL
- 3.2 avec Python
- présentation Dask
- présentation Xarray avec Dask => A faire
- présentation Future vs Map => Différentes façon de faire
- présentation SLURM vs LOCAL
- Exemples simple avec spatial
- Chapitre 4 : Perspectives / Conclusion
- Cache ? Rapidité ?
- support de cours
- quarto : possibilité de générer le doc sans kernel
- Myst : nécessite un kernel pour builder le document
- Dans un dépot séparé pour le moment
- => https://gitlab.huma-num.fr/ar-magis-hpc/workshops/atelier-hpc-notebooks/website-hpc-sageo-2025#
- Todo => CI Seb/Nicolas