Home > Contributions
Volume Calculation Pipeline for Planetary Geomorphology
Antoine Lucas
29/02/2008
__________________
Purpose:
Planetary geomorphology often requires some morphometrics parameters such as areas, distance and volume. The purpose of this paper is to introduce a simple volume calculation pipeline using GMT (Generic Mapping Tools) under UNIX. An example is here presented on ILD (Interior Layered Deposits) within Ganges Chasma, Mars.
Requierements:
- Workstation under any POSIX OS's (Unix, Linux, Mac OS X, BSD). With some modifications, MS windows would be supported.
- GMT ( > 3.0.x) installed
- DTM data in binary or ASCII format (e.g. MOLA, HRSC-DTM, HiRISE-DTM)
Contributions:
Suggestions, remarks or any improvement are always welcome.
_______________________________________________________________________
1. Introduction
ILD have been obsverved in many place in Valles Marineris structure [1]. Many processes have been proposed for the formation of such object on Mars. In order to understand the formation, volume calculation could be very usefull.
Fig.1 - Perspective view of ILD within Ganges Chasma (IR THEMIS Mosaic draped on MOLA DTM) [2]
2. Data ingestion
In our purpose, we need a topographic data. For this example we will use MOLA MEGRD Grid provided by NASA and USGS and available online (fig.2).
Fig.2 - MOLA DTMs croped on ILD area.
Currently data need to be transpose to meter (there are in degrees). MOLA data use a sphere in the referencing data with the following parameters :
A_AXIS_RADIUS = 3396.0 Km
B_AXIS_RADIUS = 3396.0 Km
C_AXIS_RADIUS = 3396.0 Km
GMT thus need to know the Martian ellipsoid params by editing a new file named "MARS-MOLA" and containing the ellipsoid (here is a sphere) data. This operation is done at the beginning of the GMT script (see below).
Afterwards, we determine boudaries contour of our feature (here ILD mount) so as to perform the volume calcultation. In our exemple we use grdvolume program provided by GMT.
This contour is viewable by uncomment "#' in the code below (in the display part). We finally choose a baseline level needed by grdvolume program.
#! /bin/csh
# GMT Script Shell
# Author: A. Lucas (lucas@ipgp.fr)
# -----------------------------------------------------------------------------
#
# Need GMT > 3.x.x.
# Runs under Unix,Linux, BSD, Mac OS X
# -----------------------------------------------------------------------------
# Data parameters:
# - Ellipsoid
echo "MARS-MOLA 2001 3396000 3396000 0" >! MARS-MOLA
gmtset ELLIPSOID=MARS-MOLA # Set Martian ellipsoid uses by MOLA
# - I/O files
set file_in=ild_mola # Input filename
set oline=outline.txt # Outline file
set file_out=ild_out # output map filename
# - Geographical params
set inc=0.0078125 # Space increment of the DTM grid
set range=309/313/-8.5/-6 # Range for data injection
# - Display params
set niv=100 # Level increment for contour lines
set t=500 # Level increment for bolding
set proj=I311/7i # Projection type (see gmt manpage)
set base=-3200 # Baseline level for volume calculation
set coord=f1a2 # Displaying boudaries
set az=120 # Azimuth for 3D view
set elev=45 # Elevation for 3D view
set color=220/220/220 # Color box setting
set xt=312.1
set yt=-8.22
set yt2=-8.32
set ulx=312
set uly=-8.15
set urx=313
set ury=-8.15
set dlx=312
set dly=-8.4
set drx=313
set dry=-8.4
set level=4000
# -----------------------------------------------------------------------------
echo "GMT script Shell"
echo "__________________"
echo " "
echo "[1/5]. Resempling data..."
xyz2grd $file_in.img -G$file_out.grd -I$inc -R$range -F -Zh -V
grd2cpt $file_out.grd -R$range -Cjet -Z > $file_out.cpt
grdgradient $file_out.grd -G$file_out.int -A0/270 -Ne0.6
echo " "
echo "[2/5]. Masking data..."
grdmask $oline -I$inc -R$range -Gmask.grd -F -N0/1/1
grdmask $oline -I$inc -R$range -Gmask2.grd -F -N0/3/3
grdsample mask.grd -I$inc -R$range -Gout.grd
grdsample mask2.grd -I$inc -R$range -Gout2.grd
mv out.grd mask.grd
mv out2.grd mask2.grd
grdmath mask.grd mask2.grd NAN = mask_tmp.grd
grdmath $file_out.grd mask_tmp.grd MUL = sub.grd
#grd2cpt sub.grd -R$range -Cjet -Z > sub.cpt
echo " "
echo "[3/5]. Volume calculation..."
set info=`grdvolume sub.grd -C$base -S`
echo " "
echo "[4/5]. Postscript output..."
grdimage $file_out.grd -I$file_out.int -C$file_out.cpt -J$proj \
-B$coord -R$range -K -Y5 >! $file_out.ps
psscale -D3.5i/-0.5i/3i/0.2ih -C$file_out.cpt -B1500/:"Elevation (m)": -O \
-K >> $file_out.ps
grdcontour $file_out.grd -C$niv -R$range -J$proj -A$t -B$coord \
-K -O -U >> $file_out.ps
psxy -R -J$proj -O -K -G$color -L -Wthicker << END >> $file_out.ps
$ulx $uly
$urx $ury
$drx $dry
$dlx $dly
END
echo $xt $yt 8 0 1 5 Area: $info[2] m@+2@+| pstext -R$range -J$proj \
-O -K >> $file_out.ps
echo $xt $yt2 8 0 1 5 Volume: $info[3] m@+3@+ | pstext -R -J$proj \
-O -K >> $file_out.ps
#grdimage sub.grd -Csub.cpt -J$proj \
#-B$coord -R$range -K -Y3 >! mask.ps
#psscale -D3.5i/-0.5i/3i/0.2ih -Csub.cpt -B1500/:"Elevation (m)": -O \
# -K >> mask.ps
#grdcontour sub.grd -C$niv -R$range -J$proj -A$t -B$coord \
#-K -O -U >> mask.ps
rm *.grd
rm *.cpt
rm *.int
echo " "
echo "[5/5]. Convert PostScript to PDF using ghostscript"
ps2pdf $file_out.ps
ps2pdf mask.ps
echo " "
echo " "
echo "Done !"
echo ""
echo " Thanks Antwan :)"
echo ""
# -----------------------------------------------------------------------------
You can download the MOLA grid and the Outline file.
This GMT script gives you this result:
Fig.3 -Topographic map from GMT.
See also
Ross Beyer provides a very efficient PERL program named mola2gmt.
References
[1] S.S. Nedell, S.W. Squyres, and D.W. Andersen. Origin and evolution of layered deposits in the Valles Marineris, Mars. Icarus,70, 409-441, 1987.
[2] A. Lucas, Géologie de Ganges Chasma, Mars, Université Paris-Sud, Mémoire de Master, Juin 2005
[3] Some Free GIS software : GRASS - gvGIS - QuantumGIS - SAGA GIS
