Blažej Bucha

Senior Lecturer at the Department of Theoretical Geodesy, Slovak University of Technology in Bratislava.

This page offers output products from some of research papers that I co-authored, most of them being related to gravity field modelling. All codes are written in Matlab and were tested in Matlab R2015b and R2018b. Anyone is encouraged to modify the routines according to their needs.

In case you have any comments, suggestions, questions or criticism, please feel free to contact me at blazej.bucha@stuba.sk.

GrafLab (GRAvity Field LABoratory) — Ultra-high-degree surface and solid spherical harmonic synthesis

GrafLab is a MATLAB-based routine for computing functionals of the geopotential up to ultra-high degrees and orders (tens of thousands and beyond). It allows:

• evaluation of 38 functionals of the geopotential up to ultra-high degrees and orders (e.g., the geoid, height anomaly, gravity anomalies/disturbances, deflections of the vertical, gravitational tensor),
• evaluation of commission error of 26 functionals using full variance-covariance matrix of spherical harmonic coefficients,
• depiction of the computed data on a map.

The source code accompanied by a manual and test data can be downloaded here.

GrafLab can be fully controlled either through GUI (see the figure below) or as a function (without GUI). The latter option is beneficial especially when it needs to be called multiple times with, for instance, different geopotential model, maximum degree or evaluation points.

Further details can be found in the paper:

Bucha, B., Janák, J., 2013. A MATLAB-based graphical user interface program for computing functionals of the geopotential up to ultra-high degrees and orders. Computers & Geosciences 56, 186-196, http://doi.org/10.1016/j.cageo.2013.03.012.

Some sample outputs from GrafLab

• Geoid undulations evaluated and decipted via GrafLab from EGM96 and DTM2006.0 up to degree 360:

This output can be reproduced by using test data and running the following command:

GrafLab('OK',3986004.415E+8,6378136.3,0,360,1,'EGM96.mat',0,0,-90,0.25,90,0,0.25,360,0,[],[],[],[],'EGM96_Geoid_to360',0,10,1,'DTM2006.mat',0,1,1,2,6,1,60,600,1)

• Besides the evaluation of gravity field quantities such as the geoid, gravity anomalies or deflections of the vertical, GrafLab can also be used to synthesize, for instance, the Earth's topography (look for the text "Some useful tips and tricks" in the list_of_changes_in_GrafLab.txt file):

This output can be reproduced by using test data and running the following command:

GrafLab('OK',1,1,0,360,[0 0 0 0 0],'DTM2006.mat',1,0,-90,0.25,90,0,0.25,360,0,[],[],[],[],'Topography',0,11,1,[],0,1,1,2,6,1,60,600,1)

• In fact, this command can be used to perform any surface spherical harmonic synthesis of the form $$f(\varphi,\lambda) = \sum_{n=0}^{n_{\max}} \sum_{m=0}^{n} \left( \bar{C}_{nm}\, \cos(m\, \lambda) + \bar{S}_{nm}\, \sin(m\, \lambda) \right) \bar{P}_{nm}(\sin\varphi)\,,$$ $$f(\varphi,\lambda) = \sum_{n=0}^{n_{\max}} \sum_{m=0}^{n} \Bigl( \bar{C}_{nm}\, \cos(m\, \lambda)$$ $$+ \bar{S}_{nm}\, \sin(m\, \lambda) \Bigr) \, \bar{P}_{nm}(\sin\varphi)$$ where $\bar{C}_{nm}$ and $\bar{S}_{nm}$ are $4\pi$-fully-normalized (real) surface spherical harmonic coefficients of the function $f$, $n$ and $m$ are spherical harmonic degree and order, respectively, $\bar{P}_{nm}(\sin\varphi)$ are the $4\pi$-fully-normalized (real) associated Legendre functions of the first-kind, and, finally, $\varphi$ and $\lambda$ are the spherical latitude and longitude, respectively. This means that GrafLab can synthesize a wide range of (real) functions given on a sphere up to ultra-high harmonic degrees (tens of thousands and beyond), provided that the respective coefficients $\bar{C}_{nm}$ and $\bar{S}_{nm}$ are available.

