MMASTER : Différence entre versions

De MicMac
Aller à : navigation, rechercher
(Processing)
Ligne 1 : Ligne 1 :
 
==Description==
 
==Description==
This tutorial is there to give the workflow necessary to compute MMASTER DEMs from ASTER L1A data. The method is described in details in an upcoming paper (will be linked here as soon as it is published).
+
This tutorial is there to give the workflow necessary to compute MMASTER DEMs from ASTER L1A data. The method is described in details in the following paper : [http://dx.doi.org/10.3390/rs9070704].
  
 
==Requirement==
 
==Requirement==
Ligne 9 : Ligne 9 :
 
==Download==
 
==Download==
 
Any AST_L1A scene can work here, you can download them from [https://earthdata.nasa.gov/]. When chosing the file format, choose GeoTiff (NOT hdf).
 
Any AST_L1A scene can work here, you can download them from [https://earthdata.nasa.gov/]. When chosing the file format, choose GeoTiff (NOT hdf).
 +
 +
Workflows bash file (for now only for Unix systems) are available here : https://github.com/luc-girod/MMASTER-workflows
  
 
==Preparing the data==
 
==Preparing the data==
 +
 +
=== One scene ===
 
Examples are given using a scene named ''AST_L1A_00302212004225011_20170204145459_24750'' and the UTM4N cartographic projection.
 
Examples are given using a scene named ''AST_L1A_00302212004225011_20170204145459_24750'' and the UTM4N cartographic projection.
  
Ligne 16 : Ligne 20 :
  
 
Then from the folder where **SCENE_FOLDER** (ex AST_L1A_00302212004225011_20170204145459_24750/) is situated run :
 
Then from the folder where **SCENE_FOLDER** (ex AST_L1A_00302212004225011_20170204145459_24750/) is situated run :
*<pre>WorkFlowMMASTER.sh -s AST_L1A_00302212004225011_20170204145459_24750 -z "4 +north"</pre>
+
*<pre>WorkFlowMMASTER_onescene.sh -s AST_L1A_00302212004225011_20170204145459_24750 -z "4 +north"</pre>
 
*extra options are possible with this file (given here with default values) :  <pre>-t 30 -f 1 -n false -c 0.7</pre>
 
*extra options are possible with this file (given here with default values) :  <pre>-t 30 -f 1 -n false -c 0.7</pre>
  
==Processing==
+
=== Batch processing ===
The following code is the '''WorkflowMMASTER.sh''' file referenced above.
+
<pre>
+
 
+
#Fixed symboles
+
Nx="_3N.xml"
+
Bx="_3B.xml"
+
Nt="_3N.tif"
+
Bt="_3B.tif"
+
Bcor="_3B.tif_corrected.tif"
+
RPC="RPC_"
+
scene_set=0
+
utm_set=0
+
# add default values for ZoomF, RESTERR, CorThr and NoCorDEM
+
ZoomF=1
+
RESTERR=30
+
CorThr=0.7
+
NoCorDEM=false
+
 
+
while getopts "s:z:c:nf:t:h" opt; do
+
case $opt in
+
h)
+
echo "Run the second step in the MMASTER processing chain."
+
echo "usage: WorkflowMMASTER.sh -s SCENENAME -z 'UTMZONE' -f ZOOMF -t RASTERRES -h"
+
echo "    -s SCENENAME: Aster scenename/folder where data is located."
+
echo "    -z UTMZONE  : UTM Zone of area of interest. Takes form 'NN +north(south)'"
+
echo "    -c CorThr  : Correlation Threshold for estimates of Z min and max (optional, default : 0.7)"
+
echo "    -n NoCorDEM : Compute DEM with the uncorrected 3B image (computing with correction as well)"
+
echo "    -f ZOOMF    : Run with different final resolution  (optional; default: 1)"
+
echo "    -t RESTERR  : Run with different terrain resolution (optional; default: 30)"
+
echo "    -h   : displays this message and exits."
+
echo " "
+
exit 0
+
;;
+
n)
+
NoCorDEM=$OPTARG
+
;;
+
s)
+
name=$OPTARG
+
scene_set=1
+
;;
+
z)
+
UTM=$OPTARG
+
utm_set=1
+
;;   
+
c)
+
CorThr=$OPTARG
+
echo "CorThr set to $CorThr"
+
;;
+
f)
+
ZoomF=$OPTARG
+
;;
+
t)
+
RESTERR=$OPTARG
+
;;
+
\?)
+
echo "RunMicMacAster.sh: Invalid option: -$OPTARG" >&2
+
exit 1
+
;;
+
:)
+
echo " WorkflowMMASTER.sh: Option -$OPTARG requires an argument." >&2
+
exit 1
+
;;
+
esac
+
done
+
 
