MMASTER : Différence entre versions

De MicMac
Aller à : navigation, rechercher
(→‎Preparing the data)
(→‎processing the scene)
 
(7 révisions intermédiaires par le même utilisateur non affichées)
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 Girod et al, 2017]
  
 
==Requirement==
 
==Requirement==
Ligne 8 : Ligne 8 :
  
 
==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/ NASA's Earthdata]. When chosing the file format, choose GeoTiff (NOT hdf).
  
−
==Preparing the data==
+
Workflows bash file (for now only for Unix systems) are available here : https://github.com/luc-girod/MMASTER-workflows . You probably want to add the folder where you have those files to your path.
 +
 
 +
==Processing a single scene==
 +
 
 +
=== Preparing the data ===
 
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.
  
 
First extract the GeoTiff L1A to a folder named SceneName/RawData (ex AST_L1A_00302212004225011_20170204145459_24750/RawData).
 
First extract the GeoTiff L1A to a folder named SceneName/RawData (ex AST_L1A_00302212004225011_20170204145459_24750/RawData).
 +
 +
=== Processing the scene ===
  
 
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 -q 5</pre>
−
 
+
−
==Processing==
+
−
The following code can be in a bash (.sh) file and directly ran.
+
−
<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
+
== Batch processing ==
−
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
+
The ''RunMicMacAster_batch.sh" bash file can run MMASTER with the same UTM zone parameter for a folder full of zipped ASTER_L1A scenes.
−
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
+
== Notes on default parameters ==
−
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>
+
Default parameters are optimized for areas of low contrast (glaciers, sandy deserts...). If you are trying to get more details in your data, the correlation window size should be decreased (default "-q 5" means a 11*11 window). This can be done with the "-q" option, we found that "-q 2" gives the best results on well textured surfaces.

Version actuelle en date du 17 novembre 2017 à 12:11

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 : Girod et al, 2017

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 NASA's Earthdata. 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 . You probably want to add the folder where you have those files to your path.

Processing a single scene

Preparing the data

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).

Processing the scene

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 -q 5

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.

Notes on default parameters

Default parameters are optimized for areas of low contrast (glaciers, sandy deserts...). If you are trying to get more details in your data, the correlation window size should be decreased (default "-q 5" means a 11*11 window). This can be done with the "-q" option, we found that "-q 2" gives the best results on well textured surfaces.

Récupérée de « http://micmac.ensg.eu/index.php?title=MMASTER&oldid=2546 »