• After some simple tricks (look for the text "Some useful tips and tricks" in the list_of_changes_in_GrafLab.txt file), GrafLab can be used to synthesize gravity field of other celestial objects, such as the Moon, Mars and others. Below is shown an example: the gravitational vector (in the local north-oriented reference frame) induced by topographic masses of the Moon. The evaluation points reside on a Brillouin sphere (out of all gravitating masses).

The gravitational vector was computed up to degree 2160 at 5 arc-min grid from the model "STU_Moon_topography_to720_gravity_to2160.gfc" (see the section on divergence-free gravity models of the Moon's gravity field below) using the following command:

GrafLab('OK',0.4902805823000e+13,0.1737999e+07,0,2160,[0 0 0 0 0],'STU_Moon_topography_to720_gravity_to2160.gfc',1,0,-90,0.08,90,0,0.08,360,11000,[],[],[],[],'Moon_gravitational_vector',0,16,3,[],0,1,1,2,6,1,60,600,1)

Computation time

Below are shown computation times required by GrafLab to synthesise the Earth's topography as represented by Earth2014.TBI2014.degree10800 up to harmonic degrees nmax = 90, 180, 360, 720, 1080, 2160, 4320, 6480, 8640 and 10800. The position of evaluation points is defined by the Gauss—Legendre quadrature (see Section Ultra-high-degree surface spherical harmonic analysis). The dimensions of the grids are (nmax+1)*(2*nmax+2) in terms of the latitude and the longitude, respectively, and Legendre's functions were evaluated after Fukushima (2012). The tests were conducted on a PC with Intel® CoreTM i7-6800K CPU, 128 GB of RAM (no more than a few GBs were needed in this example though) and a 250 GB SSD drive. The same grids and PC were used to demonstrate the computational times required by the software package to perform spherical harmonic analysis in Section Ultra-high-degree surface spherical harmonic analysis, so that the two figures can be mutually compared.

Go to top

isGrafLab (Irregular Surface GRAvity Field LABoratory) — Efficient ultra-high-degree solid spherical harmonic synthesis at regular grids residing on irregular surfaces

isGrafLab is a modified version of GrafLab that allows accurate and fast computation of functionals of the geopotential at dense grids on irregular surfaces such as the Earth's surface. It employs the highly efficient lumped coefficients approach for surface spherical harmonic synthesis of the functionals and their radial derivatives at regular surfaces (sphere or ellipsoid of revolution), and the Taylor series expansions to continue the functionals to the irregular surface. All the other options available in GrafLab, such as three approaches to compute Legendre's functions, the graphical user interface or data plotting, are also available in isGrafLab.

The Matlab source code accompanied by a manual and test data can be downloaded here.

Further details can be found in the paper:

Bucha, B., Janák, J., 2014. A MATLAB-based graphical user interface program for computing functionals of the geopotential up to ultra-high degrees and orders: Efficient computation at irregular surfaces. Computers & Geosciences 66, 219-227, http://doi.org/10.1016/j.cageo.2014.02.005.

Go to top

GOCE-only global gravity field models based on the kinematic orbit

The GOCE-only global spherical harmonic gravity field models models developed by Bucha et al. (2015) can be downloaded here (~1.8 MB).

The models were estimated by inverting the kinematic orbit of the GOCE satellite (time period from 01 November 2009 to 11 January 2010 which is equal to the R1 period used by the official ESA's models). Some further specifications:

• the models were originally derived by means of spherical radial basis functions (Shannon SRBF) and were transformed into spherical harmonics after the evaluation; the transformation is exact without any approximations,
• the acceleration approach was used to invert the kinematic orbit into gravity field models,
• the problematic near-zonal coefficients (polar gap) are stabilized by prior information obtained from our own low-degree GOCE-only models, so that the final models are still GOCE-only,
• the models are estimated up to relatively high harmonic degrees (75, 100, 130) to mitigate the aliasing, but it is not recommended to use them beyond degree, say, 70 due to the adverse signal-to-noise ratio,
• the models are provided in the gfc format as defined by ICGEM (International Centre for Global Earth Models).

The effect caused by the stabilization of the near-zonal coefficients is shown in terms of difference degree amplitudes in the figure below (Bucha et al., 2015).

When using the models, please use the following reference:

Bucha, B., Bezděk, A., Sebera, J., Janák, J., 2015. Global and regional gravity field determination from GOCE kinematic orbit by means of spherical radial basis functions. Surveys in Geophysics 36, 773-801, http://doi.org/10.1007/s10712-015-9344-0

Go to top

Divergence-free spherical harmonic models of the gravity field implied by the Moon's topographic masses

These models approximate the gravitational field induced by topographic masses of the Moon and mitigate the divergence effect of spherical harmonics in the vicinity of the Moon's topography. The models rely on the Runge—Krarup theorem and enable generally a more accurate evaluation of gravity field quantities close to the lunar topography as compared to the models from spectral gravity forward modelling. This is because the latter ones may suffer from the divergence effect when evaluating the spherical harmonic series on or below the limit sphere encompassing all the gravitating masses (that is, also on the topography). In the figure below, the two approaches are validated with respect to a third (reference) one (Bucha et al., 2019). The upper map shows the reference signal, the middle one depicts the differences related to the Runge—Krarup solution and the bottom one validates the spectral gravity forward modelling technique; note the different colour bars.

We provide in total of 12 models that approximate the gravity field implied by the lunar topographic masses. The topographic masses are expanded up to degrees 90, 180, 360 and 720, and the maximum degree of the gravity field models varies from 360 up to 2160. All models (~320 MB) are available for download here in the gfc format as defined by ICGEM. One of the models, "STU_Moon_topography_to720_gravity_to2160", can also be accessed from ICGEM, where it can be find under a shortened name "STU_MoonTopo720".

The corresponding models that may suffer from the divergence effect were derived by Hirt and Kuhn (2017) and can be downloaded here, including grids of reference gravity disturbances from spatial-domain Newtonian integration.

Please use the following reference when using the models in your works:

Bucha, B., Hirt, C., Kuhn, M., 2019. Divergence-free spherical harmonic gravity field modelling based on the Runge—Krarup theorem: a case study for the Moon. Journal of Geodesy 93, 489-513, https://doi.org/10.1007/s00190-018-1177-4

Go to top

Ultra-high-degree surface spherical harmonic analysis

These routines are designed for accurate and fast FFT-based ultra-high-degree (tens of thousands) surface spherical harmonic analysis. As their key features, the exact Gauss—Legendre quadrature is employed and Legendre's functions are evaluated after Fukushima (2012). Two versions are available:

• non-parallel, which is designed to harmonically analyze an input grid, while storing the grid, output spherical harmonic coefficients, Legendre's functions, etc., in RAM during the entire computation (requires tens of GBs of RAM when the maximum degree is equal to a few tens of thousands),
• parallel, which performs ultra-high-degree surface spherical harmonic analysis of integer powers of input function represented by spherical harmonic coefficients; this code is specifically tailored to HPC clusters, where only a few GBs of RAM per worker may be available.

As an example, shown below are accuracy levels achieved in a closed-loop test environment, where signals up to degrees 90, 180, 360, 720, 1080, 2160, 4320, 6480, 8640 and 10800 were harmonically analysed. These results were obtained by the "SHA_GL.m" routine (non-parallel computation) and the Earth's topography was used as the input signal (Earth2014.TBI2014.degree10800). The first figure shows the reference signal up to degree 10,800 and the second one reveals the accuracy achieved (recovered minus reference coefficients), both in terms of (difference) degree amplitudes (square root of degree variances). The third figure shows computation times required by the "SHA_GL.m" routine to perform the analysis. The tests were conducted on a PC with Intel® CoreTM i7-6800K CPU, 128 GB of RAM (no more than a few GBs were needed in this example though) and a 250 GB SSD drive.

When using the routines, please use the following reference:

Bucha, B., Hirt, C., Kuhn, M., 2019. Cap integration in spectral gravity forward modelling up to the full gravity tensor. Journal of Geodesy, https://doi.org/10.1007/s00190-019-01277-3.

Go to top

Molodensky's truncation coefficients for cap-modified spectral gravity forward modelling (near- and far-zone gravity effects implied by topographic masses)

Cap-modified spectral gravity forward modelling is a spectral technique to deliver near- and/or far-zone gravity effects implied by topographic masses, the shape of which is given as a surface spherical harmonic series. From the numerical point of view, one of its most difficult parts is the accurate computation of the associated Molodensky's truncation coefficients. This is because this evaluation may necessiate to extend the number of significant digits from 16 (double precision) to, say, 64 or even 256. Here, this was achieved via the ADVANPIX toolbox, which can be used in Matlab.

Some pre-computed ready-to-use values of the near-zone and far-zone truncation coefficients that were computed within the Bucha et al. study (see below) are available here (~2.5 GB). The coefficients are evaluated for

• the spherical distance of $\psi_0 = 100000\ \mathrm{m}\ /\ 6378137\ \mathrm{m}$ (100 km integration radius from the evaluation point),
• the reference sphere having the radius $R = 6378137\ \mathrm{m}$,
• the radius of the evaluation point $r = 6378137\ \mathrm{m} + 7000\ \mathrm{m}$,
• harmonic degrees $n = 0,\dots,21600$,
• topography powers $p = 1,\dots,30$,
• radial derivatives $k = 0,\dots,40$,
• the first- and second-order horizontal derivatives.

The coefficients were computed using 256 significant digits, ensuring 24-digit accuracy or better. After the evaluation, the coefficients were converted to double precision with 16 significant digits. Importantly, in some cases, the loss of significance errors may be encountered during the spherical harmonic synthesis when using the coefficients (see the reference below).

Please use the following reference when using either the codes or the coefficients:

Bucha, B., Hirt, C., Kuhn, M., 2019. Cap integration in spectral gravity forward modelling up to the full gravity tensor. Journal of Geodesy, https://doi.org/10.1007/s00190-019-01277-3.

Go to top

TView package — Single, total and cumulative viewshed computations from digital elevation models

Using an input digital elevation model given as a grid (DEM), the TView package allows to identify DEM cells that are visible from a user-defined point (single viewshed). Consider, for instance, a hilly area over some part of the Earth's surface. Based on the input DEM, the TView package can tell you which parts of the territory are visible from a particular observation point that can be located, e.g., at a summit or in a valley. The two figures below show a DEM and the computed visibility map (single viewshed).

Now imagine you want to find the sweet spot that offers the best view (that is, the location with the largest number of visible DEM cells). This can be done by evaluating the single viewshed for each cell of the DEM, followed by some simple summation scheme (see the reference below). In the figure below, the undulation shows terrain elevations over a small part of Slovakia, while the color says how much of the area can be seen from that particular location. The red color indicates an excelent view and the blue color means poor visibility. The visualization was kindly done by Alexandra Rášová in ArcScene.

The computation of total viewshed is generally a time consuming task. The TView package is therefore parallelized via the "parfor" loop, so that it can be used on HPC clusters to lower computational times. An example of the total viewshed can also be found here as a web map service. It offers a total viewshed grid over the entire Slovakia that was evaluated from SRTM v4.1 (1 arcsec spatial resolution, roughly equals to ~30 m) with an observation radius of 20 km for each cell. The grid was computed by Alexandra Rášová (Rášová, A., 2017. Substantive visibility analysis in GIS. PhD Thesis, Slovak University of Technology in Bratislava, Slovakia, 83 pp, in Slovak, Download — therein, click the "Prehliadať" button).

The TView package supports the flat-Earth, spherical and ellipsoidal approximations of the Earth's shape. Please use the following reference when using the TView package:

Rášová, A., Bucha, B. Large–scale parallel computation of total viewshed in MATLAB. Submitted to Computers and Geosciences.

Go to top

SRBF_package — Routines for gravity field modelling with spherical radial basis functions

SRBF_package is a collection of a few routines for gravity field modelling with band-limited spherical radial basis functions (SRBFs), allowing both analysis and synthesis of the gravity field. It is an enhanced version of the MATLAB-based routines used by Bucha et al. (2016). As an example of a band-limited SRBF, the figure below depicts the smoothed Shannon SRBF.

The package allows SRBF synthesis of the gravitational potential, the gravitational vector in LNOF and the gravitational tensor in LNOF. To the best of our knowledge, the equations to synthesize the gravitational tensor are novel and are more efficient than other formulae (e.g., that from Appendix of Bucha et al., 2016). This is because only three radial basis functions are needed to compute six unique elements of the gravitational tensor. This may significantly improve computational speed, as the evaluation of band-limited spherical radial basis functions (the sum over harmonic degrees) is the most time-consuming part of the synthesis.

SRBF_package is available for download here. Enclosed are also test data, a manual with a detailed description of the package and a script with two test computations to verify the correctness of the implementation. While the test examples are based on global gravity field modelling for clarity, SRBF_package was developed primarily for regional applications with a residual gravity field. This is because SRBF_package is significantly slower when compared to FFT-based SHS and SHA, both of which are designed for global applications. To improve the computational speed of the package, SRBF synthesis is parallelized over computational points via the parfor loop. Although the package was designed mainly for gravity field modelling, other applications are possible. Further details can be found in the attached PDF.

Please use the following reference when using SRBF_package:

Bucha, B., Janák, J., Papčo, J., Bezděk, A., 2016. High-resolution regional gravity field modelling in a mountainous area from terrestrial gravity data. Geophysical Journal International 207, 949-966, http://doi.org/10.1093/gji/ggw311.

Go to top

List of publications

• Bucha, B., Hirt, C., Yang, M., Kuhn, M., Rexer, M., 2019. Residual terrain modelling (RTM) in terms of the cap-modified spectral technique: RTM from a new perspective. Journal of Geodesy, https://doi.org/10.1007/s00190-019-01303-4, free full-text view-only version available

• Hirt, C., Bucha, B., Yang, M., Kuhn, M., 2019. A numerical study of residual terrain modelling (RTM) techniques and the harmonic correction using ultra-high-degree spectral gravity modelling. Journal of Geodesy 93, 1469-1486, https://doi.org/10.1007/s00190-019-01261-x, free full-text view-only version available

• Bucha, B., Hirt, C., Kuhn, M., 2019. Cap integration in spectral gravity forward modelling up to the full gravity tensor. Journal of Geodesy 93, 1707-1737, https://doi.org/10.1007/s00190-019-01277-3, free full-text view-only version available

• Hirt, C., Yang, M., Kuhn, M., Bucha, B., Kurzmann, A., Pail, R., 2019. SRTM2gravity: An Ultrahigh Resolution Global Model of Gravimetric Terrain Corrections. Geophysical Research Letters 46, 4618-4627, https://doi.org/10.1029/2019GL082521, free full-text view-only version available

• Bucha, B., Hirt, C., Kuhn, M., 2019. Divergence-free spherical harmonic gravity field modelling based on the Runge—Krarup theorem: a case study for the Moon. Journal of Geodesy 93, 489-513, https://doi.org/10.1007/s00190-018-1177-4, free full-text view-only version available

• Bucha, B., Hirt, C., Kuhn, M., 2019. Cap integration in spectral gravity forward modelling: near- and far-zone gravity effects via Molodensky's truncation coefficients. Journal of Geodesy 93, 65-83, https://doi.org/10.1007/s00190-018-1139-x, free full-text view-only version available

• Rexer, M., Hirt, C., Bucha, B., Holmes, S., 2018. Solution to the spectral filter problem of residual terrain modelling (RTM). Journal of Geodesy 92, 675-690, https://doi.org/10.1007/s00190-017-1086-y, free full-text view-only version available

• Bucha, B., Janák, J., Papčo, J., Bezděk, A., 2016. High-resolution regional gravity field modelling in a mountainous area from terrestrial gravity data. Geophysical Journal International 207, 949-966, http://doi.org/10.1093/gji/ggw311

• Bucha, B., Bezděk, A., Sebera, J., Janák, J., 2015. Global and regional gravity field determination from GOCE kinematic orbit by means of spherical radial basis functions. Surveys in Geophysics 36, 773-801, http://doi.org/10.1007/s10712-015-9344-0, free full-text view-only version available

• Bucha, B., Janák, J., 2014. A MATLAB-based graphical user interface program for computing functionals of the geopotential up to ultra-high degrees and orders: Efficient computation at irregular surfaces. Computers & Geosciences 66, 219-227, http://doi.org/10.1016/j.cageo.2014.02.005

• Bucha, B., Janák, J., 2013. A MATLAB-based graphical user interface program for computing functionals of the geopotential up to ultra-high degrees and orders. Computers & Geosciences 56, 186-196, http://doi.org/10.1016/j.cageo.2013.03.012

Preprints are freely available at my ResearchGate profile. For some of the papers, free full-text view-only versions are available (look for the text "full-text view-only" in the aforementioned references).

Go to top

Go to top