+
#Variable symboles
+
echo $name
+
echo $UTM
+
cd $name
+
pwd
+
cd RawData
+
pwd
+
mm3d SateLib ASTERGT2MM $name
+
cd ..
+
 
+
mm3d SateLib Aster2Grid $name$Bx 20 "+proj=utm +zone=$UTM +datum=WGS84 +units=m +no_defs" HMin=-500 HMax=9000 expDIMAP=1 expGrid=1
+
mm3d SateLib Aster2Grid $name$Nx 20 "+proj=utm +zone=$UTM +datum=WGS84 +units=m +no_defs" HMin=-500 HMax=9000 expDIMAP=1 expGrid=1
+
mm3d SateLib Aster2Grid "FalseColor_$name.xml" 20 "+proj=utm +zone=$UTM +datum=WGS84 +units=m +no_defs" HMin=-500 HMax=9000 expDIMAP=1 expGrid=1
+
 
+
mm3d Malt Ortho ".*$name(|_3N|_3B).tif" GRIBin ImMNT="$name(_3N|_3B).tif" MOri=GRID ZMoy=2500 ZInc=2500 ZoomF=8 ZoomI=32 ResolTerrain=30 NbVI=2 EZA=1 Regul=0.1 DefCor=$CorThr DoOrtho=0 DirMEC=MEC-Mini
+
 
+
gdalinfo -nomd -norat -noct -nofl -stats MEC-Mini/Z_Num6_DeZoom8_STD-MALT.tif > gdalinfo.txt
+
deminfo=$(grep -P 'Minimum+' gdalinfo.txt)
+
Min=$(echo $deminfo | cut -d, -f1 | tr -d ' ' | tr -d 'Minimum=' | xargs printf "%.0f")
+
Max=$(echo $deminfo | cut -d, -f2 | tr -d ' ' | tr -d 'Maximum=' | xargs printf "%.0f")
+
 
+
echo Min=$Min
+
echo Max=$Max
+
 
+
#Filter obvious error in min/max (limit to earth min/max)
+
Min=$((($Min)<-420 ? -420 : $Min))
+
Max=$((($Max)>8850 ? 8850 : $Max))
+
#next 2 lines is basically if the auto min/max function failed / DEM is really bad, happen if a lot of sea or a lot of clouds
+
Min=$((($Min)>8850 ? -420 : $Min))
+
Max=$((($Max)<-420 ? 8850 : $Max))
+
#From min/max, compute the nb of grids needed in Z and the values for ZMoy and Zinc
+
DE=$(echo $Max - $Min| bc )
+
NbLvl=$(echo $DE/200| bc )
+
NbLvl=$((($NbLvl)<10 ? 10 : $NbLvl))
+
Mean=$(echo $Max + $Min| bc )
+
Mean=$(echo $Mean/2| bc )
+
Inc=$(echo $Max - $Mean| bc | xargs printf "%.0f")
+
echo Min=$Min
+
echo Max=$Max
+
echo NbLvl=$NbLvl
+
echo Mean=$Mean
+
echo Inc=$Inc
+
echo Min Max NbLvl Mean Inc >> Stats.txt
+
echo $Min $Max $NbLvl $Mean $Inc >> Stats.txt
+
 
