# Contents of /trunk/doc/inversion/CookGravity.tex

Revision 4286 - (show annotations)
Thu Mar 7 04:28:11 2013 UTC (6 years, 5 months ago) by caltinay
File MIME type: application/x-tex
File size: 38035 byte(s)
Assorted spelling fixes.


 1 \chapter{Gravity Inversion}\label{Chp:cook:gravity inversion} 2 3 \section{Introduction} 4 5 In this part of the documentation we give an introduction on how to use the 6 \downunder module and the inversion driver functions to perform the inversion 7 of gravity and magnetic data. 8 The driver functions enable geologists and geophysicists to apply the 9 \downunder module quickly and in an easy way without requiring detailed 10 knowledge of the theory behind inversion or programming skills. 11 However, users who are interested in specializing or extending the inversion 12 capacity are referred to Part~\ref{part2} of this manual. 13 It is beyond the intention of this manual to give a detailed introduction to 14 geophysical inversion, in particular to the appropriate preprocessing of data 15 sets. 16 17 The \downunder module described here is designed to calculate estimations for 18 the 3-D distribution of density and/or susceptibility from 2-D gravity and 19 magnetic data measured in ground or airborne surveys. 20 This process is generally called inversion of geophysical data. 21 Following the standard assumption it is assumed that the data are measured as 22 perturbation of an expected gravity and/or magnetic field of the Earth. 23 In this context measured gravity and magnetic data are in fact describing 24 anomalies in the gravity and magnetic field. 25 As a consequence the inversion process provides corrections to an average 26 density (typically $2670 kg/m^3$) and susceptibility (typically $0$). 27 So in the following we will always assume that given data are anomalies and 28 therefore not in all cases explicitly use the terms gravity anomalies or 29 density corrections but just use the terms gravity and density. 30 31 In this chapter we will give a detailed introduction into usage of the driver 32 functions for inversion of gravity data. 33 In the following chapters~\ref{Chp:cook:magnetic inversion} and~\ref{Chp:cook:joint inversion} 34 we will discuss the inversion of magnetic data, and the joint inversion of 35 gravity and magnetic data using \downunder. 36 Note, that the principles introduced in this chapter apply to magnetic and 37 joint inversion so the presentation for these problem classes is kept short 38 and users interested in magnetic data only should still work through this 39 chapter on gravity data. 40 41 To run the examples discussed you need to have \escript (version 3.3.1 or 42 newer) installed on your computer. 43 Moreover, if you want to visualize the results you need to have access to a 44 data plotting software which is able to process \VTK input files, e.g. 45 \mayavi or \VisIt. 46 As \mayavi can be easily obtained and installed for most platforms the 47 tutorial includes commands to visualize output files using \mayavi . 48 However, it is pointed out that \VisIt is the preferred visualization tool for 49 \escript as it can deal with very large data sets more efficiently. 50 51 \begin{figure} 52 \centering 53 \includegraphics[width=0.7\textwidth]{QLDWestGravityDataPlot.png} 54 \caption{Gravity Anomaly Data in $mgal$ from Western Queensland, Australia 55 (file \examplefile{data/QLDWestGravity.nc}). Data obtained from Geoscience Australia.} 56 %\AZADEH{Tracey R, Bacchin M, \& Wynne P. 2008. In preparation. AAGD07: A new absolute gravity datum for Australian gravity and new standards for the Australian National Gravity Database. Exploration Geophysics.}} 57 \label{FIG:P1:GRAV:0} 58 \end{figure} 59 60 \begin{figure} 61 \centering% 62 \includegraphics[width=0.9\textwidth]{QLDGravContourMu10.png} 63 \caption{3-D contour plot of the density distribution obtained by inversion of 64 file \file{data/QLDWestGravity.nc} (with $\mu=10$). 65 Colours represent values of density where high values are represented by 66 blue and low values are represented by red.} 67 \label{FIG:P1:GRAV:1} 68 \end{figure} 69 70 \section{How does it work?} 71 The execution of the inversion is controlled by a script which, in essence, is 72 a text file and can be edited using any text editor. 73 The script contains a series of statements which are executed by an 74 interpreter which is an executable program reading the text file and executing 75 the statements line-by-line. 76 In the case of \downunder the interpreter is \Python. 77 In order to be able to process the statements in each line of the script 78 certain rules (called syntax) need to be obeyed. 79 There is a large number of online tutorials for \Python available\footnote{e.g. 80 \url{http://www.tutorialspoint.com/python} and \url{http://doc.pyschools.com}}. 81 We also refer to the \escript cook book \cite{ESCRIPTCOOKBOOK} and user's 82 guide \cite{ESCRIPT} which is in particular useful for users who like to dive 83 deeper into \downunder. 84 For this part of the manual no \Python knowledge is required but it is 85 recommended that users acquire some basic knowledge on \Python as they 86 progress in their work with \downunder. 87 88 The following script~\ref{code: gravity1}\footnote{The script is similar to 89 \examplefile{grav_netcdf.py} within the \escript example file directory.} is a 90 simple example to run an inversion for gravity data: 91 92 \begin{pyc}\label{code: gravity1} 93 \ 94 \begin{python} 95 # Header: 96 from esys.downunder import * 97 from esys.weipa import * 98 from esys.escript import unitsSI as U 99 100 # Step 1: set up domain 101 dom=DomainBuilder() 102 dom.setVerticalExtents(depth=40.*U.km, air_layer=6.*U.km, num_cells=25) 103 dom.setFractionalPadding(pad_x=0.2, pad_y=0.2) 104 dom.fixDensityBelow(depth=40.*U.km) 105 106 # Step 2: read gravity data 107 source0=NetCdfData(NetCdfData.GRAVITY, 'GravitySmall.nc') 108 dom.addSource(source0) 109 110 # Step 3: set up inversion 111 inv=GravityInversion() 112 inv.setSolverTolerance(1e-4) 113 inv.setSolverMaxIterations(50) 114 inv.setup(dom) 115 116 # Step 4: run inversion 117 inv.getCostFunction().setTradeOffFactorsModels(10.) 118 rho = inv.run() 119 120 # Step 5: write reconstructed density to file 121 saveVTK("result.vtu", density=rho) 122 \end{python} 123 \end{pyc} 124 The result, in this case the density distribution, is written to an external 125 file for further processing. You can copy and paste the text of the script 126 into a file of any name, let's say for further reference we use the file 127 name \file{grav.py}. 128 It is recommended to use the extension \file{.py} to identify the file as a 129 \Python script. 130 We will discuss the statements of the script later in this chapter. 131 132 Obviously the inversion needs to be fed with some gravity data. You can find 133 example data from western Queensland, Australia in two resolutions in the 134 \escript example directory. In this case the data are loaded from the file 135 \examplefile{GravitySmall.nc} which is given in the \netcdf file format. 136 After you have copied this file into the directory in which you have saved the 137 script \file{grav.py} you can run the program using the command line 138 \begin{verbatim} 139 run-escript grav.py 140 \end{verbatim} 141 We are running \file{grav.py} through the \escript start-up command since 142 \escript is used as a back end for the inversion algorithm\footnote{Please see 143 the \escript user's guide~\cite{ESCRIPT} on how to run your script in parallel 144 using threading and/or MPI.}. 145 Obviously it is assumed that you have an installation of \escript available on 146 your computer, see \url{https://launchpad.net/escript-finley}. 147 148 After the execution has successfully completed you will find the result file 149 \file{result.vtu} in the directory from where you have started the execution 150 of the script. 151 The file has the \VTK format and can be imported easily into many 152 visualization tools. 153 One option is the \mayavi package which is available on most platforms. 154 You can invoke the visualization using the commands 155 \begin{verbatim} 156 mayavi2 -d result.vtu -m SurfaceMap 157 \end{verbatim} 158 from the command line. 159 Figure~\ref{FIG:P1:GRAV:1} shows the result of this inversion as a contour 160 plot\footnote{These plots were generated by \VisIt using the higher resolution 161 data.}, while the gravity anomaly data is shown in Figure~\ref{FIG:P1:GRAV:0}. 162 We will discuss data later in Section~\ref{SEC:P1:GRAV:DATA}. 163 164 Let us take a closer look at the script\footnote{In \Python lines starting 165 with \# are comments and are not processed further.}. Besides the header 166 section one can separate the script into five steps: 167 \begin{enumerate} 168 \item set up domain on which the inversion problem is solved 169 \item load the data 170 \item set-up the inversion problem 171 \item run the inversion 172 \item further processing of the result. Here we write the reconstructed 173 density distribution to a file. 174 \end{enumerate} 175 In the following we will discuss the steps of the scripts in more detail. 176 Before we do this it is pointed out that the header section, following 177 \Python conventions, makes all required packages available to access within 178 the script. 179 At this point we will not discuss this in more details but emphasize that the 180 header section is a vital part of the script. 181 It is is required in each \downunder inversion script and should not be 182 altered except if additional modules are needed. 183 184 \begin{figure} 185 \centering 186 \includegraphics[width=\textwidth]{dom2D.pdf} 187 \caption{2-D domain set-up for gravity inversion} 188 \label{FIG:P1:GRAV:2} 189 \end{figure} 190 191 \section{Creating the Inversion Domain} 192 First step in Script~\ref{code: gravity1} is the definition of the domain over 193 which the inversion is performed. 194 We think of the domain as a block with orthogonal, plain faces. 195 Figure~\ref{FIG:P1:GRAV:2} shows the set-up for a two-dimensional domain 196 (see also Figure~\ref{fig:domainBuilder} for 3-D). 197 The lateral coordinates along the surface of the Earth are denoted by $x$ and 198 $y$ (only $x$-direction is used in 2-D). 199 The $z$ direction defines the vertical direction where the part above the 200 surface has positive coordinates and the subsurface negative coordinates. 201 The height of the section above the surface, which is assumed to be filled 202 with air, needs to be set by the user. 203 The inversion assumes that the density in the section is known to be 204 zero\footnote{Always keeping in mind that these are not absolute values but 205 anomalies.}. 206 The density below the surface is unknown and is calculated through the 207 inversion. The user needs to specify the depth below the surface in which the 208 density is to be calculated. 209 The lateral extension of the domain is defined by the data sets fed into the 210 inversion. 211 It is chosen large enough to cover all data sets (in case more than one is 212 used). In order to reduce the impact of the boundary a padding zone around the 213 data sets can be introduced. 214 215 \begin{figure} 216 \centering 217 \includegraphics[width=\textwidth]{dom2DFEM.pdf} 218 \caption{Cell distribution and boundary conditions for a 2-D domain} 219 \label{FIG:P1:GRAV:3} 220 \end{figure} 221 222 The reconstruction of the gravity field from an estimated density distribution 223 is the key component of the inversion process. 224 To do this \downunder uses the finite element method (FEM). 225 We need to introduce a grid over the domain, see Figure~\ref{FIG:P1:GRAV:3}. 226 The number of vertical cells is set by the user while the number of horizontal 227 cells is derived from the grid spacing of the gravity data set(s). 228 It is assumed that gravity field data given are constant across a cell. 229 To be able to reconstruct the gravity field some assumptions on the values of 230 the gravity field on the domain boundary have to be made. 231 \downunder assumes that on all faces the lateral gravity field component 232 equals zero. No assumptions on the horizontal components are 233 made\footnote{It is assumed that the gravity potential equals zero on the top 234 and bottom surface, see Section~\ref{sec:forward gravity} for details}% 235 \footnote{Most inversion codes use Green's functions over an unbounded domain 236 to reconstruct the gravity field. This approach makes the assumption that the 237 gravity field (or potential) is converging to zero when moving away from the 238 region of interest. The boundary conditions used here are stronger in the 239 sense that the lateral gravity component is enforced to be zero in a defined 240 distance of the region of interest but weaker in the sense that no constraint 241 on the horizontal component is applied.}. 242 243 In script~\ref{code: gravity1} the statement 244 \begin{verbatim} 245 dom=DomainBuilder() 246 \end{verbatim} 247 creates something like a factory to build a domain. 248 We then define the features of the domain we would like to create: 249 \begin{verbatim} 250 dom.setVerticalExtents(depth=40.*U.km, air_layer=6.*U.km, num_cells=25) 251 dom.setFractionalPadding(pad_x=0.2, pad_y=0.2) 252 \end{verbatim} 253 Here we specify the depth of the domain to $40 km$, the thickness of the air 254 layer above the surface to $6 km$ and the number of vertical cells to $25$. 255 We also introduce a lateral padding of $20 \%$ of the expansion of the gravity 256 data on each side of the data and in both lateral directions. 257 258 In some cases it can be appropriate to assume that the density below a certain 259 depth is zero\footnote{As we are in fact calculating density corrections this 260 means that the density is assumed to be the average density.}. 261 The statement 262 \begin{verbatim} 263 dom.fixDensityBelow(depth=40.*U.km) 264 \end{verbatim} 265 introduces this constraint. 266 As in the case discussed here if the depth for zero density is not less than 267 the depth of the domain no constraint at depth is applied to the density. 268 269 \downunder uses the metre-kilogram-second based International System of Units 270 (SI)\footnote{see \url{http://en.wikipedia.org/wiki/International_System_of_Units}}\index{SI}. 271 So all values must be converted to appropriate units. 272 This does not apply to geographic coordinates which in \downunder are given in 273 fractional degrees (as a floating point number) to represent longitude and 274 latitude. In the script we have used the expression 275 \begin{verbatim} 276 depth=40.*U.km 277 \end{verbatim} 278 to define the depth of the domain to $40 km$. 279 The expression \verb|U.km| denotes the unit $km$ (kilometer) and ensures 280 appropriate conversion of the value $40$ into the base unit $m$ (meter). 281 It is recommended to add units to values (where present) in order to make sure 282 that the final values handed into \downunder is given with the appropriate 283 units. 284 The physical units module of \escript, which we have imported here under the 285 name \verb|U| in the script header, defines a large number of physical units 286 and constants, please see~\cite{ESCRIPT} and~\cite{ESCRIPTONLINE}. 287 288 \section{Loading Gravitational Data}\label{SEC:P1:GRAV:DATA} 289 In practice gravity acceleration is measured in various ways, for instance by 290 airborne surveys~\cite{Telford1990a}. 291 \downunder assumes that all data supplied as input are already appropriately 292 pre-processed. In particular, corrections for 293 \begin{itemize} 294 \item free-air, to remove effects from altitude above ground; 295 \item latitude, to remove effects from ellipsoidicity of the Earth; 296 \item terrain, to remove effects from topography 297 \end{itemize} 298 must have been applied to the data. 299 In general, data prepared in such a form are called Bouguer anomalies~\cite{Telford1990a}. 300 301 To load gravity data into \downunder the data are given on a plane parallel 302 to the surface of the Earth at a constant altitude, see 303 diagram~\ref{FIG:P1:GRAV:2}. 304 The data need to be defined over a rectangular grid covering a subregion of 305 the Earth surface. 306 The grid uses a geographic coordinate system with latitudes and longitudes 307 assumed to be given in the Clarke 1866 geodetic system. 308 Figure~\ref{FIG:P1:GRAV:0} shows an example of such a data set from the 309 western parts of Queensland, Australia. 310 The data set covers a rectangular region between $140^o$ and $141^o$ east 311 and between $20^o$ and $21^o$ south. 312 Notice that latitude varies between $-90^o$ to $90^o$ where negative signs 313 refer to places in the southern hemisphere and longitude varies between 314 $-180^o$ to $180^o$ where negative signs refer to places west of Greenwich. 315 The colour at a location represents the value of the vertical Bouguer gravity 316 anomaly at this point at the surface of the Earth. 317 Values in this data set range from $-160 \; mgal$ to about $500 \; mgal$\footnote{The unit 318 $mgal$ means milli $gal$ (galileo) with $1 \; gal = 0.01 \frac{m}{sec^2}$.} 319 over a $121 \times 121$ grid. 320 321 In general, a patch of gravity data needs to be defined over a plane 322 \verb|NX| $\times$ \verb|NY| where \verb|NX| and \verb|NY| define the number 323 of grid lines in the longitude (\verb|X|) and the latitude (\verb|Y|) 324 direction, respectively. 325 The grid is spanned from an origin with spacing \verb|DELTA_X| and 326 \verb|DELTA_Y| in the longitude and the latitude direction, respectively. 327 The gravity data for all grid points need to be given as an \verb|NX| 328 $\times$ \verb|NY| array. 329 If available, measurement errors can be associated with the gravity data. 330 The values are given as an \verb|NX| $\times$ \verb|NY| array matching the 331 shape of the gravity array. 332 Note that data need not be available on every single point of the grid, see 333 Section~\ref{SEC:P1:GRAV:REMARK:DATAHOLES} for more information on this. 334 335 Currently, two data file formats are supported, namely \emph{ER Mapper Raster} 336 \cite{ERMAPPER} files and \netcdf~\cite{NETCDF} files. 337 In the examples of this chapter we use \netcdf files and refer to 338 Section~\ref{SEC:P1:GRAV:REMARK:ERMAPPER} and 339 Section~\ref{sec:ref:DataSource:ERM} for more information on using ER Mapper 340 Raster files. 341 If you have data in any other format you have the option of writing a suitable 342 reader (for advanced users, see Chapter~\ref{Chp:ref:data sources}) or, 343 assuming you are able to read the data in \Python, refer to the example 344 script \examplefile{create_netcdf.py} which shows how to create a file 345 in the \netcdf file format~\cite{NETCDF} compatible with \downunder from 346 a data array. 347 348 In script~\ref{code: gravity1} we use the statement 349 \begin{verbatim} 350 source0=NetCdfData(NetCdfData.GRAVITY, 'GravitySmall.nc') 351 \end{verbatim} 352 to load the gravity data stored in \examplefile{GravitySmall.nc} in the 353 \netcdf format. 354 Within the script the data set is now available under the name \verb|source0|. 355 We need to link the data set to the \verb|DomainBuilder| using 356 \begin{verbatim} 357 dom.addSource(source0) 358 \end{verbatim} 359 At the time the domain for the inversion is built the \verb|DomainBuilder| 360 will use the information about origin, extent and spacing along with the other 361 options provided to build an appropriate domain. 362 As at this point a flat Earth is assumed geographic coordinates used to 363 represent data in the input file are mapped to a (local) Cartesian coordinate 364 system. This is achieved by projecting the geographic coordinates into the 365 \emph{Universal Transverse Mercator} (UTM) coordinate system\footnote{See e.g. 366 \url{http://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system}}. 367 368 There are a few optional arguments you can add when constructing a data source. 369 While Section~\ref{sec:ref:DataSource} has a detailed description of all 370 arguments it is worth noting a few. 371 Firstly, it is important to know that data slices are assumed to be at altitude 372 $0$ by default. This can be easily changed though: 373 \begin{verbatim} 374 source0=NetCdfData(NetCdfData.GRAVITY, 'GravitySmall.nc', altitude=2.5*U.km) 375 \end{verbatim} 376 Another important setting is the scale or unit of the measurements. 377 The default is dependent on the data type and for gravity anomalies a scale 378 of $\frac{\mu m}{sec^2}$ (or $0.1 \; mgal$) is assumed. 379 To change the default scale to $mgal$ (which is $10^{-5} \frac{m}{sec^2}$), 380 you could use: 381 \begin{verbatim} 382 source0=NetCdfData(NetCdfData.GRAVITY, 'GravitySmall.nc', scale_factor=1e-5) 383 \end{verbatim} 384 Finally, it is possible to specify measurement errors (i.e. uncertainties) 385 alongside the data. 386 Since these can never be zero, a value of $2$ units is used if nothing else 387 has been specified. 388 The error value is assumed to be given in the same units as the data so the 389 default value translates to an error of $0.2 \; mgal$. 390 There are two possibilities to specify the error, namely by providing a 391 constant value which is applied to all data points: 392 \begin{verbatim} 393 source0=NetCdfData(NetCdfData.GRAVITY, 'GravitySmall.nc', error=1.7) 394 \end{verbatim} 395 or, if the information is available in the same \netcdf file under the name 396 \verb|errors|, provide \downunder with the variable name: 397 \begin{verbatim} 398 source0=NetCdfData(NetCdfData.GRAVITY, 'GravitySmall.nc', error="errors") 399 \end{verbatim} 400 401 It is important to keep an eye on the complexity of the inversion. 402 A good measure is the total number of cells being used. 403 Assume we have given a data set on a $20 \times 20$ grid and we add lateral 404 padding of, say, $20 \%$ to each side of the data, the lateral number of cells 405 becomes $(20\cdot 1.4)\times (20\cdot 1.4)=1.4^2\cdot 20^2\approx 2\cdot 10^2=800$. 406 If we use $20$ cells in the vertical direction we end up with a total number 407 of $800 \times 20 = 16,000$ cells. 408 This size can be easily handled by a modern desktop PC. 409 If we increase the grid size of the data to $40 \times 40$ points and use $40$ 410 cells in the vertical extent we get a total of $(2\cdot 40^2)\cdot 40=128,000$ 411 cells, a problem size which is considerably larger but can still be handled by 412 a desktop computer. 413 Taking this one step further, if the amount of data is increased to 414 $200\times 200$ points and we use $200$ cells in the vertical extent the 415 domain will contain $16,000,000$ ($16$ million) cells. 416 This scenario requires a computer with enough memory and (a) fast processor(s) 417 to run the inversion. 418 This estimate of complexity growth applies to the case where the increase 419 of data grid size is driven by an increase of resolution where it is 420 recommended to increase the vertical resolution in synch with the lateral 421 resolution. Note that if more than one data set is used the target resolution 422 will be the resolution of the finest data set (see also 423 Section~\ref{SEC:P1:GRAV:REMARK:MULTIDATA}). 424 In other cases the expansion of the region of interest drives an increase of 425 data grid size and the increase of total number of cells is less dramatic as 426 the vertical number of cells can remain constant while keeping a balanced 427 resolution in vertical and lateral direction. 428 429 \section{Setting up the Inversion and Running it} 430 We are now at step three of script~\ref{code: gravity1} in which the actual 431 inversion is set up. 432 First we create an empty inversion under the name \verb|inv|: 433 \begin{verbatim} 434 inv=GravityInversion() 435 \end{verbatim} 436 As indicated by the name we can use \verb|inv| to perform an inversion of 437 gravity data\footnote{\verb|GravityInversion| is a driver with a simplified 438 interface which is provided for convenience. See Part~\ref{part2} for more 439 details on how to write inversion scripts with more general functionality, e.g. 440 constraints.}. The inversion is an iterative process which sequentially 441 calculates updates to the density distribution in an attempt to improve the 442 match of the gravity field produced by the density distribution with the data. 443 Termination of the iteration is controlled by the tolerance which is set by 444 the user: 445 \begin{verbatim} 446 inv.setSolverTolerance(1e-4) 447 \end{verbatim} 448 Here we set the tolerance to $10^{-4}$, i.e. the iteration is terminated if 449 the maximum density correction is less than or equal to $10^{-4}$ relative to 450 the maximum value of estimated density anomaly. 451 In case the iteration does not converge a maximum number of iteration steps is 452 set: 453 \begin{verbatim} 454 inv.setSolverMaxIterations(50) 455 \end{verbatim} 456 If the maximum number of iteration steps (here $50$) is reached the iteration 457 process is aborted and an error message is printed. 458 In this case you can try to rerun the inversion with a larger value for the 459 maximum number of iteration steps. 460 If even for a very large number of iteration steps no convergence is achieved, 461 it is very likely that the inversion has not been set up properly. 462 463 The statement 464 \begin{verbatim} 465 inv.setup(dom) 466 \end{verbatim} 467 links the inversion with the domain and the data. 468 At this step -- as we are solving a gravity inversion problem -- only 469 gravitational data attached to the domain builder \verb|dom| are considered. 470 Internally a cost function $J$ is created which is minimized during the 471 inversion iteration. 472 It is a combination of a measure of the data misfit of the gravity field from 473 the given density distribution and a measure of the smoothness of the density 474 distribution. 475 The latter is often called the regularization term. 476 By default the gradient of density is used as the regularization term, see 477 also Section~\ref{SEC:P1:GRAV:REMARK:REG}. 478 Obviously, the result of the inversion is sensitive to the weighting between 479 the misfit and the regularization. 480 This trade-off factor $\mu$ for the misfit function is set by the following 481 statement: 482 \begin{verbatim} 483 inv.getCostFunction().setTradeOffFactorsModels(0.1) 484 \end{verbatim} 485 Here we set $\mu=0.1$. The statement \verb|inv.setup| must appear in the 486 script before setting the trade-off factor. 487 A small value for the trade-off factor $\mu$ will give more emphasis to the 488 regularization component and create a smoother density distribution. 489 A large value of the trade-off factor $\mu$ will emphasize the misfit more 490 and typically creates a better fit to the data and a rougher density 491 distribution. 492 It is important to keep in mind that the regularization reduces noise in the 493 date and in fact gives the problem a unique solution. 494 Consequently, the trade-off factor $\mu$ may not be chosen too large in order 495 control the noise on the solution and ensure convergence in the iteration 496 process. 497 498 We can now actually run the inversion: 499 \begin{verbatim} 500 rho = inv.run() 501 \end{verbatim} 502 The answer as calculated during the inversion is returned and can be accessed 503 under the name \verb|rho|. 504 As pointed out earlier the iteration process may fail in which case the 505 execution of the script is aborted with an error message. 506 507 \section{Taking a Look} 508 In the final step of script~\ref{code: gravity1} the calculated density 509 distribution is written to an external file. 510 A popular file format used by several visualization packages such as 511 \VisIt~\cite{VISIT} and \mayavi~\cite{MAYAVI} is the \VTK file format. 512 The result of the inversion which has been named \verb|rho| can be written to 513 the file \file{result.vtu} by adding the statement 514 \begin{verbatim} 515 saveVTK("result.vtu", density=rho) 516 \end{verbatim} 517 at the end of script. 518 The inversion solution is tagged with the name \verb|density| in the result 519 file, however any other name for the tag could be used. 520 As the format is text-based (as opposed to binary) \VTK files tend to be very 521 large and take compute time to create, in particular when it comes to large 522 numbers of cells ($>10^6$). 523 For large problems it is more efficient to use the \SILO file format~\cite{SILO}. 524 \SILO files tend to be smaller and are faster generated and read. 525 It is the preferred format to import results into the visualization program 526 \VisIt~\cite{VISIT} which is particularly suited for the visualization of 527 large data sets. 528 Inversion results can directly exported into \SILO files using the statement 529 \begin{verbatim} 530 saveSilo("result.silo", density=rho) 531 \end{verbatim} 532 replacing the \verb|saveVTK(...)| statement. 533 Similar to \VTK files the result \verb|rho| is tagged with the name 534 \verb|density| so it can be identified in the visualization program. 535 536 \begin{figure} 537 \begin{center} 538 \subfigure[$\mu=0.1$]{% 539 \label{FIG:P1:GRAV:10 MU01} 540 \includegraphics[width=0.45\textwidth]{QLDGravContourMu01.png} 541 }% 542 \subfigure[$\mu=1.$]{% 543 \label{FIG:P1:GRAV:10 MU1} 544 \includegraphics[width=0.45\textwidth]{QLDGravContourMu1.png} 545 }\\ % ------- End of the first row ----------------------% 546 \subfigure[$\mu=10.$]{% 547 \label{FIG:P1:GRAV:10 MU10} 548 \includegraphics[width=0.45\textwidth]{QLDGravContourMu10.png} 549 }% 550 \subfigure[$\mu=100.$]{% 551 \label{FIG:P1:GRAV:10 MU100} 552 \includegraphics[width=0.45\textwidth]{QLDGravContourMu100.png} 553 }\\ % ------- End of the second row ----------------------% 554 \subfigure[$\mu=1000.$]{% 555 \label{FIG:P1:GRAV:10 MU1000} 556 \includegraphics[width=0.45\textwidth]{QLDGravContourMu1000.png} 557 }% 558 \end{center} 559 \caption{3-D contour plots of gravity inversion results with data from 560 Figure~\ref{FIG:P1:GRAV:0} for various values of the model trade-off 561 factor $\mu$. Visualization has been performed in \VisIt. 562 \AZADEH{check images.}} 563 \label{FIG:P1:GRAV:10} 564 \end{figure} 565 566 \begin{figure} 567 \begin{center} 568 \subfigure[$\mu=0.1$]{% 569 \label{FIG:P1:GRAV:11 MU01} 570 \includegraphics[width=0.45\textwidth]{QLDGravDepthMu01.png} 571 }% 572 \subfigure[$\mu=1.$]{% 573 \label{FIG:P1:GRAV:11 MU1} 574 \includegraphics[width=0.45\textwidth]{QLDGravDepthMu1.png} 575 }\\ % ------- End of the first row ----------------------% 576 \subfigure[$\mu=10.$]{% 577 \label{FIG:P1:GRAV:11 MU10} 578 \includegraphics[width=0.45\textwidth]{QLDGravDepthMu10.png} 579 }% 580 \subfigure[$\mu=100.$]{% 581 \label{FIG:P1:GRAV:11 MU100} 582 \includegraphics[width=0.45\textwidth]{QLDGravDepthMu100.png} 583 }\\ % ------- End of the second row ----------------------% 584 \subfigure[$\mu=1000.$]{% 585 \label{FIG:P1:GRAV:11 MU1000} 586 \includegraphics[width=0.45\textwidth]{QLDGravDepthMu1000.png} 587 }% 588 \end{center} 589 \caption{3-D slice plots of gravity inversion results with data from 590 Figure~\ref{FIG:P1:GRAV:0} for various values of the model trade-off 591 factor $\mu$. Visualization has been performed \VisIt. 592 \AZADEH{check images.}} 593 \label{FIG:P1:GRAV:11} 594 \end{figure} 595 596 Figures~\ref{FIG:P1:GRAV:10} and~\ref{FIG:P1:GRAV:11} show two different 597 styles of visualization generated in \VisIt using the result of the inversion 598 of the gravity anomalies shown in Figure~\ref{FIG:P1:GRAV:0}. 599 The inversions have been performed with different values for the model 600 trade-off factor $\mu$. 601 The visualization shows clearly the smoothing effect of lower values for the 602 trade-off factors. 603 For larger values of the trade-off factor the density distribution becomes 604 rougher showing larger details. 605 Computational costs are significantly higher for larger trade-off factors. 606 Moreover, noise in the data has a higher impact on the result. 607 Typically several runs are required to adjust the value for the trade-off 608 factor to the datasets used. 609 610 For some analysis tools it is useful to process the results in form of 611 Comma-separated Values (\CSV)\footnote{see 612 \url{http://en.wikipedia.org/wiki/Comma-separated_values}}. 613 Such a file can be created using the statement 614 \begin{verbatim} 615 saveDataCSV("result.csv", x=rho.getFunctionSpace().getX(), density=rho) 616 \end{verbatim} 617 in the script. 618 This will create a \file{result.csv} with columns separated by a comma. 619 Each row contains the value of the density distribution and the three 620 coordinates of the corresponding location in the domain. 621 There is a header specifying the meaning of the corresponding column. 622 Notice that rows are not written in a particular order and therefore, if 623 necessary, the user has to apply appropriate sorting of the rows. 624 Columns are written in alphabetic order of their corresponding tag names. 625 For the interested reader: the statement \verb|rho.getFunctionSpace()| returns 626 the type used to store the density data \verb|rho|. 627 The \verb|getX()| method returns the coordinates of the sampling points used 628 for the particular type of representation, see~\cite{ESCRIPT} for details. 629 630 \section{Remarks} 631 632 \subsection{ER Mapper Raster Files}\label{SEC:P1:GRAV:REMARK:ERMAPPER} 633 The \downunder module can read data stored in ER Mapper Raster files. A data 634 set in this format consists of two files, a header file whose name usually ends 635 in \verb|.ers| and the actual data file which has the same filename as the 636 header file but without any file extension. 637 These files are usually produced by a commercial software package and the 638 contents can be quite diverse. 639 Therefore, it is not guaranteed that every data set is supported by \downunder 640 but the most common types of raster data should work\footnote{If your data 641 does not load please contact us through \url{https://launchpad.net/escript-finley}.}. 642 643 The interface for loading ER Mapper files is very similar to the \netcdf 644 interface described in Section~\ref{SEC:P1:GRAV:DATA}. 645 To load gravity data stored in the file pair\footnote{These files are available 646 in the example directory.} \examplefile{GravitySmall.ers} (the 647 header) and \examplefile{GravitySmall} (the data) without changing any of the 648 defaults use: 649 \begin{verbatim} 650 source0=ErMapperData(ErMapperData.GRAVITY, 'GravitySmall.ers') 651 \end{verbatim} 652 If your data set does not follow the default naming convention you can specify 653 the name of the data file explicitly: 654 \begin{verbatim} 655 source0=ErMapperData(ErMapperData.GRAVITY, 'GravitySmall.ers', 656 datafile='GravityData') 657 \end{verbatim} 658 Please note that there is no way for the reader to determine if the two files 659 really form a pair so make sure to pass the correct filenames when constructing 660 the reader object. 661 The same optional arguments explained in sections \ref{SEC:P1:GRAV:DATA} and 662 \ref{SEC:P1:GRAV:REMARK:DATAHOLES} are available for ER Mapper data sets. 663 However, due to the limitation of the file format only a constant error value 664 is supported. 665 666 \subsection{Data With Holes}\label{SEC:P1:GRAV:REMARK:DATAHOLES} 667 As described previously in this chapter input data is always given in the form 668 of a rectangular grid with constant cell size in each dimension. 669 However, there are cases when this is not necessarily the case. 670 Consider an onshore data set which includes parts of the offshore region as in 671 Figure~\ref{FIG:P1:GRAV:onshoreoffshore}. 672 The valid data in this example has a value range of about $-600$ to $600$ and 673 the inversion is to be run based on these values only, disregarding the 674 offshore region. 675 In order to achieve that, the offshore region is \emph{masked} by using a 676 constant value which is not found within the onshore area. 677 Figure~\ref{FIG:P1:GRAV:onshoreoffshore} clearly shows this masked area in 678 dark blue since a mask value of $-1000$ was used. 679 680 \begin{figure}[ht] 681 \centering 682 \includegraphics[width=0.45\textwidth]{onshore-offshore.png} 683 \caption{Plot of a rectangular gridded onshore data set that includes 684 offshore regions which have a value (here $-1000$) not found within the 685 real data (Bouguer anomalies in Tasmania, courtesy Geoscience Australia)} 686 \label{FIG:P1:GRAV:onshoreoffshore} 687 \end{figure} 688 689 The \netcdf conventions supported in \downunder include a standard way of 690 specifying such a mask value. 691 The example script \examplefile{create_netcdf.py} demonstrates how this is 692 accomplished in an easy way with any data. 693 If, for any reason, the mask value in the input file is invalid it can be 694 overridden via the \verb|null_value| argument when constructing the 695 \verb|NetCdfData| object: 696 \begin{verbatim} 697 source0=NetCdfData(NetCdfData.GRAVITY, 'data0.nc', null_value=-1000) 698 \end{verbatim} 699 In this example, all data points that have a value of $-1000$ are ignored and 700 not used in the inversion. 701 Please note that the special value \emph{NaN} (not-a-number) is sometimes used 702 for the purposes of masking in data sets. 703 Areas marked with this value are always disregarded in \downunder. 704 705 \subsection{Multiple Data Sets}\label{SEC:P1:GRAV:REMARK:MULTIDATA} 706 It is possible to run a single inversion using more than one input data set, 707 possibly in different file formats. 708 To do so, simply create the data sources and add them to the domain builder: 709 \begin{verbatim} 710 source0=NetCdfData(NetCdfData.GRAVITY, 'data0.nc') 711 source1=ErMapperData(ErMapperData.GRAVITY, 'data1.ers') 712 dom.addSource(source0) 713 dom.addSource(source1) 714 \end{verbatim} 715 However, there are some restrictions when combining data sets: 716 \begin{itemize} 717 \item Due to the coordinate transformation all data sets must be located in 718 the same UTM zone. If a single dataset crosses UTM zones only the zone 719 of the central longitude is used when projecting. 720 For example, if one data set lies mostly in zone 51 but contains areas 721 of zone 52, it is transformed using zone 51. In this case more data 722 from zone 51 can be added, but not from any other zone. 723 \item All data sets should have the same spatial resolution but this is not 724 enforced. Combining data with different resolution is currently 725 considered experimental but works best when the resolutions are 726 multiples of each other. For example if the first data set has a 727 resolution (or cell size) of $100$ metres and the second has a cell 728 size of $50$ metres then the target domain will have a cell size of 729 $50$ metres (the finer resolution) and each point of the coarse data 730 will occupy two cells (in the respective dimension). 731 \end{itemize} 732 733 \subsection{Regularization Term}\label{SEC:P1:GRAV:REMARK:REG} 734 The \verb|GravityInversion| class supports the following form for the 735 regularization: 736 \begin{equation} 737 \int w^{(0)} \cdot \rho^2 + w^{(1)}_0 \rho_{,0}^2 + w^{(1)}_1 \rho_{,1}^2 + w^{(1)}_2 \rho_{,2}^2\; dx 738 \end{equation} 739 where the integral is calculated across the entire domain. 740 $\rho$ represents the density distribution where $\rho_{,0}$ $\rho_{,1}$ and 741 $\rho_{,2}$ are the spatial derivatives of $\rho$ with respect to the two 742 lateral and the vertical direction, respectively. 743 $w^{(0)}$, $w^{(1)}_0$, $w^{(1)}_1$ and $w^{(1)}_2$ are weighting 744 factors\footnote{A more general form, e.g. spatially variable values for the 745 weighting factors, is supported, see Part~\ref{part2}}. 746 By default these are $w^{(0)}=0$, $w^{(1)}_0=w^{(1)}_1=w^{(1)}_2=1$. 747 Other weighting factors can be set in the inversion set-up. 748 For instance to set $w^{(0)}=10$, $w^{(1)}_0=w^{(1)}_1=0$ and $w^{(1)}_2=100$ 749 use the statement: 750 \begin{verbatim} 751 inv.setup(dom, w0=10, w1=[0,0,100]) 752 \end{verbatim} 753 It is pointed out that the weighting factors are rescaled in order to improve 754 numerical stability. Therefore the relative size of the weighting factors is 755 relevant and using 756 \begin{verbatim} 757 inv.setup(dom, w0=0.1, w1=[0,0,1]) 758 \end{verbatim} 759 would lead to the same regularization as the statement above. 760