MMASTER : Différence entre versions

Version du 31 mai 2017 à 15:51

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

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 [1]. When chosing the file format, choose GeoTiff (NOT hdf).

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

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

• WorkFlowMMASTER.sh -s AST_L1A_00302212004225011_20170204145459_24750 -z "4 +north"

• extra options (given here with default values) :

  -t 30 -f 1 -n false -c 0.7

Processing

The following code can be in a bash (.sh) file and directly ran.

#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."