+
#Re compute RPCs with updated min/max
+
mm3d SateLib Aster2Grid $name$Bx $NbLvl "+proj=utm +zone=$UTM +datum=WGS84 +units=m +no_defs" HMin=$Min HMax=$Max expDIMAP=1 expGrid=1
+
mm3d SateLib Aster2Grid $name$Nx $NbLvl "+proj=utm +zone=$UTM +datum=WGS84 +units=m +no_defs" HMin=$Min HMax=$Max expDIMAP=1 expGrid=1
+
mm3d SateLib Aster2Grid "FalseColor_$name.xml" $NbLvl "+proj=utm +zone=$UTM +datum=WGS84 +units=m +no_defs" HMin=$Min HMax=$Max expDIMAP=1 expGrid=1
+
 
+
mm3d MMTestOrient $name$Bt $name$Nt GRIBin PB=1 MOri=GRID ZoomF=1 ZInc=$Inc ZMoy=$Mean
+
 
+
 
+
 
+
# if we want to compute the uncorrected DEM
+
if [ "$NoCorDEM" = true ]; then #check variable name!
+
mm3d Malt Ortho ".*$name(|_3N|_3B).tif" GRIBin ImMNT="$name(_3N|_3B).tif" ImOrtho="FalseColor_$name.tif" MOri=GRID ZInc=$Inc ZMoy=$Mean ZoomF=1 ZoomI=32 ResolTerrain=30 NbVI=2 EZA=1 DefCor=0 Regul=0.1 ResolOrtho=2 DirMEC=MEC-NoCor
+
fi
+
 
+
#Applying correction to the 3B image
+
mm3d SateLib ApplyParralaxCor $name$Bt GeoI-Px/Px2_Num16_DeZoom1_Geom-Im.tif FitASTER=1 ExportFitASTER=1
+
mkdir ImOrig
+
mv $name$Bt ImOrig/$name$Bt
+
mv $name$Bcor $name$Bt
+
 
+
# Correlation with corrected image
+
mm3d Malt Ortho ".*$name(|_3N|_3B).tif" GRIBin ImMNT="$name(_3N|_3B).tif" ImOrtho="FalseColor_$name.tif" MOri=GRID ZInc=$Inc ZMoy=$Mean ZoomF=$ZoomF ZoomI=32 ResolTerrain=$RESTERR NbVI=2 EZA=1 DefCor=0 Regul=0.1 ResolOrtho=2
+
 
+
 
+
mm3d GrShade MEC-Malt/Z_Num9_DeZoom1_STD-MALT.tif ModeOmbre=IgnE
+
mm3d Tawny Ortho-MEC-Malt/ RadiomEgal=0
+
mm3d Nuage2Ply MEC-Malt/NuageImProf_STD-MALT_Etape_9.xml Out=$name.ply Attr=Ortho-MEC-Malt/Orthophotomosaic.tif
+
 
+
#Postprocessing to have easy to use geotifs
+
mkdir -p OUTPUT
+
outdir=$(pwd)/OUTPUT
+
cd MEC-Malt
+
finalimgs=($(ls Z_Num*_DeZoom1_STD-MALT.tif))
+
finalmsks=($(ls AutoMask_STD-MALT_Num*.tif))
+
finalcors=($(ls Correl_STD-MALT_Num*.tif))
+
lastimg=${finalimgs[-1]}
+
lastmsk=${finalmsks[-1]}
+
lastcor=${finalcors[-1]}
+
# strip the extension
+
laststr="${lastimg%.*}"
+
maskstr="${lastmsk%.*}"
+
corrstr="${lastcor%.*}"
+
 
