====== IGN BdOrtho ======

Des versions de la BdOrtho sont accessibles en [[/glossaire/TMS|TMS]] via https://data.geopf.fr/tms/1.0.0/

Il y a aussi des version à télécharger sur <del>https://geoservices.ign.fr/telechargement-api/BDORTHO</del> le lien est maintenant https://cartes.gouv.fr/rechercher-une-donnee/dataset/IGNF_BD-ORTHO -> penser à scroller !!

Pour récupérer automatiquement ou détecter les nouveautés:
  * https://data.geopf.fr/telechargement/resource/BDORTHO?format=JP2-E080&page=2
  * https://data.geopf.fr/telechargement/resource/BDORTHO?zone=D036&format=JP2-E080&page=1
  * la réponse (Atom) contient les entrées pour la pagination: ''page'', ''pagesize'', ''pagecount'', ''totalentries'

===== Indre-et-Loire (37) =====

  * <del>https://geoservices.ign.fr/telechargement-api/BDORTHO?zone=D037</del>
  * 20 cm en RVB au 01/01/2025 BDORTHO_2-0_RVB-0M20_JP2-E080_LAMB93_D037_2025-01-01
  * 7 fichiers 7z contenant des images au format JP2/JPEG2000 pour 28.3 Go
  * Décompressés on obtient 302 images JP2 de 94 Mo chacune soit ~29 Go
    * 1 image découpée 640 pixels fait ~1520 fichiers pour ~70 Mo
    * au final 459 342 fichiers pour 24 Go

===== Outils =====

P'tit script Python pour découper les JP2 en JPG 640 pixels pour Yolo (//utilisé pour [[/informatique/ai_lm/ai_vision#panneaux_solaires_photovoltaiques|AI Vision panneaux photovoltaiques]]//) :

📢 À noter que c'est plusieurs heures de traitement, ~> 4h selon les départements et les dossiers des tuiles font entre 20 et 30 Go.

<code python>
"""
Installation:
pip3 install rasterio pillow numpy tqdm
"""
import rasterio
import os
import glob
import shutil
from rasterio.windows import Window
from PIL import Image
import numpy as np
from tqdm import tqdm
from concurrent.futures import ThreadPoolExecutor, as_completed
import argparse
import sys

#os.environ["GDAL_CACHEMAX"] = "24000000000"

def decouper_fichier_jp2(fichier, dossier_sortie, taille=640):
    nom_base = os.path.splitext(os.path.basename(fichier))[0]
    sous_dossier = os.path.join(dossier_sortie, nom_base)
    os.makedirs(sous_dossier, exist_ok=True)
    # Copie du fichier .tab
    fichier_tab = fichier.replace('.jp2', '.tab')
    if os.path.exists(fichier_tab):
        shutil.copy(fichier_tab, os.path.join(sous_dossier, f"{nom_base}.tab"))
    # Découpe le fichier JP2 en tuiles JPG
    # optimize only for write "blockxsize=taille*100, blockysize=taille*100"
    # GDAL possède un cache interne (GDAL_CACHEMAX) qui garde les blocs décompressés en mémoire.
    # Si un bloc est relu pour deux tuiles adjacentes, il est servi depuis le cache sans décompression supplémentaire.
    # Par défaut ce cache est à 5% de la RAM, mais on peut l'augmenter via :
    #os.environ["GDAL_CACHEMAX"] = "12000000000"
    #with rasterio.Env(CPL_DEBUG=True, GDAL_CACHEMAX=12000000000):
    with rasterio.open(fichier) as src:
        h, w = src.shape
        nb_tuiles_h = h // taille
        nb_tuiles_w = w // taille
        total_tuiles = nb_tuiles_h * nb_tuiles_w
        # Update progress bar every 100 tiles instead of every tile
        update_interval = 50
        tile_count = 0
        # réserve un buffer pour limiter une allocation à chaque tuile.
        buffer = np.empty((src.count, taille, taille), dtype=np.uint8)
        with tqdm(total=total_tuiles, desc=f"Découpe {nom_base}") as pbar:
            for y in range(0, h, taille):
                for x in range(0, w, taille):
                    window = Window(x, y, taille, taille)
                    if y + taille <= h and x + taille <= w:
                        tile = src.read(window=window, out=buffer)
                        tile_rgb = np.moveaxis(tile, 0, -1).astype(np.uint8)
                        img = Image.fromarray(tile_rgb)
                        img.save(f"{sous_dossier}/{nom_base}_x{x}_y{y}.jpg")
                        tile_count += 1
                        if tile_count % update_interval == 0:
                            pbar.update(update_interval)
            pbar.update(tile_count % update_interval)

def decouper_jp2_en_tuiles_parallele(dossier_entree, dossier_sortie, taille=640, nb_threads=4):
    os.makedirs(dossier_sortie, exist_ok=True)
    fichiers = glob.glob(os.path.join(dossier_entree, "*.jp2"))
    #fichiers = glob.glob(os.path.join(dossier_entree, "37-2025-0505-6715-LA93-0M20-E080.jp2"))
    with tqdm(total=len(fichiers), desc=f"Dossier {os.path.basename(dossier_entree)}", unit="image") as pbar:
        with ThreadPoolExecutor(max_workers=nb_threads) as executor:
            futures = [executor.submit(decouper_fichier_jp2, fichier, dossier_sortie, taille) for fichier in fichiers]
            for future in as_completed(futures):
                future.result()  # Attend la fin de chaque tâche
                pbar.update(1)

if __name__ == "__main__":

    parser = argparse.ArgumentParser(description="Découp les images BDOrtho IGN en tuiles de 640 pixels.")
    parser.add_argument("folder_in", help="Dossier des images à traiter (obligatoire)")
    parser.add_argument("folder_out", help="Dossier de sorties pour les tuiles découpées (obligatoire)")
    parser.add_argument("--threads", type=int, default=4, help="Combien de threads pour la découpe (défaut: 4).")
    args = parser.parse_args()

    if not os.path.isdir(args.folder_in):
        print(f"Error: Folder '{args.folder_in}' not found.")
        sys.exit(1)

    decouper_jp2_en_tuiles_parallele(
        #"/media/cyrille/Seagate Backup Plus Drive/BDORTHO_2-0_RVB-0M20_JP2-E080_LAMB93_D032_2025-01-01/ORTHOHR/1_DONNEES_LIVRAISON_2025-10-00030/OHR_RVB_0M20_JP2-E080_LAMB93_D32-2025/",
        #"/mnt/nas/IGN_BDOrtho/BDORTHO_2-0_RVB-0M20_JP2-E080_LAMB93_D005_2025-01-01/ORTHOHR/1_DONNEES_LIVRAISON_2026-01-00129/OHR_RVB_0M20_JP2-E080_LAMB93_D05-2025/",
        args.folder_in,
        #"/mnt/nas/IGN_BDOrtho/BDORTHO_tuiles-640_05/",
        args.folder_out,
        nb_threads=args.threads,
    )
</code>