forked from luc-girod/MicMacWorkflowsByLucGirod
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathSat-Pleiades-SPOT.sh
executable file
·322 lines (297 loc) · 11.6 KB
/
Sat-Pleiades-SPOT.sh
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
#!/bin/bash
# Workflow MICMAC for Pleiades and SPOT 6/7 images
#
# Modified from Luc Girod https://github.com/luc-girod/MicMacWorkflowsByLucGirod/blob/master/DroneNadir.sh
# Jules Fleury, SIGEO/CEREGE
# 08/2018
# Large rewrite on 2022-11-14 by Luc Girod (see git commits)
# add default values
EXTIM=TIF # Image file extension
EXTRPC=XML # RPC file extension
PREFIM=IMG # Prefix for image name
PREFRPC=RPC # Prefix for RPC name
DEG=0 # Degree of the polynomial
PROJ_set=false
# necessary to check in this file that the UTM zone is OK
use_Schnaps=false
SH= #postfix for homologous points in case of schnaps use
wait_for_mask=false
ZOOM=2 # Final zoom
gresol_set=false # ground resolution set
RESSIZE=10000 # RESOLUTION OF SUBSAMPLED IMAGE FOR TAPIOCA, FOR FULL IMAGE USE -1
orthob=0 #boolean for creation of orthophotomosaic
EPSG=32632 # Coordinate system EPSG code, !!!!!!!! be coherent with the WGS84toUTM.xml file !!!!!!
ResolOrtho=1
DoOrtho=0
DEMInit="None"
NamePrefix="SatPleiades"
SkipValid=false
echo "
********************************************
*** MicMac workflow for ***
*** Pleiades/SPOT6-7 satellites images ***
*** without GCPs ***
********************************************
"
#input arguments
while getopts "e:f:p:q:d:u:v:r:smz:g:o:i:n:bkh" opt; do
case $opt in
h)
echo " "
echo "Usage: $0 -e TIF -f XML -p IMG -q RPC -d 0 -c WGS84toUTM.xml -r 10000 -s -m -z 2 -g 1 -o -a 32632"
echo " -e EXTIM : image file extension (JPG, jpg, TIF, tif, png..., default=$EXTIM)."
echo " -f EXTRPC : RPC file extension (XML, xml, ..., default=$EXTRPC)."
echo " -p PREFIM : Prefix for image name (default=$PREFIM)"
echo " -q PREFRPC : Prefix for RPC name (default=$PREFRPC)"
echo " -d DEG : Degree of the polynomial (default=$DEG)"
echo " -u UTMZONE : UTM Zone of area of interest. Takes form 'NN +north(south)'"
echo " -v PROJ : PROJ.4 string for coordinate system of output (use if not UTM)"
echo " -r RESSIZE : Resolution of the subsampled image for tapioca, (default=$RESSIZE, for full image use -1)"
echo " -s : Use 'Schnaps' optimized homologous points (default=$use_Schnaps)"
echo " -m : Pause for Mask before correlation (default=$wait_for_mask)"
echo " -z ZOOM : Zoom Level (default=$ZOOM)"
echo " -g GRESOL : Output Ground resolution (in meters)(if not set, will be defined automatically)"
echo " -o : 0 - no Ortho, 1 - Ortho using all provided images, 2 - Use _P for geometry and _MS for Ortho (default=$orthob)"
echo " -i DEMInit : Name of initialization DEM (without suffix, must have a MicMac style XML descriptor as well)"
echo " -n NamePrefix : Prefix name for output (default=SatPleiades)"
echo " -b BringData : Run GatherAirbusImages.sh to get data in folder (default=false)"
echo " -k SkipValid : Skip the validation of parameters, JUST DO IT (default=false)"
echo " -h : displays this message and exits."
echo " "
exit 0
;;
e)
EXTIM=$OPTARG
;;
f)
EXTRPC=$OPTARG
;;
p)
PREFIM=$OPTARG
;;
q)
PREFRPC=$OPTARG
;;
d)
DEG=$OPTARG
;;
c)
CHSYSXML=$OPTARG
;;
r)
RESSIZE=$OPTARG
;;
s)
use_Schnaps=true
;;
m)
wait_for_mask=true
;;
z)
ZOOM=$OPTARG
;;
u)
PROJ="+proj=utm +zone=$OPTARG +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
proj_set=true
;;
v)
PROJ=$OPTARG
proj_set=true
;;
g)
gresol_set=true
GRESOL=$OPTARG
;;
i)
DEMInit=$OPTARG
;;
n)
NamePrefix=$OPTARG
;;
o)
orthob=$OPTARG
if [ "$orthob" = 1 ]; then
ImOrtho="$PREFIM(.*).$EXTIM"
ImMNT="$PREFIM(.*).$EXTIM"
ResolOrtho=1
DoOrtho=1
elif [ "$orthob" = 2 ]; then
ImOrtho="$PREFIM(.*_MS.*).$EXTIM"
ImMNT="$PREFIM(.*_P.*).$EXTIM"
ResolOrtho=0.25
DoOrtho=1
fi
;;
b)
# Run the script to get the files from folder structure
GatherAirbusImages.sh
;;
k)
SkipValid=true
;;
\?)
echo "Script : Invalid option: -$OPTARG" >&1
exit 1
;;
:)
echo "Script.sh: Option -$OPTARG requires an argument." >&1
exit 1
;;
esac
done
#check arguments and choose to continue or not
if [ "$proj_set" = false ]; then
echo "Projection system not set (use option -u or -v, -h more more info)"
exit 1
fi
selection=
until [ "$selection" = "1" -o "$SkipValid" = true ]; do
echo "
CHECK (carefully) PARAMETERS
- Image file extension : $EXTIM
- RPC file extension : $EXTRPC
- Prefix of image file name : $PREFIM
- Prefix of RPC file name : $PREFRPC
- Degree of polynomial : $DEG
- Target projection system : $PROJ
- Tapioca resolution : $RESSIZE
- Use Schnaps : $use_Schnaps
- Pause for mask : $wait_for_mask
- ZoomF : $ZOOM
- Output GSD : $GRESOL
- Orthophotomosaic type : $orthob
- Ininialization DEM : $DEMInit
- Prefix name for output : $NamePrefix
"
echo "
CHOOSE BETWEEN
1 - Continue with these parameters
0 - Exit program
2 - Help
"
echo -n "Enter selection: "
read selection
echo ""
case $selection in
1 ) echo "Let's process now" ; continue ;;
0 ) exit ;;
2 ) echo "
For help use : ./SatPleiades.sh -h \n" >&1
exit 1 ;;
* ) echo "
Only 0 or 1 are valid choices
For help use : ./SatPleiades.sh -h \n" >&1
exit 1 ;;
esac
done
# MICMAC PROCESSING #################################
if [ "$use_Schnaps" = true ]; then
echo "Using Schnaps!"
SH="_mini"
else
echo "Not using Schnaps!"
SH=""
fi
if [ "$orthob" = 0 ]; then
ImOrtho="$PREFIM(.*).$EXTIM"
ImMNT="$PREFIM(.*).$EXTIM"
fi
#Create SysPROJ.xml
rm SysPROJ.xml
echo "<SystemeCoord> " >> SysPROJ.xml
echo " <BSC> " >> SysPROJ.xml
echo " <TypeCoord> eTC_Proj4 </TypeCoord> " >> SysPROJ.xml
echo " <AuxR> 1 </AuxR> " >> SysPROJ.xml
echo " <AuxR> 1 </AuxR> " >> SysPROJ.xml
echo " <AuxR> 1 </AuxR> " >> SysPROJ.xml
echo " <AuxStr> "$PROJ" </AuxStr> " >> SysPROJ.xml
echo " " >> SysPROJ.xml
echo " </BSC> " >> SysPROJ.xml
echo "</SystemeCoord> " >> SysPROJ.xml
#convert RPC info from nominal to MicMac format
#(specify the degree of your polynomial + working coordinate system)
mm3d Convert2GenBundle "$PREFIM(.*).$EXTIM" "$PREFRPC\$1.$EXTRPC" RPC-d$DEG Degre=$DEG ChSys=SysPROJ.xml
MaltOri=RPC-d$DEG
if [ "$DEG" != 0 ]; then
#Find Tie points using all images
mm3d Tapioca All "$PREFIM(.*).$EXTIM" $RESSIZE
if [ "$use_Schnaps" = true ]; then
#filter TiePoints (better distribution, avoid clogging)
mm3d Schnaps .*$EXTIM MoveBadImgs=1
fi
#Bundle adjustment, compensation
mm3d Campari "$PREFIM(.*).$EXTIM" RPC-d$DEG RPC-d$DEG-adj SH=$SH
MaltOri=RPC-d$DEG-adj
fi
#HERE, MASKING COULD BE DONE!!!
if [ "$wait_for_mask" = true ]; then
mm3d Tarama $ImMNT $MaltOri
read -rsp $'Do the masking and Press any key to continue...\n' -n1 key
fi
#Correlation into DEM
if [ "$DEMInit" != "None" ]; then
cp $DEMInit* .
DEMInit=${DEMInit##*/}
if [ "$gresol_set" = true ]; then
mm3d Malt Ortho "$PREFIM(.*).$EXTIM" $MaltOri EZA=1 ZoomF=$ZOOM VSND=-9999 DefCor=0 Spatial=1 MaxFlow=1 ImOrtho=$ImOrtho ImMNT=$ImMNT DoOrtho=$DoOrtho ResolOrtho=$ResolOrtho DEMInitIMG=$DEMInit.tif DEMInitXML=$DEMInit.xml ZoomI=4 ResolTerrain=$GRESOL
gdal_calc.py -A MEC-Malt/Correl_STD-MALT_Num_4.tif --outfile=MEC-Malt/AutoMask_STD-MALT_Num_5.tif --calc="A>100"
else
mm3d Malt Ortho "$PREFIM(.*).$EXTIM" $MaltOri EZA=1 ZoomF=$ZOOM VSND=-9999 DefCor=0 Spatial=1 MaxFlow=1 ImOrtho=$ImOrtho ImMNT=$ImMNT DoOrtho=$DoOrtho ResolOrtho=$ResolOrtho DEMInitIMG=$DEMInit.tif DEMInitXML=$DEMInit.xml ZoomI=4
gdal_calc.py -A MEC-Malt/Correl_STD-MALT_Num_4.tif --outfile=MEC-Malt/AutoMask_STD-MALT_Num_5.tif --calc="A>100"
fi
else
if [ "$gresol_set" = true ]; then
mm3d Malt Ortho "$PREFIM(.*).$EXTIM" $MaltOri EZA=1 ZoomF=$ZOOM VSND=-9999 DefCor=0 Spatial=1 MaxFlow=1 ImOrtho=$ImOrtho ImMNT=$ImMNT DoOrtho=$DoOrtho ResolOrtho=$ResolOrtho ResolTerrain=$GRESOL
gdal_calc.py -A MEC-Malt/Correl_STD-MALT_Num_7.tif --outfile=MEC-Malt/AutoMask_STD-MALT_Num_7.tif --calc="A>100"
else
mm3d Malt Ortho "$PREFIM(.*).$EXTIM" $MaltOri EZA=1 ZoomF=$ZOOM VSND=-9999 DefCor=0 Spatial=1 MaxFlow=1 ImOrtho=$ImOrtho ImMNT=$ImMNT DoOrtho=$DoOrtho ResolOrtho=$ResolOrtho
gdal_calc.py -A MEC-Malt/Correl_STD-MALT_Num_7.tif --outfile=MEC-Malt/AutoMask_STD-MALT_Num_7.tif --calc="A>100"
fi
fi
#Merge orthophotos to create Orthomosaic
if [ $DoOrtho -eq 1 ]; then
mm3d Tawny Ortho-MEC-Malt RadiomEgal=0
fi
#Post Processing ######################################
echo "
********************************************
*** Post-processing ***
********************************************
"
mkdir OUTPUT
cd MEC-Malt
lastDEM=($(find . -regex '.*Z_Num[0-9]*_DeZoom[0-9]*_STD-MALT.tif' | sort -r))
lastNuageImProf=($(find . -regex '.*NuageImProf.*[0-9].xml' | sort -r))
lastmsk=($(find . -regex '.*AutoMask_STD-MALT_Num_[0-9]*\.tif' | sort -r))
lastcor=($(find . -regex '.*Correl_STD-MALT_Num_[0-9]*\.tif' | sort -r))
lastDEMstr="${lastDEM%.*}"
lastcorstr="${lastcor%.*}"
lastmskstr="${lastmsk%.*}"
cp $lastDEMstr.tfw $lastcorstr.tfw
cp $lastDEMstr.tfw $lastmskstr.tfw
# Converting MicMac output DEM to files with masked areas as nodata
echo "gdal_translate -a_srs "$PROJ" $lastDEM tmp_geo.tif"
gdal_translate -a_srs "$PROJ" $lastDEM tmp_geo.tif
echo "gdal_translate -a_srs "$PROJ" -a_nodata 0 $lastmsk tmp_msk.tif"
gdal_translate -a_srs "$PROJ" -a_nodata 0 $lastmsk tmp_msk.tif
echo "gdal_calc.py -A tmp_msk.tif -B tmp_geo.tif --outfile=../OUTPUT/"${NamePrefix}"_DEM_MICMAC.tif --calc=\"B*(A>0)\" --NoDataValue=-9999"
gdal_calc.py -A tmp_msk.tif -B tmp_geo.tif --outfile=../OUTPUT/${NamePrefix}_DEM_MICMAC.tif --calc="B*(A>0)" --NoDataValue=-9999
rm tmp_geo.tif tmp_msk.tif
cd ..
gdal_translate -a_srs "$PROJ" MEC-Malt/$lastcor OUTPUT/${NamePrefix}_CORR_MICMAC.tif -co COMPRESS=DEFLATE
# export Ortho
if [ $DoOrtho -eq 1 ]; then
if [ -f "./Ortho-MEC-Malt/Orthophotomosaic_Tile_0_0.tif" ]
then
cd Ortho-MEC-Malt
mosaic_micmac_tiles.py -filename Orthophotomosaic
cd ..
fi
gdal_translate -a_nodata 0 -a_srs "$PROJ" Ortho-MEC-Malt/Orthophotomosaic.tif OUTPUT/${NamePrefix}_ORTHOMOSAIC_MICMAC.tif -co COMPRESS=DEFLATE
fi
echo "
********************************************
*** Finished ***
*** Results are in OUTPUT folder ***
********************************************
"