+
cp -v $laststr.tfw $maskstr.tfw
+
cp -v $laststr.tfw $corrstr.tfw
+
# now, assign the CRS we got to the mask, dem, and apply.
+
echo "Georeferencing correlation mask"
+
gdal_translate -a_nodata 0 -a_srs "+proj=utm +zone=$UTM +ellps=WGS84 +datum=WGS84 +units=m +no_defs" $lastcor $dir\_CORR.tif
+
echo "Creating temporary georeferenced DEM"
+
gdal_translate -a_srs "+proj=utm +zone=$UTM +ellps=WGS84 +datum=WGS84 +units=m +no_defs" $lastimg tmp_geo.tif
+
echo "Creating temporary georeferenced Mask"
+
gdal_translate -a_srs "+proj=utm +zone=$UTM +ellps=WGS84 +datum=WGS84 +units=m +no_defs" -a_nodata 0 $lastmsk tmp_msk.tif
+
cd ../
+
if [ -d "Ortho-MEC-Malt" ]; then
+
cd Ortho-MEC-Malt
+
echo "Creating double size correlation mask for ortho"
+
gdal_translate -tr 15 15 -a_srs "+proj=utm +zone=$UTM +ellps=WGS84 +datum=WGS84 +units=m +no_defs" -a_nodata 0 ../MEC-Malt/$lastmsk tmp_mskDouble.tif
+
echo "Creating temporary georeferenced ortho"
+
gdal_translate -a_srs "+proj=utm +zone=$UTM +ellps=WGS84 +datum=WGS84 +units=m +no_defs" Orthophotomosaic.tif tmp_V123.tif
+
cd ../
+
fi
+
cd MEC-Malt
+
# apply the mask
+
echo "Applying mask to georeferenced DEM"
+
gdal_calc.py -A tmp_msk.tif -B tmp_geo.tif --outfile=$dir\_Z.tif --calc="B*(A>0)" --NoDataValue=-9999
+
cp -v $dir\_Z.tif $outdir/$datestr #might be good to code orig. wd here.
+
gdaldem hillshade $dir\_Z.tif $outdir/$datestr/$dir\_HS.tif
+
gdal_calc.py -A $dir\_CORR.tif --outfile=$outdir/$datestr/$dir\_CORR.tif --calc="((A.astype(float)-127)/128)*100" --NoDataValue=-9999
+
#cp -v $dir\_CORR.tif $outdir/$datestr #might be good to code orig. wd here.
+
rm -v tmp_msk.tif tmp_geo.tif
+
rm -v $dir\_Z.tif $dir\_CORR.tif
+
cd ../
+
if [ -d "Ortho-MEC-Malt" ]; then
+
cd Ortho-MEC-Malt
+
gdal_calc.py -B tmp_mskDouble.tif -A tmp_V123.tif --outfile=$dir\_V123.tif --calc="((A!=255)*(A+1)+(A==255)*A)*(B>0)" --NoDataValue=0 --allBands=A
+
#Expression complicated to solve real 0 values not being NoData and 255 no being +1-ed to 0
+
rm -v tmp_V123.tif tmp_mskDouble.tif
+
cp -v $dir\_V123.tif $outdir/$datestr
+
rm -v $dir\_V123.tif
+
cd ../
+
fi
+
echo "MMASTER processing is complete."
+
  
</pre>
+
The ''RunMicMacAster_batch.sh" bash file can run MMASTER with the same UTM zone parameter for a folder full of zipped ASTER_L1A scenes.

Version du 6 novembre 2017 à 14:32

Description

This tutorial is there to give the workflow necessary to compute MMASTER DEMs from ASTER L1A data. The method is described in details in the following paper : [1].

Requirement

Part of the algorithms used in MMASTER rely on the ALGLIB library. Since it is licensed under the GNU GPL license, it cannot be part of the standard MicMac release (under CeCILL-B license). For that reason, it is necessary to compile MicMac from source, getting the IncludeALGLIB branch of the Git repository.

You also need to have Gdal installed on your system.

Download

Any AST_L1A scene can work here, you can download them from [2]. When chosing the file format, choose GeoTiff (NOT hdf).

Workflows bash file (for now only for Unix systems) are available here : https://github.com/luc-girod/MMASTER-workflows

Preparing the data

One scene

Examples are given using a scene named AST_L1A_00302212004225011_20170204145459_24750 and the UTM4N cartographic projection.

First extract the GeoTiff L1A to a folder named SceneName/RawData (ex AST_L1A_00302212004225011_20170204145459_24750/RawData).

Then from the folder where **SCENE_FOLDER** (ex AST_L1A_00302212004225011_20170204145459_24750/) is situated run :

  • WorkFlowMMASTER_onescene.sh -s AST_L1A_00302212004225011_20170204145459_24750 -z "4 +north"
  • extra options are possible with this file (given here with default values) :
    -t 30 -f 1 -n false -c 0.7

Batch processing

The RunMicMacAster_batch.sh" bash file can run MMASTER with the same UTM zone parameter for a folder full of zipped ASTER_L1A scenes.