Expérimentation du comptage de bouses dans un champ autour de Gramat.
Protocole en cours de rédaction
- Prise de vue DJI Mini 4 Pro sur RC2
- Orthomosaïque générée par RéalityCapture
- Script IA Claude pour Qgis
Pour les collègues relecteurs (biologie, aménagement) : ce document décrit une démarche testée les 21-23/10/2026 sur une prairie pâturée de 0,71 ha, du plan de vol jusqu'au comptage automatique des bouses. Il est soumis à votre avis avant que le protocole ne soit figé pour les BTSA GPN. Les points sur lesquels votre avis est attendu sont réunis à la section 6.
Convention : les cadres « Figure n – image à coller ici » indiquent l'emplacement de chaque copie d'écran.
1. Contexte et chaîne de travail
La mission drone se prépare entièrement dans QGIS : le script produit un fichier KMZ que la RC 2 exécute tel quel. Testée le 21/10/2026 (essai 2 points) puis le 22/10/2026 (mission complète, 286 photos), la chaîne fonctionne de bout en bout.
| Élément | Valeur |
|---|---|
| Objet du test | Décompte des bouses sur placettes, comparaison sol / drone |
| Site | Prairie pâturée, causse de Gramat (Lot), près de la route de Roussières |
| Zone de vol | zone_vol, 0,71 ha (≈ 167 × 56 m), EPSG:2154 |
| Drone | DJI Mini 4 Pro + radiocommande RC 2 (application DJI Fly) |
| Logiciels | QGIS 3.34 (console Python), RealityScan (orthophoto) |
| Dossier de travail | C:\Users\pierr\Downloads\Test Drone\ |
| Projet QGIS | Projet test drone.qgz |
La chaîne, dans l'ordre :
- QGIS : dessiner la zone de vol, placer le point de décollage, tirer les placettes.
- QGIS :
mission_drone_dji.pycalcule les waypoints, corrige les hauteurs par le relief (RGE ALTI de l'IGN) et écrit le KMZ à partir d'une mission « leurre » créée dans DJI Fly. - RC 2 : le KMZ remplace le fichier du leurre (même nom, même dossier) ; DJI Fly affiche la mission.
- Vol : décollage au point prévu, prise de vue automatique toutes les 2 s, lancement de la mission.
- QGIS :
trier_photos_vol.pyclasse les photos (orthophoto / placettes / écartées) d'après leur hauteur réelle au-dessus du sol. - RealityScan : orthophoto calculée avec les seules photos de la grille.
- À venir : comptage automatique des bouses sur toute l'orthophoto, puis contrôle de fiabilité (
compter_bouses_ortho.py).
Pourquoi passer par un leurre : DJI Fly n'a pas de bouton d'import pour le Mini 4 Pro. Le script recopie les réglages propres au drone contenus dans le leurre (modèle, mode de hauteur, action de fin) et ne remplace que les waypoints.
2. Phase 1 – Préparation et vol d'essai (21/10/2026)
L'essai à 2 points a validé l'import du KMZ, la prise de vue à intervalle, la nacelle à −90° et la correction des hauteurs. Il a aussi révélé la règle la plus importante : décoller exactement au point prévu.
Étape 1 – Dessiner la zone de vol et la passer en EPSG:2154
zone_volest dessinée sur l'orthophoto IGN, en retrait des houppiers.- Le projet était en EPSG:3857 (celui des flux IGN). En 3857, à 45° N, les « mètres » sont étirés de × 1,41 : un tampon de 10 m n'en ferait que 7.
- Conversion par clic droit > Exporter > Sauvegarder les entités sous… > SCR EPSG:2154, dans
couches\zone_vol.gpkg. Ne pas utiliser « Définir la projection », qui change l'étiquette sans recalculer les coordonnées. - Contrôle : emprise X ≈ 598 000, Y ≈ 6 414 100 (coordonnées Lambert 93 plausibles).
zone_vol dessinée, encore en EPSG:3857
Propriétés de la copie : EPSG:2154, emprise en Lambert 93
Étape 2 – Créer la couche du point de décollage
- Couche > Créer une couche > Nouvelle couche GeoPackage, dans le même fichier
zone_vol.gpkg, tabledecollage, type Point, EPSG:2154, aucun champ. - Piège évité : la première saisie mettait
zone_voldans « Nom de la table », ce qui aurait écrasé la zone. Toujours choisir « Ajouter une nouvelle couche », jamais « Écraser ».
Saisie corrigée : table decollage, sans champ
Étape 3 – Créer la mission leurre sur la RC 2
- Dans DJI Fly, créer une mission de 2 waypoints quelconques et l'enregistrer.
- Brancher la RC 2 au PC en USB.
- Copier le
.kmzle plus récent deAndroid/data/dji.go.v5/files/waypoint/<ID>/sur le PC, sans le modifier, et en garder une copie.
Le leurre contient droneEnumValue 68 (Mini 4 Pro), des hauteurs relatives au décollage et un retour automatique en fin de mission.
Étape 4 – Générer la mission d'essai
mission_drone_dji.py, mode Essai 2 points : 6 waypoints, 238 m, 1 min environ. La prairie est 6,5 m plus bas que le point de décollage prévu : A est à 23,5 m et B à 25,2 m sur la RC, soit 30 m au-dessus de l'herbe.
Message de fin du script, mission d'essai
Trajet d'essai : 0 et 5 superposés au-dessus du décollage, A et B au centre
Étape 5 – Charger la mission sur la RC 2
- Sauvegarder le KMZ d'origine du leurre.
- Copier le KMZ généré à sa place, même nom, même dossier.
- Fermer complètement DJI Fly, le rouvrir, ouvrir la mission : 6 points affichés, import validé.
Étape 6 – Vol d'essai et contrôle
Déroulé : décollage, stationnaire à 3 m, intervalle photo 2 s lancé à la main, puis lancement de la mission. 90 photos, analysées avec lire_hauteurs_photos_dji.py et RealityScan.
| Contrôle | Résultat |
|---|---|
| Intervalle 2 s maintenu pendant la mission | oui, aucune coupure |
| Montée verticale à 40 m | oui |
| Nacelle à −90° au waypoint 0 | oui |
| Hauteur en A (prévu 23,5 m) | 23,7 m |
| Longueur A → B (prévu 20 m) | 19,8 m |
| Retour et atterrissage automatiques | oui |
| Décollage au point prévu | non : ≈ 47 m plus au sud-ouest |
Le sol au décollage réel (390,0 m) est 2,3 m plus bas qu'au point prévu (392,3 m). Toutes les hauteurs étant comptées depuis le vrai décollage, le drone a volé à 27,9 m au-dessus de l'herbe au lieu de 30 m.
Règle retenue : décoller exactement au point decollage (à 2-3 m près), ou déplacer le point dans QGIS et régénérer la mission.
Vol d'essai dans RealityScan : trajectoire des caméras (descente, A → B, remontée)
3. Phase 2 – Placettes et mission complète (22/10/2026)
Huit placettes de 10 × 10 m ont été tirées (graine 2026), puis la mission complète générée : 54 waypoints, environ 1 450 m et 10 min, une seule batterie.
Étape 7 – Déplacer le point de décollage
Le point decollage est placé au lieu de décollage réel du 21/10, dont on sait qu'il convient : X = 598120,02 ; Y = 6414116,21 (coin sud-est de la zone, sol à 389,7 m). Le déplacement se fait par la console Python (script en annexe), pour avoir des coordonnées exactes. Contrôle : la couche contient un seul point.
Point de décollage déplacé au coin sud-est de zone_vol
Étape 8 – Tirer les placettes
tirage_placettes_generique.py avec les paramètres suivants :
| Paramètre | Valeur |
|---|---|
| Zone de tirage | zone_vol (0,71 ha), sans tampon |
| Placettes demandées / placées | 10 / 8 (P01 à P08, du nord au sud) |
| Forme, taille | carré de 10 × 10 m (100 m²) |
| Espacement entre bords | 10 m, soit 24,1 m exigés entre centres |
| Distance minimale obtenue entre centres | 24,3 m |
| Taux d'échantillonnage | 11,34 % |
| Graine | 2026 |
| Fichier retenu | placettes_20261022_0955.gpkg |
Un premier tirage (9 h 51) a été annulé : la valeur proposée de 4 m avait été gardée pour le côté du carré. Il a été relancé avec la même graine, et le fichier renommé placettes_20261022_0951_ANNULE.gpkg. Seules 8 placettes ont pu être placées à cause de l'espacement : le résultat est conservé, sans nouveau tirage.
Tirage annulé : carrés de 4 m
Tirage retenu : 8 carrés de 10 m et message de fin
Dossier Test Drone : tirage retenu et tirage annulé conservés
Étape 9 – Générer la mission complète
mission_drone_dji.py, mode Mission complète, avec zone_vol, decollage et centres. Résultat dans mission_20261022_1039 :
- 54 waypoints, environ 1 450 m, environ 10 min, environ 289 photos ;
- 6 lignes espacées de 10,2 m, les lignes extérieures sur les bords de la zone ;
- relief de −14,8 m (ouest) à +3,0 m (est) par rapport au décollage, d'où des hauteurs RC de 6,1 à 40 m pour un vol à 30 m (grille) et 15 m (placettes) au-dessus de l'herbe.
Mission complète dans QGIS : trajet, waypoints étiquetés, message de fin
Lecture des hauteurs : hauteur RC = hauteur au-dessus de l'herbe + (sol au waypoint − sol au décollage).
| Phase | WP | Rôle | Hauteur RC (m) | Au-dessus de l'herbe (m) | Sol / décollage (m) |
|---|---|---|---|---|---|
| Montée | 0 | Montée verticale au-dessus du décollage | 40,0 | 40,0 | 0,0 |
| 1 | Approche, 2 m avant la ligne 1 | 40,0 | 40,6 | −0,6 | |
| Grille | 2 → 7 | L1 est → ouest | 29,1 → 18,4 | 30 | −0,9 → −11,6 |
| 8 → 14 | L2 ouest → est | 16,8 → 29,8 | 30 | −13,2 → −0,2 | |
| 15 → 21 | L3 est → ouest | 30,5 → 15,2 | 30 | +0,5 → −14,8 | |
| 22 → 28 | L4 ouest → est | 15,3 → 31,5 | 30 | −14,7 → +1,5 | |
| 29 → 35 | L5 est → ouest | 32,5 → 15,8 | 30 | +2,5 → −14,2 | |
| 36 → 42 | L6 ouest → est | 16,6 → 33,0 | 30 | −13,4 → +3,0 | |
| Placettes | 43 | Approche de P01 | 32,6 | 30,0 | +2,6 |
| 44 | P01 | 17,6 | 15 | +2,6 | |
| 45 | P04 | 14,9 | 15 | −0,1 | |
| 46 | P07 | 12,9 | 15 | −2,1 | |
| 47 | P06 | 11,2 | 15 | −3,8 | |
| 48 | P03 | 10,0 | 15 | −5,0 | |
| 49 | P08 | 7,4 | 15 | −7,6 | |
| 50 | P05 | 6,1 | 15 | −8,9 | |
| 51 | P02 | 13,0 | 15 | −2,0 | |
| Retour | 52 | Remontée sur place | 40,0 | 41,7 | −1,7 |
| 53 | Au-dessus du décollage, puis retour automatique | 40,0 | 40,0 | 0,0 |
Le drone s'arrête en bout de ligne et passe sans s'arrêter sur les points intermédiaires. Il enchaîne les placettes à 15 m sans remonter, les trajets restant à plus de 5 m du bord de la zone.
Points de vigilance en vol :
- Les 6,1 m affichés sur P05 correspondent à 15 m au-dessus de l'herbe : le sol y est 8,9 m plus bas que le décollage.
- Les extrémités ouest (WP 7, 8, 21, 22, 35, 36) sont les points les plus loin (≈ 170 m) et les plus bas : garder le drone à vue.
- Les extrémités est débordent de 10 m vers la haie de la route : vérifier que les arbres ne dépassent pas ≈ 20 m.
- Entre deux points d'une ligne (30 m), le vol est en ligne droite : une bosse entre les deux n'est pas suivie.
Étape 10 – Charger la mission sur la RC 2
Même procédure qu'à l'étape 5 : KMZ copié à la place du leurre, DJI Fly redémarré, 54 points affichés.
4. Phase 3 – Vol complet, tri des photos, orthophoto (22/10/2026)
Le vol complet a produit 286 photos (environ 289 prévues). L'orthophoto est en cours de calcul dans RealityScan à partir de 193 photos de grille.
Étape 11 – Vol complet
- Vol vers 14 h 08 – 14 h 18, décollage au point
decollage. - 286 photos dans
DJI_001\22 oct\. - Les photos font environ 10 à 11 Mo : l'appareil était en 48 MP (12 MP prévus). L'image est 2 fois plus fine (environ 2,7 mm/pixel à 15 m, contre 5,4 mm en 12 MP), un avantage pour distinguer les classes de bouses, mais l'orthophoto est plus longue à calculer.
Dossier des photos du vol complet : 286 éléments
Étape 12 – Trier les photos
Avec 15 m de pente, la hauteur affichée sur la RC ne suffit plus à trier : une photo de la grille prise à l'ouest (≈ 15 m sur la RC) a la même hauteur qu'une photo de placette prise à l'est (≈ 17 m). trier_photos_vol.py calcule donc pour chaque photo sa hauteur réelle au-dessus du sol et tient compte de sa position :
| Classe | Conditions |
|---|---|
| ortho | nacelle verticale, ≈ 30 m au-dessus du sol (± 4 m), à moins de 15 m de zone_vol |
| placette | nacelle verticale, ≈ 15 m au-dessus du sol (± 4 m), à moins de 8 m d'un centre |
| écartée | décollage, montées et descentes, transits, atterrissage |
Les photos retenues sont copiées dans ortho\ et placettes\ (les originaux ne bougent pas). Les photos de placettes sont renommées, par exemple P05_MEILLEURE_DJI_….JPG.
Étape 13 – Orthophoto dans RealityScan
- Import du seul dossier
ortho\: 193 images. - Zone de reconstruction (« Verify Reconstruction Region ») : le cadre englobe les 6 lignes de vol, de la lisière ouest à la route de Roussières.
- Calcul : alignement, reconstruction, simplification et texture terminés ; projection orthographique en cours (≈ 23 min au total).
RealityScan : zone de reconstruction et positions des 193 photos de grille
RealityScan : calcul de l'orthophoto en cours
Étape 14 – Convertir l'orthophoto pour QGIS (23/10/2026)
L'export de RealityScan (OrthoRealityScan.tiff, 2 Go) s'ouvrait dans QGIS mais était très lent, et les pyramides étaient refusées. gdalinfo en a donné trois causes :
| Défaut de l'export | Conséquence |
|---|---|
| Coordonnées en WGS 84 (degrés), pas en Lambert 93 | reprojection à la volée de 537 millions de pixels |
| Image stockée tournée de 180° (origine au coin sud-est) | retournement à chaque lecture ; pyramides refusées |
| Bandes de 35 311 × 256 pixels, non tuilées ni compressées | lecture de lignes entières pour le moindre zoom |
Correction dans OSGeo4W Shell (installé avec QGIS), QGIS fermé :
cd /d "C:\Users\pierr\Downloads\Test Drone\couches"
gdalinfo OrthoRealityScan.tiff
gdalwarp -t_srs EPSG:2154 -tr 0.01 0.01 -tap -r average -multi -wo NUM_THREADS=ALL_CPUS -co TILED=YES -co COMPRESS=DEFLATE -co BIGTIFF=IF_SAFER OrthoRealityScan.tiff ortho_tmp.tif
gdal_translate -b 1 -b 2 -b 3 -mask 4 -co TILED=YES -co COMPRESS=JPEG -co PHOTOMETRIC=YCBCR -co JPEG_QUALITY=85 --config GDAL_TIFF_INTERNAL_MASK YES ortho_tmp.tif ortho_20261022.tif
gdaladdo -r average ortho_20261022.tif 2 4 8 16 32 64
gdalinfo ortho_20261022.tif
del ortho_tmp.tif
| Commande | Rôle |
|---|---|
gdalinfo |
diagnostic : système de coordonnées, taille de pixel, bandes, structure |
gdalwarp |
reprojette en Lambert 93, remet le nord en haut, rééchantillonne à 1 cm (moyenne), tuile le fichier |
gdal_translate |
compression JPEG qualité 85 ; la bande alpha devient un masque (bords transparents) |
gdaladdo |
aperçus réduits (1/2 à 1/64) pour un affichage immédiat |
Résultat, couches\ortho_20261022.tif : EPSG:2154, pixel de 1 cm, 22 273 × 9 894 pixels (222,7 × 98,9 m), tuiles de 256 × 256, 6 niveaux d'aperçus, 78 Mo au lieu de 2 Go (27 fois plus léger). Pourquoi 1 cm : une bouse d'environ 25 cm couvre encore environ 25 pixels de large, et le comptage fin se fait sur les photos de placettes (≈ 2,7 mm/pixel). Les exports d'origine vont dans couches\archives_ortho\.
Pour la prochaine fois : convertir systématiquement l'export RealityScan avec ces commandes, ou régler l'export en EPSG:2154 avec une taille de pixel de 0,01 m si l'option est trouvée dans RealityScan.
À faire
- [ ] Lire le bilan de
trier_photos_vol.py: hauteur de la grille ≈ 30 m partout, chaque placette avec ses photos, cap ≈ 0°, écart au centre < 2 m. - [ ] Exporter l'orthophoto (GeoTIFF, EPSG:2154) et l'ouvrir dans QGIS : fait le 23/10, voir l'étape 14.
- [ ] Écrire
compter_bouses_ortho.pyet le régler sur 4 placettes, puis le valider sur les 4 autres en le comparant au comptage au sol.
5. Phase 4 – Comptage automatique des bouses sur l'orthophoto (23/10/2026)
Un premier comptage sur toute la zone, sans placettes, donne 3 453 taches, soit environ 4 900 bouses/ha, couvrant 1,1 % de la surface. Ce résultat n'est pas encore validé : un contrôle visuel sur 5 carrés tirés au hasard doit mesurer sa fiabilité (étape 19).
Étape 15 – Principe de la méthode
Le script compter_bouses_ortho.py reprend l'idée d'une palette de couleurs et de pixels jointifs :
- Palette par l'exemple : des points posés sur l'orthophoto, sur des bouses (fraîche, croûtée, ancienne) et sur ce qui n'en est pas (herbe, herbe sèche, sol nu, autre). Le script lit leur couleur moyenne sur 5 × 5 pixels, en couleurs Lab (proches de la perception humaine).
- Classement : chaque pixel prend la classe de l'exemple le plus proche en couleur. Trop loin de tous les exemples, il n'est pas classé bouse.
- Ouverture puis fermeture : suppression des grains isolés, puis comblement des fentes et des trous d'une bouse.
- Nettoyage : suppression des taches plus petites que la surface minimale.
- Pixels jointifs : chaque groupe de pixels « bouse » qui se touchent, y compris en diagonale, devient une bouse.
- Filtres : surface minimale, compacité (une trace allongée est rejetée), part minimale de vrais pixels bouse ; une tache très grande est marquée « à vérifier ».
- Classe de chaque bouse : celle de la majorité de ses pixels.
Sorties : couche de points bouses_auto, contours bouses_contours, raster classes.tif (classe de chaque pixel), fichier parametres.txt (traçabilité).
Étape 16 – Constituer la palette d'exemples
Le vol complet a eu lieu par ciel couvert : il n'y a pas d'ombre portée, et une même bouse a la même couleur partout. C'est favorable au comptage, mais la palette n'est valable que pour ce vol.
| Classe | 1re version | Version utilisée (14 h 06) |
|---|---|---|
| bouse_fraiche | 9 | 12 |
| bouse_croutee | 5 | 6 |
| bouse_ancienne | 8 | 10 |
| herbe | 4 | 16 |
| herbe_seche | 4 | 18 |
| sol_nu | 8 | 8 |
| ombre | 1 | 1 |
| autre | 18 | 20 |
| Total | 57 | 91 |
La première version manquait d'exemples d'herbe : l'herbe couvre plus de 90 % de la prairie, et ses nuances doivent toutes être représentées pour ne pas être confondues avec des bouses.
Couche
echantillons_bouses sur l'orthophoto : vue d'ensemble de la répartition des 91 exemples, colorés par classe
Zoom au 1:50 : un exemple de chaque classe (bouse fraîche, croûtée, ancienne, herbe, herbe sèche, litière « autre »)
Étape 17 – Premier comptage (23/10, 14 h 06)
Réglages : surface minimale 150 cm², compacité minimale 0,25, écart de couleur maximal 25, fermeture 2 px, sans ouverture ni filtre de part minimale. Durée : 285 s.
| Résultat | Valeur |
|---|---|
| Taches retenues | 3 453, soit environ 4 900/ha |
| Dont bouse fraîche / croûtée / ancienne | 382 / 129 / 2 183 |
| Indéterminées (surtout des trous comblés) | 759 |
| À vérifier (> 2 500 cm²) | 1 |
| Écartées | 29 trop petites, 6 578 trop allongées |
| Surface médiane d'une tache | 225 cm² (80 % des taches entre 160 et 450 cm²) |
| Surface couverte | environ 79 m², soit 1,1 % de la zone |
Ce qui est plausible : 1,1 % de surface couverte et environ 290 journées de pâturage d'une vache (12 bouses par jour), sachant qu'une bouse ancienne reste visible plusieurs mois. Les filtres ont éliminé la pointe ouest, entièrement classée bouse au niveau des pixels.
Ce qui interroge : 759 taches indéterminées, une majorité de bouses anciennes (la classe que la couleur sépare le moins bien de la litière sèche) et une répartition très régulière, alors que les bouses se concentrent habituellement près des zones de repos, des haies et des passages.
Couche
bouses_auto sur l'orthophoto : vue d'ensemble de la zone de vol, points colorés par classe
Zoom au 1:50 sur une zone typique : orthophoto avec
bouses_contours, pour juger à l'œil des détections justes et faussesExtraits du raster classes.tif (6 × 6 m) avant filtrage. Vert : rien ; brun foncé : fraîche ; brun : croûtée ; beige : ancienne ; magenta : trous comblés par la fermeture.
classes.tif, centre de la prairie : grains beiges épars réunis par la fermeture
classes.tif, pointe ouest : végétation brun-beige presque entièrement classée bouse, puis écartée par les filtres
Étape 18 – Réglages corrigés (version 2 du script)
Ces réglages limitent l'agrégation de grains de litière en fausses bouses. La version 2 n'a pas encore tourné ; elle sera comparée à la première sur les mêmes carrés de contrôle.
| Réglage | Version 1 | Version 2 | Effet attendu |
|---|---|---|---|
| Ouverture (nouvelle) | — | 1 px | supprime les grains isolés avant le comblement |
| Fermeture | 2 px | 1 px | soude moins les grains voisins |
| Part minimale de vrais pixels (nouveau filtre) | — | 0,6 | rejette les taches faites surtout de trous comblés |
| Écart de couleur maximal | 25 | 15 | ne garde que les pixels très proches d'un exemple de bouse |
| Surface minimale | 150 cm² | 250 cm² | ≈ une bouse de 18 cm de diamètre |
| Compacité minimale | 0,25 | 0,35 | rejette les formes déchiquetées |
Une simulation sur classes.tif, sans le nouvel écart de couleur, donne environ 700 taches (environ 1 000/ha) au lieu d'environ 6 700 avant filtres de forme.
Même zoom au 1:50 que l'étape 17, avec les contours de la version 1 et de la version 2 superposés (deux couleurs)
Étape 19 – Contrôle de fiabilité (à faire)
- Tirer 5 carrés de 10 × 10 m dans
zone_volavectirage_placettes_generique.py, graine 2027. - Pour chaque carré, au 1:50 sur l'orthophoto, compter : les vrais positifs (VP, un point sur une vraie bouse), les faux positifs (FP, un point sans bouse dessous) et les manquées (FN, une bouse sans point).
- Calculer la précision = VP / (VP + FP) et le rappel = VP / (VP + FN), pour la version 1 puis pour la version 2.
| Carré | VP | FP | FN | Précision | Rappel |
|---|---|---|---|---|---|
| C01 | |||||
| C02 | |||||
| C03 | |||||
| C04 | |||||
| C05 | |||||
| Total |
Règle de décision proposée : une version est retenue si la précision et le rappel dépassent tous deux 80 %. La densité corrigée s'obtient alors par densité brute × précision / rappel.
FIGURE 23 – IMAGE À COLLER ICI
Les 5 carrés de contrôle sur l'orthophoto (vue d'ensemble, étiquettes C01 à C05)
Les 5 carrés de contrôle sur l'orthophoto (vue d'ensemble, étiquettes C01 à C05)
FIGURE 24 – IMAGE À COLLER ICI
Un carré de contrôle au 1:50 avec VP, FP et FN repérés (exemple de comptage)
Un carré de contrôle au 1:50 avec VP, FP et FN repérés (exemple de comptage)
6. Points soumis à l'avis des collègues
- Validation sans comptage au sol : le contrôle de l'étape 19 compare le script à un comptage à l'œil sur l'orthophoto, pas à un comptage sur le terrain. Est-ce suffisant pour un usage pédagogique, ou faut-il des placettes comptées au sol ?
- Seuil de 80 % de précision et de rappel : ce niveau vous paraît-il adapté ?
- Classes d'âge des bouses : à 1 cm/pixel, la bouse ancienne se distingue mal de la litière. Faut-il garder trois classes, ou regrouper en « fraîche et croûtée » et « ancienne (faible confiance) » ?
- Interprétation biologique : une densité de l'ordre de 1 000 à 5 000 bouses/ha est-elle cohérente avec le chargement et la durée de pâturage connus de cette parcelle ?
- Conditions de vol : imposer un ciel couvert, ou le soleil haut, pour que la palette reste valable ?
- Recul aux lisières : la zone de vol a été tracée en retrait des houppiers, sans tampon supplémentaire ; les bouses de lisière sont donc exclues. Est-ce acceptable pour la question posée ?
7. Les scripts et leurs paramètres expliqués
Chaque valeur proposée découle de trois contraintes : la taille de l'image au sol, la précision des capteurs (GPS, baromètre, MNT) et la sécurité autour des haies. Les tableaux ci-dessous donnent, pour chaque fenêtre de saisie, la valeur utilisée, la raison et l'effet d'une modification.
Repère : ce que couvre une photo du Mini 4 Pro
Objectif équivalent 24 mm, image 4:3, nacelle à −90° : l'emprise au sol vaut environ 1,44 × la hauteur dans le grand côté et 1,08 × la hauteur dans le petit côté.
| Hauteur au-dessus du sol | Emprise au sol | Taille d'un pixel en 12 MP | Taille d'un pixel en 48 MP |
|---|---|---|---|
| 15 m | ≈ 21,6 × 16,2 m | ≈ 5,4 mm | ≈ 2,7 mm |
| 30 m | ≈ 43 × 32 m | ≈ 1,1 cm | ≈ 5,4 mm |
Pendant les lignes, le nez suit la trajectoire : le petit côté (32 m à 30 m) est dans le sens de la marche, le grand côté (43 m) en travers.
mission_drone_dji.py – génération de la mission
| Fenêtre | Valeur | Pourquoi cette valeur | Si on la change |
|---|---|---|---|
| Type de mission | Essai 2 points, puis Mission complète | L'essai teste l'import et les réglages en 1 min, sans risque, avant de lancer 10 min de vol | — |
| Zone / Décollage / Placettes | zone_vol / decollage / centres |
Les noms par défaut des couches du projet. Avec « (aucune) » pour les placettes, on obtient l'orthophoto seule | Une autre couche de placettes = une autre série de photos verticales |
| Altitude du sol | Service IGN | RGE ALTI à 1 m, précis à ≈ 0,5 m, sans téléchargement. Une dalle MNT locale sert de secours sans internet | Sans correction du relief, le drone aurait volé à 45 m au-dessus du sol à l'ouest au lieu de 30 m |
| KMZ leurre | fichier copié de la RC 2 | Il contient les réglages propres au Mini 4 Pro, que DJI Fly exige ; le KMZ produit garde son nom | Un leurre d'un autre drone donne une mission refusée |
| Hauteur de sécurité | 40 m | Les arbres des haies font jusqu'à ≈ 20-25 m : 40 m laisse 15 m de marge pour la montée, les transits et le retour | Moins = risque de toucher la haie ; plus = montées plus longues, limite légale à 120 m |
| Hauteur de l'orthophoto | 30 m au-dessus du sol | Emprise de 43 × 32 m : la zone de 56 m de large se couvre en 6 lignes, avec des pixels de 0,5 à 1 cm | Plus bas = image plus fine mais plus de lignes et de photos ; plus haut = moins de détail sur les bouses |
| Vitesse sur les lignes | 3 m/s | Avec une photo toutes les 2 s, une photo tous les 6 m, soit (32 − 6) / 32 ≈ 80 % de recouvrement avant. Le flou de bougé reste ≈ 1 pixel (6 mm en 1/500 s) | Plus vite = moins de recouvrement, risque de trous dans l'orthophoto |
| Écart entre lignes | 12,5 m | (43 − 12,5) / 43 ≈ 70 % de recouvrement latéral, le minimum courant en photogrammétrie. Le script le réduit pour que les lignes extérieures tombent sur les bords (10,2 m ici) | Plus grand = moins de lignes mais une orthophoto fragile ; plus petit = plus de temps de vol |
| Débord hors zone | 10 m | Les bords d'une orthophoto sont mal reconstruits faute de recouvrement. Le débord ne s'applique qu'en bout de ligne, pas sur les côtés (les haies) | Plus = survol des haies ; 0 = bords flous |
| Points intermédiaires | tous les 30 m | Le drone suit le relief point par point. Sur une pente régulière (10 % ici), 30 m suffisent ; cela donne 54 waypoints au total | Plus serré = relief mieux suivi mais plus de waypoints ; 0 = extrémités seules, à éviter sur pente |
| Hauteur des placettes | 15 m au-dessus du sol | Emprise de 21,6 × 16,2 m : le carré de 10 m tient avec ≈ 3 m de marge de chaque côté, plus que la dérive GPS du drone (1 à 3 m) | 10 m : marge de 0,2 m, le carré sort de l'image ; 20 m : moins de détail |
| Stationnaire par placette | 6 s | 3 photos à 2 s ; on garde la plus nette et la plus centrée | Moins = risque d'une seule photo floue ; plus = vol plus long |
Réglages fixes dans le code : 2 m de décalage pour les montées et descentes (DJI Fly accepte mal deux waypoints superposés ; 2 m donne une descente presque verticale) et 5 m de marge au bord de la zone pour enchaîner les placettes sans remonter (on reste au-dessus de l'herbe, loin des branches).
tirage_placettes_generique.py – tirage des placettes
| Fenêtre | Valeur | Pourquoi cette valeur | Si on la change |
|---|---|---|---|
| Zone d'étude | zone_vol (prévu : zone_tiree) |
Le carré doit être entier dans la zone. Un tampon négatif éloigne des lisières (ombre, effet lisière) ; ici, la zone était déjà tracée en retrait | Sans retrait, une placette peut toucher la haie |
| Nombre de placettes | 10 (8 placées) | Une pré-étude demande 10 à 12 placettes pour estimer la variabilité, puis n = s² / (CV² × ȳ²) | Trop peu = variance mal estimée |
| Forme | Carré | Facile à repérer au pas et à la boussole sur le terrain, facile à découper sur la photo | Un cercle se trace au pas autour d'un point, mais se repère moins bien sur l'image |
| Côté du carré | 10 m (la valeur proposée est 4) | Les bouses sont rares et dispersées : avec 100 m², on évite trop de placettes à 0 bouse. Le carré de 4 m servait aux ligneux | 4 m = beaucoup de zéros, comptage peu informatif (erreur du tirage annulé) |
| Espacement entre bords | 10 m | Évite deux placettes sur la même photo et rend les comptages indépendants. Le script exige 24,1 m entre centres (2 demi-diagonales de 7,07 m + 10 m), pour que l'écart tienne même en diagonale | Plus grand = moins de placettes possibles (c'est ce qui a limité à 8) |
| Graine | 2026 | Même graine et mêmes paramètres = même tirage : n'importe qui peut le refaire et le vérifier | Changer la graine pour « un plus joli tirage » biaise l'échantillon : interdit |
trier_photos_vol.py – tri des photos
| Fenêtre | Valeur | Pourquoi cette valeur | Si on la change |
|---|---|---|---|
| Dossier des photos | DJI_001\22 oct |
Un dossier par vol : le tri et les copies restent rangés avec leurs originaux | — |
| Hauteur de l'orthophoto | 30 m | La hauteur visée par la mission | Doit être celle de la mission, sinon aucune photo n'est classée ortho |
| Hauteur des placettes | 15 m | Idem | Idem |
| Altitude du sol | Service IGN | Même source que pour la mission : les deux calculs sont cohérents | Une autre source peut décaler les hauteurs de quelques dizaines de cm |
| Copie | Oui | RealityScan reçoit le seul dossier ortho\, et les placettes sont nommées par numéro. Les originaux ne bougent pas |
Non = tri seulement dans le CSV et la couche QGIS |
Réglages fixes dans le code :
- Nacelle ≤ −85° : la nacelle annonce −90° à 1-2° près, et les photos de transit sont à 0° ou à −13/−40° (retour automatique).
- Tolérance de ± 4 m sur la hauteur : elle couvre le baromètre (± 0,5 m), le MNT (± 0,5 m) et le GPS sur la pente, tout en excluant les montées et descentes entre 15 et 30 m.
- 15 m autour de la zone pour l'orthophoto : les 10 m de débord plus la marge du GPS.
- 8 m autour d'un centre pour les placettes : on capte toutes les photos du stationnaire malgré la dérive, puis on retient comme « meilleure » la plus proche du centre.
compter_bouses_ortho.py – comptage automatique (version 2)
| Fenêtre | Valeur | Pourquoi cette valeur | Si on la change |
|---|---|---|---|
| Orthophoto / Zone | ortho_20261022 / zone_vol |
orthophoto convertie (Lambert 93, 1 cm, nord en haut) ; le script refuse une image tournée | — |
| Surface minimale | 250 cm² | une bouse fraîche fait 300 à 700 cm² ; en dessous, ce sont des brins, des cailloux ou de la litière | plus bas = plus de faux positifs ; plus haut = petites bouses anciennes perdues |
| Surface « à vérifier » | 2 500 cm² | au-delà de 50 × 50 cm, ce sont probablement plusieurs bouses collées | n'enlève rien, signale seulement |
| Compacité minimale | 0,35 | 4π × surface / périmètre² vaut environ 0,7 pour un disque dessiné en pixels, et moins de 0,1 pour une trace de roue | plus haut = bouses éclatées rejetées |
| Écart de couleur maximal | 15 (ΔE Lab) | 15 écarts « juste perceptibles » environ : seuls les pixels très proches d'un exemple sont gardés | plus haut = plus de litière prise pour des bouses |
| Ouverture | 1 px | efface les grains de moins de 3 pixels (3 cm) de large avant le comblement | 0 = grains conservés puis soudés en fausses bouses |
| Fermeture | 1 px | comble les fentes de 2 cm des croûtes | 2 px ou plus = grains voisins soudés |
| Part minimale de vrais pixels | 0,6 | au moins 60 % de la tache doit être de la couleur d'une bouse, pas du comblement | plus bas = taches indéterminées conservées |
Réglages fixes dans le code : fenêtre de 5 × 5 pixels autour de chaque exemple (5 cm, moyenne de la couleur sans mélanger avec l'herbe voisine), traitement par tuiles de 2 048 pixels (mémoire), pixels jointifs y compris en diagonale (connexité 8).
lire_hauteurs_photos_dji.py – lecture rapide des métadonnées
Il a servi au vol d'essai. Les réglages se font en tête de script : DOSSIER (le dossier des photos), NACELLE_MIN = -85, H_ORTHO = (23, 26) (hauteurs RC attendues entre A et B), H_PLACETTE = None. Il ne compare que la hauteur RC : il convient pour un terrain plat. Sur pente, utiliser trier_photos_vol.py.
Petits scripts de console
Ils sont donnés en annexe et n'ont pas de paramètre à saisir :
- déplacer le point de décollage aux coordonnées voulues ;
- calculer la hauteur réelle au-dessus de l'herbe après un décollage décalé ;
- exporter la table des waypoints en CSV.
8. Écarts au protocole et choix retenus
| Point | Prévu | Réalisé | Raison |
|---|---|---|---|
| Recul aux lisières | tampon de −10 m (règle de Gramat), puis −2 m | 0 m (tirage dans zone_vol) |
zone_vol déjà tracée en retrait des houppiers ; pas de 3ᵉ tirage |
| Nombre de placettes | 10 (pré-étude 10-12) | 8 | Espacement de 24,1 m entre centres |
| Tirage | un seul | 1 annulé (côté de 4 m), 1 retenu, même graine | Erreur de saisie, pas une recherche de répartition |
| Point de décollage | sur la route, au nord-est | coin sud-est, lieu du décollage réel du 21/10 | Hauteurs calculées depuis le vrai décollage |
| Hauteur des photos de placettes | 10 m | 15 m | Le carré de 10 m reste dans l'image malgré la dérive GPS |
| Déclenchement des photos | une action photo par waypoint | intervalle de 2 s lancé à la main | 54 waypoints au lieu de ≈ 300 ; fonctionnement vérifié à l'essai |
| Résolution | 12 MP | 48 MP | Réglage de l'appareil ; image plus fine, calcul plus long |
| Hauteurs | constantes sur la RC | corrigées par le relief (RGE ALTI) | DJI Fly compte depuis le décollage ; 18 m de dénivelé sur la zone |
| Survol des haies | — | WP0 et transits à 40 m ; lignes extérieures sur les bords ; débord en bout de ligne seulement | Sécurité |
| Retour automatique | — | ≥ 40 m dans DJI Fly | Haies |
9. Annexes – code des scripts
Voici les versions utilisées les 21 et 22/10/2026. Les fichiers de référence se trouvent aussi dans les documents du projet (claude/…). Chaque script s'exécute dans QGIS : Console Python > Afficher l'éditeur > Ouvrir > Exécuter.
A1. mission_drone_dji.py
Afficher le code (538 lignes)
# -*- coding: utf-8 -*-
"""
Mission drone DJI (Mini 4 Pro + RC 2 / DJI Fly) generee depuis QGIS
====================================================================
A executer dans QGIS : Extensions > Console Python > Afficher l'editeur >
ouvrir ce fichier > Executer.
Prerequis
- Projet QGIS ENREGISTRE, en EPSG:2154.
- Couche de polygones 'zone_vol' (zone a cartographier).
- Couche de points 'decollage' (UN point, sur la route).
- Facultatif : couche de points 'centres' (placettes tirees) -> photo
verticale en vol stationnaire au-dessus de chaque placette.
- Le .kmz d'une mission 'leurre' creee dans DJI Fly et copiee depuis la
RC 2 (Android/data/dji.go.v5/files/waypoint/<ID>/<ID>.kmz). Le script en
garde les reglages propres au drone et ne remplace que les waypoints.
Deroulement de la mission produite
WP0 a la verticale du decollage (montee droite au-dessus de la haie)
-> grille d'orthophoto (lignes paralleles au grand axe de la zone,
debordant de la zone, hauteur constante AU-DESSUS DU SOL)
-> placettes : descente, vol stationnaire nez au nord
-> remontee, retour au-dessus du decollage, retour automatique (RTH).
Nacelle a -90 deg pendant toute la mission. Les photos sont prises par
l'intervalle automatique (ex. 2 s) lance a la main AVANT la mission.
Les hauteurs sont corrigees par le relief : altitude du sol lue sur le
service d'altimetrie IGN (RGE ALTI, internet) ou sur une couche raster MNT
chargee dans le projet. Dans le KMZ, les hauteurs sont relatives au point
de decollage (c'est ce qu'attend DJI Fly).
Sorties (dossier du projet, sous-dossier mission_AAAAMMJJ_HHMM) :
- <meme nom que le leurre>.kmz : a copier sur la RC 2 a la place du leurre
- mission.gpkg : couches 'waypoints' et 'trajet' (controle visuel)
"""
import os
import re
import json
import math
import zipfile
import time
from datetime import datetime
from qgis.core import (QgsProject, QgsVectorLayer, QgsRasterLayer, QgsField, QgsFeature,
QgsGeometry, QgsPointXY, QgsWkbTypes, QgsUnitTypes,
QgsCoordinateReferenceSystem, QgsCoordinateTransform,
QgsVectorFileWriter, QgsBlockingNetworkRequest,
QgsPalLayerSettings, QgsVectorLayerSimpleLabeling)
from qgis.PyQt.QtCore import QVariant, QUrl
from qgis.PyQt.QtNetwork import QNetworkRequest
from qgis.PyQt.QtWidgets import QInputDialog, QMessageBox, QFileDialog
try:
PARENT = iface.mainWindow() # noqa: F821 (iface existe dans la console QGIS)
except NameError:
PARENT = None
AUCUNE = '(aucune)'
IGN = 'Service altimétrie IGN (internet)'
URL_IGN = ('https://data.geopf.fr/altimetrie/1.0/calcul/alti/rest/elevation.json'
'?lon=%s&lat=%s&resource=ign_rge_alti_wld&delimiter=|&zonly=true')
DECALAGE_MONTEE = 2.0 # m : decalage horizontal des points de montee / descente
MARGE_HAIE = 5.0 # m : vol bas entre placettes seulement a > 5 m du bord de la zone
class Annule(Exception):
pass
# ---------------------------------------------------------------- saisies
def demander_reel(titre, texte, defaut, mini, maxi, decimales=1):
v, ok = QInputDialog.getDouble(PARENT, titre, texte, defaut, mini, maxi, decimales)
if not ok:
raise Annule()
return v
def demander_choix(titre, texte, choix, defaut=0):
v, ok = QInputDialog.getItem(PARENT, titre, texte, choix, defaut, False)
if not ok:
raise Annule()
return v
def message(titre, texte, erreur=False):
if PARENT is not None:
(QMessageBox.critical if erreur else QMessageBox.information)(PARENT, titre, texte)
print(('ERREUR - ' if erreur else '') + titre + ' : ' + texte)
def choisir_couche(titre, texte, couches, nom_defaut, facultative=False):
noms = [l.name() for l in couches]
if facultative:
noms = [AUCUNE] + noms
defaut = noms.index(nom_defaut) if nom_defaut in noms else 0
nom = demander_choix(titre, texte, noms, defaut)
if nom == AUCUNE:
return None
return couches[[l.name() for l in couches].index(nom)]
# ---------------------------------------------------------------- altitudes du sol
def altitudes_ign(pts_wgs):
z = []
for i in range(0, len(pts_wgs), 50):
lot = pts_wgs[i:i + 50]
url = URL_IGN % ('|'.join('%.7f' % p.x() for p in lot), '|'.join('%.7f' % p.y() for p in lot))
req = QgsBlockingNetworkRequest()
if req.get(QNetworkRequest(QUrl(url))) != QgsBlockingNetworkRequest.NoError:
raise RuntimeError('service IGN injoignable : ' + req.errorMessage())
rep = json.loads(bytes(req.reply().content()).decode('utf-8'))
for v in rep['elevations']:
z.append(float(v['z'] if isinstance(v, dict) else v))
return z
def altitudes_raster(couche, pts, crs_pts, contexte):
tr = QgsCoordinateTransform(crs_pts, couche.crs(), contexte)
z = []
for p in pts:
v, ok = couche.dataProvider().sample(tr.transform(p), 1)
z.append(float(v) if ok else -99999.0)
return z
# ---------------------------------------------------------------- geometrie
def tourner(p, c, ang):
"""Rotation de p autour de c d'un angle ang (radians, sens trigonometrique)."""
dx, dy = p.x() - c.x(), p.y() - c.y()
ca, sa = math.cos(ang), math.sin(ang)
return QgsPointXY(c.x() + dx * ca - dy * sa, c.y() + dx * sa + dy * ca)
def tourner_geom(g, c, ang):
parties = g.asMultiPolygon() if g.isMultipart() else [g.asPolygon()]
return QgsGeometry.fromMultiPolygonXY([[[tourner(p, c, ang) for p in anneau] for anneau in poly]
for poly in parties])
def axe_principal(g):
"""Angle (radians) du plus grand cote du rectangle englobant oriente minimal."""
boite = g.orientedMinimumBoundingBox()[0]
pts = boite.asPolygon()[0]
best, ang = -1, 0.0
for a, b in zip(pts[:-1], pts[1:]):
d = math.hypot(b.x() - a.x(), b.y() - a.y())
if d > best:
best, ang = d, math.atan2(b.y() - a.y(), b.x() - a.x())
return ang
def grille(zone, debord, ecart, pas, depart):
"""Lignes de vol en serpentin. Renvoie une liste de lignes, chaque ligne = liste de points."""
zb = zone.buffer(debord, 16)
c = zb.centroid().asPoint()
ang = axe_principal(zone)
loc = tourner_geom(zb, c, -ang) # grand axe horizontal
bz = tourner_geom(zone, c, -ang).boundingBox()
bb = loc.boundingBox()
# lignes extremes SUR les bords de la zone (pas au-dessus des haies laterales),
# ecart reel <= ecart demande ; le debord s'applique en bout de ligne
n = max(2, int(math.ceil(bz.height() / ecart)) + 1)
pas_y = bz.height() / (n - 1)
lignes = []
for k in range(n):
y = bz.yMinimum() + 0.01 + k * (pas_y - 0.02 / (n - 1))
seg = QgsGeometry.fromPolylineXY([QgsPointXY(bb.xMinimum() - 1, y), QgsPointXY(bb.xMaximum() + 1, y)])
inter = seg.intersection(loc)
if inter.isEmpty():
continue
x1, x2 = inter.boundingBox().xMinimum(), inter.boundingBox().xMaximum()
nb = max(1, int(math.ceil((x2 - x1) / pas))) if pas > 0 else 1
lignes.append([QgsPointXY(x1 + (x2 - x1) * i / nb, y) for i in range(nb + 1)])
# serpentin
for k in range(1, len(lignes), 2):
lignes[k].reverse()
lignes = [[tourner(p, c, ang) for p in l] for l in lignes]
# 4 parcours possibles du meme serpentin : on commence au coin le plus proche du decollage
variantes = [lignes,
[list(reversed(l)) for l in lignes],
[list(reversed(l)) for l in reversed(lignes)],
[list(l) for l in reversed(lignes)]]
lignes = min(variantes, key=lambda v: v[0][0].distance(depart))
return lignes, math.degrees(ang), pas_y
def vers(a, b, d):
"""Point a d metres de a, en direction de b."""
L = a.distance(b)
if L < 1e-6:
return QgsPointXY(a.x(), a.y() - d)
return QgsPointXY(a.x() + (b.x() - a.x()) * d / L, a.y() + (b.y() - a.y()) * d / L)
# ---------------------------------------------------------------- WPML
def xml_waypoint(i, wp, lon, lat):
if wp['cap'] is None:
cap = ('<wpml:waypointHeadingMode>followWayline</wpml:waypointHeadingMode>'
'<wpml:waypointHeadingAngle>0</wpml:waypointHeadingAngle>')
else:
cap = ('<wpml:waypointHeadingMode>smoothTransition</wpml:waypointHeadingMode>'
'<wpml:waypointHeadingAngle>%d</wpml:waypointHeadingAngle>' % wp['cap'])
virage = ('toPointAndPassWithContinuityCurvature' if wp['passage']
else 'toPointAndStopWithDiscontinuityCurvature')
actions = ''
if i == 0:
actions += ('''
<wpml:actionGroup>
<wpml:actionGroupId>{g}</wpml:actionGroupId>
<wpml:actionGroupStartIndex>0</wpml:actionGroupStartIndex>
<wpml:actionGroupEndIndex>0</wpml:actionGroupEndIndex>
<wpml:actionGroupMode>parallel</wpml:actionGroupMode>
<wpml:actionTrigger><wpml:actionTriggerType>reachPoint</wpml:actionTriggerType></wpml:actionTrigger>
<wpml:action>
<wpml:actionId>0</wpml:actionId>
<wpml:actionActuatorFunc>gimbalRotate</wpml:actionActuatorFunc>
<wpml:actionActuatorFuncParam>
<wpml:gimbalHeadingYawBase>aircraft</wpml:gimbalHeadingYawBase>
<wpml:gimbalRotateMode>absoluteAngle</wpml:gimbalRotateMode>
<wpml:gimbalPitchRotateEnable>1</wpml:gimbalPitchRotateEnable>
<wpml:gimbalPitchRotateAngle>-90</wpml:gimbalPitchRotateAngle>
<wpml:gimbalRollRotateEnable>0</wpml:gimbalRollRotateEnable>
<wpml:gimbalRollRotateAngle>0</wpml:gimbalRollRotateAngle>
<wpml:gimbalYawRotateEnable>0</wpml:gimbalYawRotateEnable>
<wpml:gimbalYawRotateAngle>0</wpml:gimbalYawRotateAngle>
<wpml:gimbalRotateTimeEnable>0</wpml:gimbalRotateTimeEnable>
<wpml:gimbalRotateTime>0</wpml:gimbalRotateTime>
<wpml:payloadPositionIndex>0</wpml:payloadPositionIndex>
</wpml:actionActuatorFuncParam>
</wpml:action>
</wpml:actionGroup>''').format(g=wp['groupe'])
if wp['stationnaire'] > 0:
actions += ('''
<wpml:actionGroup>
<wpml:actionGroupId>{g}</wpml:actionGroupId>
<wpml:actionGroupStartIndex>{i}</wpml:actionGroupStartIndex>
<wpml:actionGroupEndIndex>{i}</wpml:actionGroupEndIndex>
<wpml:actionGroupMode>sequence</wpml:actionGroupMode>
<wpml:actionTrigger><wpml:actionTriggerType>reachPoint</wpml:actionTriggerType></wpml:actionTrigger>
<wpml:action>
<wpml:actionId>0</wpml:actionId>
<wpml:actionActuatorFunc>hover</wpml:actionActuatorFunc>
<wpml:actionActuatorFuncParam>
<wpml:hoverTime>{t}</wpml:hoverTime>
</wpml:actionActuatorFuncParam>
</wpml:action>
</wpml:actionGroup>''').format(g=wp['groupe'] + 1, i=i, t=int(round(wp['stationnaire'])))
return '''
<Placemark>
<Point><coordinates>{lon:.10f},{lat:.10f}</coordinates></Point>
<wpml:index>{i}</wpml:index>
<wpml:executeHeight>{h:.1f}</wpml:executeHeight>
<wpml:waypointSpeed>{v:.1f}</wpml:waypointSpeed>
<wpml:waypointHeadingParam>
{cap}
<wpml:waypointPoiPoint>0.000000,0.000000,0.000000</wpml:waypointPoiPoint>
<wpml:waypointHeadingAngleEnable>1</wpml:waypointHeadingAngleEnable>
<wpml:waypointHeadingPathMode>followBadArc</wpml:waypointHeadingPathMode>
<wpml:waypointHeadingPoiIndex>0</wpml:waypointHeadingPoiIndex>
</wpml:waypointHeadingParam>
<wpml:waypointTurnParam>
<wpml:waypointTurnMode>{virage}</wpml:waypointTurnMode>
<wpml:waypointTurnDampingDist>0</wpml:waypointTurnDampingDist>
</wpml:waypointTurnParam>
<wpml:useStraightLine>1</wpml:useStraightLine>{actions}
<wpml:waypointGimbalHeadingParam>
<wpml:waypointGimbalPitchAngle>-90</wpml:waypointGimbalPitchAngle>
<wpml:waypointGimbalYawAngle>0</wpml:waypointGimbalYawAngle>
</wpml:waypointGimbalHeadingParam>
</Placemark>'''.format(lon=lon, lat=lat, i=i, h=wp['h_rel'], v=wp['vitesse'], cap=cap,
virage=virage, actions=actions)
# ---------------------------------------------------------------- programme
def generer():
projet = QgsProject.instance()
dossier = projet.homePath()
if not dossier:
message('Projet non enregistré', "Enregistrez d'abord le projet, puis relancez.", True)
return
ctx = projet.transformContext()
vect = [l for l in projet.mapLayers().values() if isinstance(l, QgsVectorLayer)]
polys = [l for l in vect if l.geometryType() == QgsWkbTypes.PolygonGeometry]
points = [l for l in vect if l.geometryType() == QgsWkbTypes.PointGeometry]
rasters = [l for l in projet.mapLayers().values()
if isinstance(l, QgsRasterLayer) and l.providerType() == 'gdal' and l.bandCount() == 1]
# 1. choix
mode = demander_choix('Mission', 'Type de mission :',
['Mission complète (orthophoto + placettes)',
'Essai 2 points à 30 m (sans descente)'], 0)
essai = mode.startswith('Essai')
if not polys or not points:
message('Couches manquantes', "Il faut une couche de polygones 'zone_vol' et une couche "
"de points 'decollage'.", True)
return
zone_l = choisir_couche('Zone', 'Couche de la zone à cartographier :', polys, 'zone_vol')
dec_l = choisir_couche('Décollage', 'Couche du point de décollage :', points, 'decollage')
for l in (zone_l, dec_l):
if l.crs().isGeographic() or l.crs().mapUnits() != QgsUnitTypes.DistanceMeters:
message('Projection', "La couche '%s' doit être en EPSG:2154 (mètres)." % l.name(), True)
return
crs = zone_l.crs()
zone = QgsGeometry.unaryUnion([f.geometry() for f in zone_l.getFeatures()])
decs = [f.geometry() for f in dec_l.getFeatures()]
if len(decs) != 1:
message('Décollage', "La couche '%s' doit contenir UN seul point (elle en a %d)."
% (dec_l.name(), len(decs)), True)
return
g_dec = QgsGeometry(decs[0])
if dec_l.crs() != crs:
g_dec.transform(QgsCoordinateTransform(dec_l.crs(), crs, ctx))
p_dec = g_dec.asPoint()
pl_l = None
if not essai:
autres = [l for l in points if l is not dec_l]
if autres:
pl_l = choisir_couche('Placettes', 'Couche des centres de placettes\n'
'(« (aucune) » = orthophoto seule) :', autres, 'centres', True)
source_z = demander_choix('Altitude du sol', "Source de l'altitude du sol :",
[IGN] + [r.name() for r in rasters], 0)
leurre, _ = QFileDialog.getOpenFileName(PARENT, 'KMZ de la mission leurre (copié depuis la RC 2)',
dossier, 'Mission DJI (*.kmz)')
if not leurre:
raise Annule()
with zipfile.ZipFile(leurre) as z:
tpl = z.read('wpmz/template.kml').decode('utf-8')
wpl = z.read('wpmz/waylines.wpml').decode('utf-8')
if 'droneEnumValue' not in wpl or '<Folder>' not in wpl:
message('Leurre', "Ce fichier ne ressemble pas à une mission DJI Fly.", True)
return
# 2. parametres
h_wp0 = demander_reel('Hauteur de sécurité', 'Hauteur du WP0 au-dessus du décollage et des '
'transits (m) :\n(au-dessus des arbres de la haie)', 40, 15, 120)
h_ortho = demander_reel('Orthophoto', 'Hauteur de vol au-dessus du SOL (m) :', 30, 10, 120)
v_ortho = demander_reel('Orthophoto', 'Vitesse sur les lignes (m/s) :\n'
'(3 m/s + photo toutes les 2 s = une photo tous les 6 m)', 3.0, 1, 10)
if not essai:
ecart = demander_reel('Orthophoto', 'Écart entre lignes (m) :\n(12,5 m à 30 m ≈ 70 % de '
'recouvrement latéral)', 12.5, 3, 100)
debord = demander_reel('Orthophoto', 'Débord de la grille hors de la zone (m) :', 10, 0, 50)
pas = demander_reel('Relief', 'Distance entre points intermédiaires sur les lignes (m) :\n'
'(suivi du relief ; 0 = extrémités seulement)', 30, 0, 500)
if pl_l is not None:
h_pl = demander_reel('Placettes', 'Hauteur des photos de placettes au-dessus du SOL (m) :',
15, 5, 60)
t_pl = demander_reel('Placettes', 'Durée du vol stationnaire par placette (s) :\n'
'(6 s ≈ 3 photos à 2 s)', 6, 1, 60, 0)
# 3. trajet (x, y en EPSG du projet ; h = hauteur voulue au-dessus du sol, ou 'rel' fixe)
wps = []
def ajouter(p, h_sol=None, h_rel=None, type_='transit', vitesse=None, passage=False,
cap=None, stationnaire=0, etiquette=''):
wps.append({'p': QgsPointXY(p), 'h_sol': h_sol, 'h_rel_fixe': h_rel, 'type': type_,
'vitesse': vitesse or v_ortho, 'passage': passage, 'cap': cap,
'stationnaire': stationnaire, 'etiquette': etiquette})
ajouter(p_dec, h_rel=h_wp0, type_='WP0', etiquette='WP0')
orient = 0.0
if essai:
c = zone.centroid().asPoint()
a = QgsPointXY(c.x() - 10, c.y())
b = QgsPointXY(c.x() + 10, c.y())
ajouter(vers(a, p_dec, DECALAGE_MONTEE), h_rel=h_wp0, type_='transit')
ajouter(a, h_sol=h_ortho, type_='essai', etiquette='A')
ajouter(b, h_sol=h_ortho, type_='essai', etiquette='B')
else:
lignes, orient, ecart_reel = grille(zone, debord, ecart, pas, p_dec)
ajouter(vers(lignes[0][0], p_dec, DECALAGE_MONTEE), h_rel=h_wp0, type_='transit')
for k, l in enumerate(lignes, 1):
for j, p in enumerate(l):
ext = j in (0, len(l) - 1)
ajouter(p, h_sol=h_ortho, type_='ortho', passage=not ext,
etiquette=('L%d%s' % (k, 'a' if j == 0 else 'b')) if ext else '')
if pl_l is not None:
g_pl = []
tr = QgsCoordinateTransform(pl_l.crs(), crs, ctx) if pl_l.crs() != crs else None
champ = 'id_placette' if 'id_placette' in pl_l.fields().names() else None
for f in pl_l.getFeatures():
g = QgsGeometry(f.geometry())
if tr:
g.transform(tr)
g_pl.append((str(f[champ]) if champ else 'P%d' % f.id(), g.asPoint()))
interieur = zone.buffer(-MARGE_HAIE, 16)
precedent, bas = wps[-1]['p'], False
while g_pl: # plus proche voisin
k = min(range(len(g_pl)), key=lambda i: g_pl[i][1].distance(precedent))
nom, p = g_pl.pop(k)
seg = QgsGeometry.fromPolylineXY([precedent, p])
if not (bas and interieur.contains(seg)):
ajouter(vers(p, precedent, DECALAGE_MONTEE), h_sol=h_ortho, type_='transit')
ajouter(p, h_sol=h_pl, type_='placette', cap=0, stationnaire=t_pl, etiquette=nom)
precedent, bas = p, True
# remontee puis retour au-dessus du decollage
ajouter(vers(wps[-1]['p'], p_dec, DECALAGE_MONTEE), h_rel=h_wp0, type_='transit')
ajouter(p_dec, h_rel=h_wp0, type_='retour', etiquette='fin')
# 4. altitudes du sol
wgs = QgsCoordinateReferenceSystem('EPSG:4326')
vers_wgs = QgsCoordinateTransform(crs, wgs, ctx)
pts = [p_dec] + [w['p'] for w in wps]
try:
if source_z == IGN:
z = altitudes_ign([vers_wgs.transform(p) for p in pts])
else:
z = altitudes_raster(rasters[[r.name() for r in rasters].index(source_z)], pts, crs, ctx)
except Exception as e:
message('Altitude du sol', "Lecture impossible (%s).\nCharger une dalle MNT (.tif) dans le "
"projet et la choisir comme source." % e, True)
return
if min(z) < -1000:
message('Altitude du sol', "Pas d'altitude pour certains points (hors de la dalle MNT ?).", True)
return
z0 = z[0]
for w, zz in zip(wps, z[1:]):
w['z_sol'] = zz
if w['h_rel_fixe'] is not None:
# transit : jamais plus bas que h_ortho au-dessus du sol local
w['h_rel'] = max(w['h_rel_fixe'], zz - z0 + h_ortho)
else:
w['h_rel'] = w['h_sol'] + zz - z0
w['h_sur_sol'] = w['h_rel'] - (zz - z0)
# 5. bilan
dist, duree = 0.0, 0.0
for a, b in zip(wps[:-1], wps[1:]):
d = math.sqrt(a['p'].sqrDist(b['p']) + (a['h_rel'] - b['h_rel']) ** 2)
dist += d
duree += d / b['vitesse'] + b['stationnaire'] + (0 if b['passage'] else 2)
# 6. ecriture du KMZ
g = 1
for w in wps:
w['groupe'] = g
g += 2
lonlat = [vers_wgs.transform(w['p']) for w in wps]
corps = ''.join(xml_waypoint(i, w, ll.x(), ll.y()) for i, (w, ll) in enumerate(zip(wps, lonlat)))
entete = re.search(r'^(.*?<Folder>)', wpl, re.S).group(1)
entete = re.sub(r'<wpml:globalTransitionalSpeed>[^<]*<',
'<wpml:globalTransitionalSpeed>%.1f<' % v_ortho, entete)
dossier_folder = re.search(r'<Folder>(.*?)<Placemark>', wpl, re.S).group(1)
dossier_folder = re.sub(r'<wpml:distance>[^<]*<', '<wpml:distance>%d<' % dist, dossier_folder)
dossier_folder = re.sub(r'<wpml:duration>[^<]*<', '<wpml:duration>%d<' % duree, dossier_folder)
dossier_folder = re.sub(r'<wpml:autoFlightSpeed>[^<]*<',
'<wpml:autoFlightSpeed>%.1f<' % v_ortho, dossier_folder)
if 'relativeToStartPoint' not in dossier_folder:
message('Leurre', "Le leurre n'est pas en hauteurs relatives au décollage : recréez-le "
"dans DJI Fly sans modifier le réglage d'altitude.", True)
return
waylines = entete[:-len('<Folder>')] + '<Folder>' + dossier_folder + corps + \
'\n </Folder>\n </Document>\n</kml>\n'
maintenant = str(int(time.time() * 1000))
template = re.sub(r'<wpml:updateTime>[^<]*<', '<wpml:updateTime>%s<' % maintenant, tpl)
sortie = os.path.join(dossier, 'mission_%s' % datetime.now().strftime('%Y%m%d_%H%M'))
os.makedirs(sortie, exist_ok=True)
kmz = os.path.join(sortie, os.path.basename(leurre))
with zipfile.ZipFile(kmz, 'w', zipfile.ZIP_DEFLATED) as zf:
zf.writestr('wpmz/template.kml', template)
zf.writestr('wpmz/waylines.wpml', waylines)
# 7. couches de controle
champs = [QgsField('index', QVariant.Int), QgsField('type', QVariant.String),
QgsField('etiquette', QVariant.String), QgsField('h_rel_decollage', QVariant.Double),
QgsField('h_sur_sol', QVariant.Double), QgsField('z_sol', QVariant.Double),
QgsField('vitesse', QVariant.Double), QgsField('stationnaire_s', QVariant.Double)]
cw = QgsVectorLayer('Point?crs=%s' % crs.authid(), 'waypoints', 'memory')
ct = QgsVectorLayer('LineString?crs=%s' % crs.authid(), 'trajet', 'memory')
cw.dataProvider().addAttributes(champs)
cw.updateFields()
feats = []
for i, w in enumerate(wps):
f = QgsFeature(cw.fields())
f.setAttributes([i, w['type'], w['etiquette'], round(w['h_rel'], 1), round(w['h_sur_sol'], 1),
round(w['z_sol'], 2), w['vitesse'], w['stationnaire']])
f.setGeometry(QgsGeometry.fromPointXY(w['p']))
feats.append(f)
cw.dataProvider().addFeatures(feats)
ft = QgsFeature()
ft.setGeometry(QgsGeometry.fromPolylineXY([w['p'] for w in wps]))
ct.dataProvider().addFeatures([ft])
gpkg = os.path.join(sortie, 'mission.gpkg')
premier = True
for lyr in (ct, cw):
opt = QgsVectorFileWriter.SaveVectorOptions()
opt.driverName = 'GPKG'
opt.layerName = lyr.name()
opt.actionOnExistingFile = (QgsVectorFileWriter.CreateOrOverwriteFile if premier
else QgsVectorFileWriter.CreateOrOverwriteLayer)
QgsVectorFileWriter.writeAsVectorFormatV3(lyr, gpkg, ctx, opt)
premier = False
for nom in ('trajet', 'waypoints'):
couche = QgsVectorLayer('%s|layername=%s' % (gpkg, nom), 'mission_' + nom, 'ogr')
if nom == 'waypoints':
r = QgsPalLayerSettings()
r.fieldName = "\"index\" || ' (' || \"h_rel_decollage\" || ' m)'"
r.isExpression = True
r.enabled = True
couche.setLabeling(QgsVectorLayerSimpleLabeling(r))
couche.setLabelsEnabled(True)
projet.addMapLayer(couche)
# 8. message final
hs = [w['h_rel'] for w in wps]
alertes = []
if len(wps) > 99:
alertes.append('%d waypoints : au-delà de 99, vérifier que DJI Fly accepte la mission '
'(sinon augmenter la distance entre points intermédiaires).' % len(wps))
if min(z[1:]) - z0 < -10 or max(z[1:]) - z0 > 10:
alertes.append('Relief : le sol varie de %.0f à %+.0f m par rapport au décollage.'
% (min(z[1:]) - z0, max(z[1:]) - z0))
if duree > 15 * 60:
alertes.append('Durée estimée > 15 min : prévoir deux batteries ou réduire la zone.')
lignes_msg = [
'Mission %s : %d waypoints.' % ('d\'essai' if essai else 'complète', len(wps)),
'Longueur ≈ %.0f m ; durée ≈ %.0f min ; ≈ %d photos à 2 s.' % (dist, duree / 60, duree / 2),
'Hauteurs relatives au décollage : %.1f à %.1f m.' % (min(hs), max(hs)),
'Altitude du sol au décollage : %.1f m (%s).' % (z0, source_z)]
if not essai:
lignes_msg.append('%d lignes parallèles au grand axe de la zone, écart réel %.1f m.'
% (len(lignes), ecart_reel))
lignes_msg += alertes
lignes_msg += ['', 'KMZ : %s' % kmz,
'Copier ce fichier sur la RC 2 à la place du leurre, MÊME NOM, même dossier.',
"Vérifier le trajet et les hauteurs (couche mission_waypoints) avant le vol."]
message('Mission générée', '\n'.join(lignes_msg))
try:
generer()
except Annule:
print('Génération annulée.')
A2. tirage_placettes_generique.py
Afficher le code (288 lignes)
# -*- coding: utf-8 -*-
"""
Tirage aleatoire de placettes dans une zone d'etude - script generique
=======================================================================
A executer dans QGIS : Extensions > Console Python > Afficher l'editeur >
ouvrir ce fichier > Executer.
Prerequis
- Projet QGIS ENREGISTRE (les resultats sont ecrits dans son dossier).
- Une couche de polygones dessinee par l'utilisateur, appelee par defaut
'zone_tiree', dans un SCR en metres (ex. EPSG:2154 Lambert 93).
Un polygone = une strate. Si la couche a un champ 'strate', son contenu
sert de nom de strate ; sinon les strates sont numerotees S1, S2...
Le script demande, par des fenetres de saisie :
- la couche zone ;
- le nombre de placettes ;
- la forme (carre ou cercle) et la dimension (cote ou rayon, en m) ;
- la distance minimale entre les BORDS de deux placettes (m) ;
- la repartition entre strates (si plusieurs polygones) ;
- la graine aleatoire (meme graine = meme tirage).
Chaque placette est tiree ENTIEREMENT a l'interieur de la zone (pas seulement
son centre). Resultat : placettes_AAAAMMJJ_HHMM.gpkg avec deux couches,
'centres' (points) et 'placettes' (carres ou cercles), numerotees P01, P02...
du nord au sud.
"""
import os
import math
import random
from datetime import datetime
from qgis.core import (QgsProject, QgsVectorLayer, QgsField, QgsFeature, QgsGeometry,
QgsPointXY, QgsWkbTypes, QgsUnitTypes, QgsVectorFileWriter,
QgsSpatialIndex, QgsPalLayerSettings, QgsVectorLayerSimpleLabeling)
from qgis.PyQt.QtCore import QVariant
from qgis.PyQt.QtWidgets import QInputDialog, QMessageBox
try:
PARENT = iface.mainWindow() # noqa: F821 (iface existe dans la console QGIS)
except NameError:
PARENT = None
NOM_ZONE_DEFAUT = 'zone_tiree'
ESSAIS_PAR_PLACETTE = 2000
class Annule(Exception):
pass
# ---------------------------------------------------------------- saisies
def demander_entier(titre, texte, defaut, mini, maxi):
v, ok = QInputDialog.getInt(PARENT, titre, texte, defaut, mini, maxi, 1)
if not ok:
raise Annule()
return v
def demander_reel(titre, texte, defaut, mini, maxi, decimales=1):
v, ok = QInputDialog.getDouble(PARENT, titre, texte, defaut, mini, maxi, decimales)
if not ok:
raise Annule()
return v
def demander_choix(titre, texte, choix, defaut=0):
v, ok = QInputDialog.getItem(PARENT, titre, texte, choix, defaut, False)
if not ok:
raise Annule()
return v
def message(titre, texte, erreur=False):
if PARENT is not None:
(QMessageBox.critical if erreur else QMessageBox.information)(PARENT, titre, texte)
print(('ERREUR - ' if erreur else '') + titre + ' : ' + texte)
# ---------------------------------------------------------------- geometrie
def forme_placette(x, y, forme, dim):
"""Carre de cote dim, ou cercle de rayon dim, centre en (x, y)."""
if forme == 'Carré':
h = dim / 2.0
return QgsGeometry.fromPolygonXY([[QgsPointXY(x - h, y - h), QgsPointXY(x + h, y - h),
QgsPointXY(x + h, y + h), QgsPointXY(x - h, y + h),
QgsPointXY(x - h, y - h)]])
return QgsGeometry.fromPointXY(QgsPointXY(x, y)).buffer(dim, 36)
def repartir(effectif, surfaces, mode):
"""Nombre de placettes par strate (methode des plus forts restes)."""
k = len(surfaces)
if k == 1:
return [effectif]
if mode.startswith('Égale'):
base = [effectif // k] * k
for i in range(effectif - sum(base)):
base[i] += 1
return base
total = sum(surfaces)
brut = [effectif * s / total for s in surfaces]
base = [int(math.floor(b)) for b in brut]
restes = sorted(range(k), key=lambda i: brut[i] - base[i], reverse=True)
for i in restes[:effectif - sum(base)]:
base[i] += 1
return base
# ---------------------------------------------------------------- programme
def tirer():
projet = QgsProject.instance()
dossier = projet.homePath()
if not dossier:
message('Projet non enregistré', "Enregistre d'abord le projet (Projet > Enregistrer), puis relance.", True)
return
# 1. couche zone
polys = [l for l in projet.mapLayers().values()
if isinstance(l, QgsVectorLayer) and l.geometryType() == QgsWkbTypes.PolygonGeometry]
if not polys:
message('Aucune zone', "Dessine d'abord une couche de polygones 'zone_tiree'.", True)
return
noms = [l.name() for l in polys]
defaut = noms.index(NOM_ZONE_DEFAUT) if NOM_ZONE_DEFAUT in noms else 0
zone = polys[noms.index(demander_choix('Zone d\'étude', 'Couche de polygones où tirer les placettes :',
noms, defaut))]
if zone.crs().isGeographic() or zone.crs().mapUnits() != QgsUnitTypes.DistanceMeters:
message('Projection', "La couche '%s' doit être dans un SCR en mètres (ex. EPSG:2154), "
"pas en %s. La reprojeter d'abord." % (zone.name(), zone.crs().authid()), True)
return
strates = []
champ_strate = 'strate' if 'strate' in zone.fields().names() else None
for i, f in enumerate(zone.getFeatures(), 1):
g = f.geometry()
if g is None or g.isEmpty():
continue
nom = str(f[champ_strate]) if champ_strate and f[champ_strate] not in (None, '') else 'S%d' % i
strates.append({'nom': nom, 'geom': g, 'surface': g.area()})
if not strates:
message('Zone vide', "La couche '%s' ne contient aucun polygone." % zone.name(), True)
return
surface_totale = sum(s['surface'] for s in strates)
# 2. parametres
n = demander_entier('Nombre de placettes', 'Combien de placettes tirer au hasard ?\n'
'(zone : %.2f ha, %d polygone(s))' % (surface_totale / 1e4, len(strates)), 12, 1, 500)
forme = demander_choix('Forme', 'Forme des placettes :', ['Carré', 'Cercle'], 0)
if forme == 'Carré':
dim = demander_reel('Dimension', 'Côté du carré (m) :', 4.0, 0.5, 500.0)
demi_diag = dim * math.sqrt(2) / 2
surface_placette = dim * dim
else:
dim = demander_reel('Dimension', 'Rayon du cercle (m) :', 2.26, 0.2, 250.0, 2)
demi_diag = dim
surface_placette = math.pi * dim * dim
ecart = demander_reel('Espacement', 'Distance minimale entre les BORDS de deux placettes (m) :\n'
'(conseil : au moins la taille d\'une placette)',
round(max(dim, 5.0), 1), 0.0, 1000.0)
mode = 'Unique'
if len(strates) > 1:
mode = demander_choix('Répartition', '%d polygones (strates). Répartir les placettes :' % len(strates),
['Proportionnelle à la surface', 'Égale entre strates'], 0)
graine = demander_entier('Graine aléatoire', 'Graine (même graine = même tirage) :',
datetime.now().year, 0, 999999)
# distance minimale entre centres garantissant l'ecart entre bords (cas le plus defavorable)
dist_centres = 2 * demi_diag + ecart
taux = 100.0 * n * surface_placette / surface_totale
if taux > 30:
message('Attention', 'Les placettes couvriraient %.0f %% de la zone : tirage probablement '
'impossible. Réduis le nombre, la taille ou l\'espacement.' % taux)
# 3. tirage
rng = random.Random(graine)
effectifs = repartir(n, [s['surface'] for s in strates], mode)
tires = [] # (x, y, nom_strate)
index = QgsSpatialIndex()
bilan = []
for s, k in zip(strates, effectifs):
g = s['geom']
# zone ou un CENTRE peut tomber pour que la placette reste entiere dans le polygone
coeur = g.buffer(-demi_diag, 16) if forme == 'Carré' else g.buffer(-dim, 16)
if coeur.isEmpty():
bilan.append('%s : polygone trop étroit pour une placette de cette taille' % s['nom'])
continue
moteur = QgsGeometry.createGeometryEngine(coeur.constGet())
moteur.prepareGeometry()
bb = coeur.boundingBox()
places = 0
essais = 0
while places < k and essais < k * ESSAIS_PAR_PLACETTE:
essais += 1
x = rng.uniform(bb.xMinimum(), bb.xMaximum())
y = rng.uniform(bb.yMinimum(), bb.yMaximum())
p = QgsGeometry.fromPointXY(QgsPointXY(x, y))
if not moteur.contains(p.constGet()):
continue
voisins = index.nearestNeighbor(QgsPointXY(x, y), 1)
if voisins:
vx, vy, _ = tires[voisins[0]]
if math.hypot(x - vx, y - vy) < dist_centres:
continue
fid = len(tires)
tires.append((x, y, s['nom']))
f = QgsFeature(fid)
f.setGeometry(p)
index.addFeature(f)
places += 1
bilan.append('%s : %d / %d placettes' % (s['nom'], places, k))
if not tires:
message('Tirage impossible', '\n'.join(bilan), True)
return
# 4. numerotation du nord au sud et creation des couches
ordre = sorted(range(len(tires)), key=lambda i: -tires[i][1])
crs = zone.crs().authid()
champs = [QgsField('id_placette', QVariant.String), QgsField('strate', QVariant.String),
QgsField('X', QVariant.Double), QgsField('Y', QVariant.Double),
QgsField('forme', QVariant.String), QgsField('dim_m', QVariant.Double),
QgsField('surface_m2', QVariant.Double), QgsField('graine', QVariant.Int)]
centres = QgsVectorLayer('Point?crs=%s' % crs, 'centres', 'memory')
placettes = QgsVectorLayer('Polygon?crs=%s' % crs, 'placettes', 'memory')
for lyr in (centres, placettes):
lyr.dataProvider().addAttributes(champs)
lyr.updateFields()
fc, fp = [], []
for rang, i in enumerate(ordre, 1):
x, y, nom = tires[i]
attrs = ['P%02d' % rang if len(tires) < 100 else 'P%03d' % rang, nom, round(x, 1), round(y, 1),
forme, dim, round(surface_placette, 2), graine]
a = QgsFeature(centres.fields()); a.setAttributes(attrs)
a.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(x, y))); fc.append(a)
b = QgsFeature(placettes.fields()); b.setAttributes(attrs)
b.setGeometry(forme_placette(x, y, forme, dim)); fp.append(b)
centres.dataProvider().addFeatures(fc)
placettes.dataProvider().addFeatures(fp)
# 5. enregistrement dans un GeoPackage date, dans le dossier du projet
chemin = os.path.join(dossier, 'placettes_%s.gpkg' % datetime.now().strftime('%Y%m%d_%H%M'))
premier = True
for lyr in (centres, placettes):
opt = QgsVectorFileWriter.SaveVectorOptions()
opt.driverName = 'GPKG'
opt.layerName = lyr.name()
opt.actionOnExistingFile = (QgsVectorFileWriter.CreateOrOverwriteFile if premier
else QgsVectorFileWriter.CreateOrOverwriteLayer)
QgsVectorFileWriter.writeAsVectorFormatV3(lyr, chemin, projet.transformContext(), opt)
premier = False
for nom in ('placettes', 'centres'):
couche = QgsVectorLayer('%s|layername=%s' % (chemin, nom), nom, 'ogr')
if nom == 'centres':
reglage = QgsPalLayerSettings()
reglage.fieldName = 'id_placette'
reglage.enabled = True
couche.setLabeling(QgsVectorLayerSimpleLabeling(reglage))
couche.setLabelsEnabled(True)
projet.addMapLayer(couche)
# 6. controle : distance entre centres voisins
mini = None
for i, (x, y, _) in enumerate(tires):
for j, (x2, y2, _) in enumerate(tires):
if i < j:
d = math.hypot(x - x2, y - y2)
mini = d if mini is None else min(mini, d)
lignes = ['%d placettes tirées sur %d demandées (graine %d).' % (len(tires), n, graine),
'Placette : %s de %s m (%.1f m²).' % (forme.lower(), ('%g' % dim) if forme == 'Carré'
else 'rayon %g' % dim, surface_placette),
'Zone : %.2f ha ; taux d\'échantillonnage : %.2f %%.' % (surface_totale / 1e4,
100.0 * len(tires) * surface_placette / surface_totale)]
if mini is not None:
lignes.append('Distance minimale entre centres : %.1f m (exigée : %.1f m).' % (mini, dist_centres))
lignes += bilan
lignes.append('Fichier : %s' % chemin)
if len(tires) < n:
lignes.append("Toutes les placettes n'ont pas pu être placées : réduire le nombre, "
"la taille ou l'espacement, et le noter dans le protocole.")
message('Tirage terminé', '\n'.join(lignes))
try:
tirer()
except Annule:
print('Tirage annulé.')
A3. trier_photos_vol.py
Afficher le code (300 lignes)
# -*- coding: utf-8 -*-
r"""
Tri des photos d'une mission drone DJI (orthophoto / placettes / ecartees)
==========================================================================
A executer dans QGIS : Console Python > Afficher l'editeur > ouvrir > Executer.
Prerequis (projet QGIS ouvert, EPSG:2154) :
- couche de points 'decollage' (point de decollage REEL de la mission) ;
- couche de polygones 'zone_vol' ;
- couche de points 'centres' (placettes, champ 'id_placette') - facultative ;
- un dossier de photos DJI (JPG) ;
- internet (service d'altimetrie IGN) OU une couche raster MNT dans le projet.
Pour chaque photo, le script lit dans les metadonnees DJI : position GPS,
hauteur relative au decollage (barometre), inclinaison de nacelle, cap. Il lit
l'altitude du sol sous la photo (RGE ALTI) et calcule la HAUTEUR REELLE
AU-DESSUS DU SOL :
h_sol = h_rel + (sol au decollage - sol sous la photo)
Classement :
ortho : nacelle verticale, h_sol ~ hauteur ortho, a moins de 15 m de la zone
placette : nacelle verticale, h_sol ~ hauteur placette, a moins de 8 m d'un centre
ecartee : decollage, transits, montees / descentes, atterrissage...
Pour chaque placette, la photo la plus proche du centre est marquee 'meilleure'.
Sorties (dans le dossier des photos) :
- tri_photos.csv + couche de points 'tri_photos' (couleur par usage) ;
- facultatif : COPIE des photos dans les sous-dossiers ortho\ et placettes\
(les originaux ne sont pas modifies), pour RealityScan / WebODM.
"""
import os
import re
import csv
import json
import math
import shutil
from qgis.core import (QgsProject, QgsVectorLayer, QgsRasterLayer, QgsGeometry, QgsPointXY,
QgsWkbTypes, QgsCoordinateReferenceSystem, QgsCoordinateTransform,
QgsBlockingNetworkRequest, QgsCategorizedSymbolRenderer, QgsRendererCategory,
QgsMarkerSymbol)
from qgis.PyQt.QtCore import QUrl
from qgis.PyQt.QtNetwork import QNetworkRequest
from qgis.PyQt.QtWidgets import QInputDialog, QMessageBox, QFileDialog
try:
PARENT = iface.mainWindow() # noqa: F821
except NameError:
PARENT = None
IGN = 'Service altimétrie IGN (internet)'
URL_IGN = ('https://data.geopf.fr/altimetrie/1.0/calcul/alti/rest/elevation.json'
'?lon=%s&lat=%s&resource=ign_rge_alti_wld&delimiter=|&zonly=true')
NACELLE_MIN = -85.0 # nacelle verticale si pitch <= -85 deg
TOLERANCE_H = 4.0 # m autour de la hauteur visee
DIST_ZONE = 15.0 # m : photo ortho a moins de 15 m de la zone de vol
DIST_CENTRE = 8.0 # m : photo de placette a moins de 8 m du centre
CHAMPS_XMP = ['RelativeAltitude', 'GimbalPitchDegree', 'FlightYawDegree', 'GpsLatitude', 'GpsLongitude']
class Annule(Exception):
pass
def demander_reel(titre, texte, defaut, mini, maxi):
v, ok = QInputDialog.getDouble(PARENT, titre, texte, defaut, mini, maxi, 1)
if not ok:
raise Annule()
return v
def demander_choix(titre, texte, choix, defaut=0):
v, ok = QInputDialog.getItem(PARENT, titre, texte, choix, defaut, False)
if not ok:
raise Annule()
return v
def message(titre, texte, erreur=False):
if PARENT is not None:
(QMessageBox.critical if erreur else QMessageBox.information)(PARENT, titre, texte)
print(('ERREUR - ' if erreur else '') + titre + ' : ' + texte)
def couche(nom, type_geom, obligatoire=True):
for l in QgsProject.instance().mapLayersByName(nom):
if isinstance(l, QgsVectorLayer) and l.geometryType() == type_geom:
return l
if obligatoire:
raise RuntimeError("couche '%s' introuvable dans le projet" % nom)
return None
def lire_xmp(chemin):
with open(chemin, 'rb') as f:
data = f.read(200000)
d, f_ = data.find(b'<x:xmpmeta'), data.find(b'</x:xmpmeta>')
if d < 0 or f_ < 0:
return {}
xmp = data[d:f_].decode('utf-8', 'ignore')
res = {}
for c in CHAMPS_XMP:
m = re.search(r'drone-dji:%s\s*=\s*"([^"]*)"' % c, xmp) or \
re.search(r'<drone-dji:%s>([^<]*)<' % c, xmp)
if m:
try:
res[c] = float(m.group(1).replace('+', ''))
except ValueError:
pass
return res
def altitudes_ign(pts_wgs):
z = []
for i in range(0, len(pts_wgs), 50):
lot = pts_wgs[i:i + 50]
url = URL_IGN % ('|'.join('%.7f' % p.x() for p in lot), '|'.join('%.7f' % p.y() for p in lot))
req = QgsBlockingNetworkRequest()
if req.get(QNetworkRequest(QUrl(url))) != QgsBlockingNetworkRequest.NoError:
raise RuntimeError('service IGN injoignable : ' + req.errorMessage())
rep = json.loads(bytes(req.reply().content()).decode('utf-8'))
z += [float(v['z'] if isinstance(v, dict) else v) for v in rep['elevations']]
return z
def altitudes_raster(r, pts_wgs, ctx):
tr = QgsCoordinateTransform(QgsCoordinateReferenceSystem('EPSG:4326'), r.crs(), ctx)
z = []
for p in pts_wgs:
v, ok = r.dataProvider().sample(tr.transform(p), 1)
z.append(float(v) if ok else -99999.0)
return z
def trier():
projet = QgsProject.instance()
ctx = projet.transformContext()
dec_l = couche('decollage', QgsWkbTypes.PointGeometry)
zone_l = couche('zone_vol', QgsWkbTypes.PolygonGeometry)
cen_l = couche('centres', QgsWkbTypes.PointGeometry, obligatoire=False)
crs = zone_l.crs()
wgs = QgsCoordinateReferenceSystem('EPSG:4326')
vers_l93 = QgsCoordinateTransform(wgs, crs, ctx)
vers_wgs = QgsCoordinateTransform(crs, wgs, ctx)
zone = QgsGeometry.unaryUnion([f.geometry() for f in zone_l.getFeatures()])
g_dec = QgsGeometry(next(dec_l.getFeatures()).geometry())
if dec_l.crs() != crs:
g_dec.transform(QgsCoordinateTransform(dec_l.crs(), crs, ctx))
centres = []
if cen_l is not None:
tr_c = QgsCoordinateTransform(cen_l.crs(), crs, ctx) if cen_l.crs() != crs else None
champ = 'id_placette' if 'id_placette' in cen_l.fields().names() else None
for f in cen_l.getFeatures():
g = QgsGeometry(f.geometry())
if tr_c:
g.transform(tr_c)
centres.append((str(f[champ]) if champ else 'P%d' % f.id(), g.asPoint()))
dossier = QFileDialog.getExistingDirectory(PARENT, 'Dossier des photos de la mission', projet.homePath())
if not dossier:
raise Annule()
h_ortho = demander_reel('Hauteurs visées', 'Hauteur ortho au-dessus du SOL (m) :', 30, 5, 120)
h_pl = demander_reel('Hauteurs visées', 'Hauteur des photos de placettes au-dessus du SOL (m) :', 15, 3, 60)
rasters = [l for l in projet.mapLayers().values()
if isinstance(l, QgsRasterLayer) and l.providerType() == 'gdal' and l.bandCount() == 1]
source = demander_choix('Altitude du sol', "Source de l'altitude du sol :",
[IGN] + [r.name() for r in rasters], 0)
copier = demander_choix('Copie', 'Copier les photos dans ortho\\ et placettes\\ ?\n'
'(les originaux ne sont pas modifiés)', ['Oui', 'Non'], 0) == 'Oui'
# 1. lecture des photos
photos = []
for nom in sorted(os.listdir(dossier)):
if not nom.lower().endswith(('.jpg', '.jpeg')):
continue
x = lire_xmp(os.path.join(dossier, nom))
if 'GpsLatitude' not in x or 'GpsLongitude' not in x or 'RelativeAltitude' not in x:
photos.append({'photo': nom, 'ok': False})
continue
photos.append({'photo': nom, 'ok': True, 'h_rel': x['RelativeAltitude'],
'nacelle': x.get('GimbalPitchDegree'), 'cap': x.get('FlightYawDegree'),
'wgs': QgsPointXY(x['GpsLongitude'], x['GpsLatitude'])})
valides = [p for p in photos if p['ok']]
if not valides:
message('Aucune photo', 'Aucune photo DJI avec position GPS dans ce dossier.', True)
return
print('%d photos, %d avec position GPS : lecture du sol...' % (len(photos), len(valides)))
# 2. altitude du sol
pts = [vers_wgs.transform(g_dec.asPoint())] + [p['wgs'] for p in valides]
try:
z = altitudes_ign(pts) if source == IGN else \
altitudes_raster(rasters[[r.name() for r in rasters].index(source)], pts, ctx)
except Exception as e:
message('Altitude du sol', 'Lecture impossible (%s).' % e, True)
return
z0 = z[0]
# 3. classement
for p, zz in zip(valides, z[1:]):
p['l93'] = vers_l93.transform(p['wgs'])
p['z_sol'] = zz
p['h_sol'] = p['h_rel'] + z0 - zz
g = QgsGeometry.fromPointXY(p['l93'])
p['dist_zone'] = 0.0 if zone.contains(g) else zone.distance(g)
p['placette'], p['dist_centre'] = '', None
if centres:
nom, c = min(centres, key=lambda t: t[1].distance(p['l93']))
p['placette_proche'], p['dist_centre'] = nom, c.distance(p['l93'])
verticale = p['nacelle'] is not None and p['nacelle'] <= NACELLE_MIN
if verticale and abs(p['h_sol'] - h_pl) <= TOLERANCE_H and \
p['dist_centre'] is not None and p['dist_centre'] <= DIST_CENTRE:
p['usage'], p['placette'] = 'placette', p['placette_proche']
elif verticale and abs(p['h_sol'] - h_ortho) <= TOLERANCE_H and p['dist_zone'] <= DIST_ZONE:
p['usage'] = 'ortho'
else:
p['usage'] = 'ecartee'
for p in photos:
if not p['ok']:
p['usage'] = 'ecartee'
# meilleure photo par placette = la plus proche du centre
for nom, _ in centres:
sel = [p for p in valides if p['placette'] == nom]
if sel:
min(sel, key=lambda p: p['dist_centre'])['meilleure'] = True
# 4. CSV
champs = ['photo', 'usage', 'placette', 'meilleure', 'h_rel', 'z_sol', 'h_sol', 'nacelle', 'cap',
'dist_centre', 'dist_zone', 'x_l93', 'y_l93']
sortie = os.path.join(dossier, 'tri_photos.csv')
with open(sortie, 'w', newline='', encoding='utf-8') as f:
w = csv.writer(f, delimiter=';')
w.writerow(champs)
for p in photos:
if not p['ok']:
w.writerow([p['photo'], 'ecartee'] + [''] * (len(champs) - 2))
continue
r1 = lambda v: '' if v is None else round(v, 1)
w.writerow([p['photo'], p['usage'], p['placette'], 'oui' if p.get('meilleure') else '',
r1(p['h_rel']), round(p['z_sol'], 2), r1(p['h_sol']), r1(p['nacelle']), r1(p['cap']),
r1(p['dist_centre']), r1(p['dist_zone']),
round(p['l93'].x(), 2), round(p['l93'].y(), 2)])
# 5. couche de points
uri = ('file:///%s?delimiter=;&xField=x_l93&yField=y_l93&crs=%s'
% (sortie.replace('\\', '/'), crs.authid()))
lyr = QgsVectorLayer(uri, 'tri_photos', 'delimitedtext')
if lyr.isValid():
cats = []
for val, coul in (('ortho', '#2b83ba'), ('placette', '#d7191c'), ('ecartee', '#bdbdbd')):
s = QgsMarkerSymbol.createSimple({'name': 'circle', 'color': coul, 'size': '2'})
cats.append(QgsRendererCategory(val, s, val))
lyr.setRenderer(QgsCategorizedSymbolRenderer('usage', cats))
projet.addMapLayer(lyr)
# 6. copie
if copier:
for u in ('ortho', 'placettes'):
os.makedirs(os.path.join(dossier, u), exist_ok=True)
for p in valides:
if p['usage'] == 'ortho':
shutil.copy2(os.path.join(dossier, p['photo']), os.path.join(dossier, 'ortho', p['photo']))
elif p['usage'] == 'placette':
cible = '%s_%s%s' % (p['placette'], 'MEILLEURE_' if p.get('meilleure') else '', p['photo'])
shutil.copy2(os.path.join(dossier, p['photo']), os.path.join(dossier, 'placettes', cible))
# 7. bilan
n = {u: sum(1 for p in photos if p['usage'] == u) for u in ('ortho', 'placette', 'ecartee')}
lignes = ['%d photos : %d ortho, %d placettes, %d écartées.' % (len(photos), n['ortho'], n['placette'],
n['ecartee']),
'Sol au décollage : %.1f m.' % z0]
o = [p['h_sol'] for p in valides if p['usage'] == 'ortho']
if o:
lignes.append('Ortho : hauteur au-dessus du sol %.1f à %.1f m (moyenne %.1f).'
% (min(o), max(o), sum(o) / len(o)))
if centres:
lignes.append('')
lignes.append('Placette : nb photos | meilleure : écart au centre, h_sol, cap')
for nom, _ in sorted(centres):
sel = [p for p in valides if p['placette'] == nom]
best = [p for p in sel if p.get('meilleure')]
if best:
b = best[0]
lignes.append('%s : %d | %.1f m, %.1f m, cap %.0f°' % (nom, len(sel), b['dist_centre'],
b['h_sol'], b['cap'] or 0))
else:
lignes.append('%s : AUCUNE photo' % nom)
lignes += ['', 'CSV : %s' % sortie]
if copier:
lignes.append('Copies : sous-dossiers ortho\\ et placettes\\')
message('Tri des photos', '\n'.join(lignes))
try:
trier()
except Annule:
print('Tri annulé.')
except RuntimeError as e:
message('Tri des photos', str(e), True)
A4. lire_hauteurs_photos_dji.py
Utilisable dans Thonny ou dans la console QGIS. À réserver au terrain plat ; sur pente, utiliser A3.
Afficher le code (128 lignes)
# -*- coding: utf-8 -*-
"""
Lecture des metadonnees de vol des photos DJI (Mini 4 Pro)
==========================================================
Utilisable dans Thonny OU dans la console Python de QGIS (Ouvrir > Executer).
Chaque photo DJI contient, dans ses metadonnees XMP, des donnees que
l'Explorateur Windows n'affiche pas :
- RelativeAltitude : hauteur par rapport au point de decollage (barometre)
- AbsoluteAltitude : altitude GPS
- GimbalPitchDegree : inclinaison de la nacelle (-90 = verticale)
- FlightYawDegree : cap du drone (0 = nord)
- GpsLatitude / GpsLongitude
Le script lit toutes les photos d'un dossier et ecrit photos_vol.csv dans ce
dossier, avec une colonne 'usage' :
ortho : nacelle verticale, hauteur proche de la hauteur ortho
placette : nacelle verticale, hauteur proche de la hauteur placettes
ecartee : decollage, transit, nacelle non verticale...
Dans QGIS, la couche de points 'photos_vol' est ajoutee automatiquement.
ATTENTION : les hauteurs sont relatives au DECOLLAGE, comme dans la mission.
Pour une mission sur relief, comparer a la colonne h_rel_decollage de la
couche mission_waypoints, pas a 30 m.
"""
import os
import re
import csv
# ------------------------------------------------------------- a regler
DOSSIER = r'C:\Users\pierr\Downloads\Test Drone\DJI_001'
NACELLE_MIN = -85.0 # nacelle consideree verticale si pitch <= -85 deg
# hauteurs relatives au decollage attendues (m) et tolerance
H_ORTHO = (23.0, 26.0) # essai : A = 23,5 m, B = 25,2 m
H_PLACETTE = None # ex. (8.0, 11.0) pour la mission complete
# ----------------------------------------------------------------------
CHAMPS = ['RelativeAltitude', 'AbsoluteAltitude', 'GimbalPitchDegree',
'GimbalYawDegree', 'FlightYawDegree', 'GpsLatitude', 'GpsLongitude']
def lire_xmp(chemin):
with open(chemin, 'rb') as f:
data = f.read(200000) # le XMP est en tete de fichier
debut = data.find(b'<x:xmpmeta')
fin = data.find(b'</x:xmpmeta>')
if debut < 0 or fin < 0:
return {}
xmp = data[debut:fin].decode('utf-8', 'ignore')
res = {}
for c in CHAMPS:
m = re.search(r'drone-dji:%s\s*=\s*"([^"]*)"' % c, xmp) or \
re.search(r'<drone-dji:%s>([^<]*)<' % c, xmp)
if m:
try:
res[c] = float(m.group(1).replace('+', ''))
except ValueError:
pass
return res
def dans(v, borne):
return borne is not None and v is not None and borne[0] <= v <= borne[1]
def lire_dossier(dossier):
lignes = []
for nom in sorted(os.listdir(dossier)):
if not nom.lower().endswith(('.jpg', '.jpeg')):
continue
x = lire_xmp(os.path.join(dossier, nom))
h, p = x.get('RelativeAltitude'), x.get('GimbalPitchDegree')
m = re.search(r'DJI_(\d{8})(\d{6})', nom)
heure = '%s:%s:%s' % (m.group(2)[:2], m.group(2)[2:4], m.group(2)[4:]) if m else ''
if p is None or p > NACELLE_MIN:
usage = 'ecartee'
elif dans(h, H_PLACETTE):
usage = 'placette'
elif dans(h, H_ORTHO):
usage = 'ortho'
else:
usage = 'ecartee'
lignes.append({'photo': nom, 'heure': heure, 'h_rel': h, 'alt_gps': x.get('AbsoluteAltitude'),
'nacelle': p, 'cap': x.get('FlightYawDegree'),
'lat': x.get('GpsLatitude'), 'lon': x.get('GpsLongitude'), 'usage': usage})
return lignes
def main():
lignes = lire_dossier(DOSSIER)
if not lignes:
print('Aucune photo JPG dans', DOSSIER)
return
sortie = os.path.join(DOSSIER, 'photos_vol.csv')
with open(sortie, 'w', newline='', encoding='utf-8') as f:
w = csv.DictWriter(f, fieldnames=list(lignes[0].keys()), delimiter=';')
w.writeheader()
w.writerows(lignes)
print('%-36s %8s %8s %8s %6s %s' % ('photo', 'heure', 'h_rel', 'nacelle', 'cap', 'usage'))
for l in lignes:
print('%-36s %8s %8s %8s %6s %s' % (l['photo'], l['heure'],
'' if l['h_rel'] is None else '%.1f' % l['h_rel'],
'' if l['nacelle'] is None else '%.0f' % l['nacelle'],
'' if l['cap'] is None else '%.0f' % l['cap'], l['usage']))
hs = [l['h_rel'] for l in lignes if l['h_rel'] is not None]
for u in ('ortho', 'placette', 'ecartee'):
sel = [l['h_rel'] for l in lignes if l['usage'] == u and l['h_rel'] is not None]
if sel:
print('%-9s : %3d photos, h_rel %.1f a %.1f m (moyenne %.1f)'
% (u, len(sel), min(sel), max(sel), sum(sel) / len(sel)))
if hs:
print('Toutes : %d photos, h_rel %.1f a %.1f m' % (len(lignes), min(hs), max(hs)))
print('CSV :', sortie)
try: # dans QGIS : couche de points
from qgis.core import QgsVectorLayer, QgsProject
uri = ('file:///%s?delimiter=;&xField=lon&yField=lat&crs=EPSG:4326'
% sortie.replace('\\', '/'))
couche = QgsVectorLayer(uri, 'photos_vol', 'delimitedtext')
if couche.isValid():
QgsProject.instance().addMapLayer(couche)
except ImportError:
pass
main()
A5. Petits scripts de console
Ils se collent directement dans la console Python de QGIS.
Déplacer le point de décollage (la couche garde un seul point) :
from qgis.core import QgsProject, QgsGeometry, QgsPointXY, QgsFeature
c = QgsProject.instance().mapLayersByName('decollage')[0]
c.startEditing()
for f in c.getFeatures():
c.deleteFeature(f.id())
p = QgsFeature(c.fields())
p.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(598120.0177, 6414116.2058)))
c.addFeature(p)
c.commitChanges()
print(c.featureCount(), 'point :', next(c.getFeatures()).geometry().asPoint())
Hauteur réelle au-dessus de l'herbe après un décollage décalé (vol d'essai du 21/10) :
import json
from qgis.core import QgsBlockingNetworkRequest
from qgis.PyQt.QtNetwork import QNetworkRequest
from qgis.PyQt.QtCore import QUrl
url = ('https://data.geopf.fr/altimetrie/1.0/calcul/alti/rest/elevation.json'
'?lon=1.7112378|1.7115045&lat=44.8187268|44.8191028'
'&resource=ign_rge_alti_wld&delimiter=|&zonly=true')
r = QgsBlockingNetworkRequest(); r.get(QNetworkRequest(QUrl(url)))
z = json.loads(bytes(r.reply().content()))['elevations']
print('sol décollage réel : %.1f m | sol point prévu : %.1f m' % (z[0], z[1]))
print('hauteur réelle au-dessus de l\'herbe en A ≈ %.1f m' % (30.2 + z[0] - z[1]))
Exporter la table des waypoints en CSV (en filtrant sur le dossier de la mission, car le projet contient deux couches mission_waypoints) :
import csv, os
c = [l for l in QgsProject.instance().mapLayersByName('mission_waypoints')
if 'mission_20261022_1039' in l.source()][0]
champs = ['index', 'type', 'etiquette', 'h_rel_decollage', 'h_sur_sol', 'z_sol', 'vitesse', 'stationnaire_s']
lignes = sorted(([f[k] for k in champs] for f in c.getFeatures()), key=lambda l: l[0])
print(len(lignes), 'waypoints')
sortie = os.path.join(os.path.dirname(c.source().split('|')[0]), 'waypoints.csv')
with open(sortie, 'w', newline='', encoding='utf-8') as f:
w = csv.writer(f, delimiter=';'); w.writerow(champs); w.writerows(lignes)
print('CSV :', sortie)
Aucun commentaire:
Enregistrer un commentaire
Tout commentaire nous engage ;)