Segments in Mojo

After implementing points in Mojo, with some assistance from ChatGPT and Phind, we now turn our attention to segments, which represent the next fundamental geometric structure preceding polylines (or simply lines).

A segment consists of a start point and an end point, rendering it an oriented entity. It is implemented as a Mojo structure, that implements a few methods, such as returning the segment start and end points, calculating the total or horizontal lenght, and so on.

To facilitate point copying when returning a segment's start or end point (via the start_point and end_point methods), the __copyinit__ method was introduced to the Point structure: 

  fn __copyinit__(inout self, existing: Point): self.coords = existing.coords 

 This method allows for efficient copying of point coordinates, which is crucial when working with segments. 

Moving on to the Segment structure, apart from incorporating several basic methods (such as dx, midx, length, horizontal_length, and others), two significant methods were added for calculating the azimuth and plunge of a Segment instance. These concepts are extensively utilized in geological contexts. 

The azimuth refers to the angle, within the horizontal plane, between the horizontal projection of the segment and the Y-axis. This angle is measured clockwise, commencing from the Y-axis, spanning from 0° to 360°. 

The plunge represents the dip of the segment in the vertical plane. It signifies the vertical angle between the horizontal plane and the segment. Its range extends from -90° to +90°, where positive values denote a downward dip, and negative values signify an upward dip. 

Both azimuth and plunge are expressed as Float64 values, but their values may remain undefined under specific circumstances. The azimuth becomes undefined when the segment is vertical or possesses zero length (i.e., when the start and end points coincide). Similarly, the plunge remains undefined when the segment's length is zero. 

To address these edge cases, both the azimuth() and plunge() methods return an Optional[Float64], with Optional being imported from the standard library's collections module. If the value is invalid, None is returned; otherwise, the valid value is returned. 

Notably, the valid return can also be encapsulated into an Optional — both approaches are considered valid. When presenting the result, the value must be extracted using the or_else method, where the input to this method serves as a placeholder for missing data. 

In the example code, a conventional no-data value frequently employed in Geographic Information Systems (GIS) was utilized by defining an alias: 

   alias NULL_ORIENTATION = -999999999.99999 

 This alias is subsequently applied as follows:

   print("Segment plunge:", segment.plunge().or_else(NULL_ORIENTATION)) 

In this scenario, when the result is valid, it will be printed; conversely, if the result is invalid, the no-data value defined in NULL_ORIENTATION will be displayed. 

By implementing these features, the Segment structure becomes more robust and versatile, capable of handling various geometric calculations essential in geological applications and beyond. 

 

The code is available at: https://gitlab.com/mauroalberti/kira

 

An example of defining a point structure in Mojo

Mojo is a relatively new language (introduced in 2022) designed for parallel computing. It is compatible with SIMD (Single Instruction Multiple Data) architectures and is intentionally interoperable with Python. SIMD allows for the application of vectorized operations, resulting in significant speedups in processing.

Currently, Mojo is available only for Ubuntu and macOS.

In this example, we will define a point geometry and implement it using Mojo.

In this context, using SIMD to store point coordinates is advantageous because it enables the application of vectorized operations to those coordinates.

To define a SIMD variable, it’s necessary to specify both the data type and the number of stored values. Here is an example of initializing a SIMD variable:

var vec = SIMD[DType.int8, 4](1, 2, 3, 4)

In Mojo, you cannot use Float32 or Float64 directly to define the variable type. Instead, you need to use DType.float32 or DType.float64 for floating-point values.

Another requirement for SIMD variables is that the number of elements must be a power of two. Therefore, it is not possible to use SIMD[DType.float64, 3] to store just the x, y, and z coordinates of a point. To meet the minimum size requirement of 4 elements, you can add a time variable (t) of type Float64 to the point coordinates. This avoids wasting memory space.

The point coordinates are stored as follows:

var coords: SIMD[DType.float64, 4]

Similar to Python, a class-like structure in Mojo is initialized using the __init__ method. To indicate that the instance is mutable, the self parameter is marked as inout:

fn __init__(inout self, x: Float64, y: Float64, z: Float64, t: Float64 = 0.0)

Since the values are stored in a SIMD variable, the elements are accessed by their index:

fn x(self) -> Float64:
	return self.coords[0]
    

To calculate the distance between two points, you cannot use the sum() or sqrt() functions directly on a SIMD variable. Instead, you must extract the scalar values by index and perform the scalar operations manually.

The square root function is imported from the math module like this:

from math import sqrt

Here is the complete implementation of the point structure:

from math import sqrt

struct Point :

    var coords: SIMD[DType.float64, 4]

    fn __init__(inout self, x: Float64, y: Float64, z: Float64, t: Float64 = 0.0):

        self.coords = SIMD[DType.float64, 4](x, y, z, t)

    fn x(self) -> Float64:
        return self.coords[0]

    fn y(self) -> Float64:
        return self.coords[1]

    fn z(self) -> Float64:
        return self.coords[2]

    fn t(self) -> Float64:
        return self.coords[3]

    fn distance_to(self, other: Point) -> Float64:

        var delta = self.coords - other.coords
        var sum_of_squares = delta[0]*delta[0] + delta[1]*delta[1] + delta[2]*delta[2];
        return sqrt(sum_of_squares[0])
    
    

And here’s an example of how to use this structure in a main function:

fn main():

    var pt1 = Point(1.0, 2.0, 3.0, 1222.34)
    var pt2 = Point(2.0, 3.0, 3.0, 22.34)
    print(pt1.distance_to(pt2))
      

You can also view this example on GitLab:

Mojo Point Structure Example

geogst, a new Python module for Structural Geology

 

geogst is a new Python module for structural geology available on PyPi, making it easily installable via the classic command pip install geogst (or python -m pip install geogst).

This module was primarily developed to facilitate the creation of stereonets for geological data, determining the intersection between geological planes and topographic surfaces (expressed through Digital Elevation Models, DEMs), and generating the skeletons of geological profiles. Among its main features, the module allows for the calculation of topographic profiles, adding geological attitudes, and determining the intersections of geological traces with these profiles.

 

Additionally, geogst includes example geospatial datasets that can be used to explore its functionalities, such as creating profiles directly within Jupyter Notebooks.

Here are a few examples of Jupyter Notebooks for creating geological profiles:

• Geologic profiles map animations - the Timpa San Lorenzo structure (Southern Italy):
  Geologic profiles with animation - Timpa San Lorenzo

• Geologic profiles - the 1997 seismic sequence in Colfiorito (Central Italy):
  Seismological profiles - Colfiorito

This module will form the foundation for the modules in the qgSurf plugin for QGIS. Currently, qgSurf directly includes geogst as a submodule, so no separate installation is required. In future versions of qgSurf (initially experimental), the geogst module will be automatically installed if it is not already present.

Since QGIS does not allow the upload of compiled code in Python modules, a key advantage of using geogst as a separate module in qgSurf is the potential to incorporate compiled code in Fortran, C++, and Rust, significantly improving the speed and efficiency of these tools.

Conclusion and call to action: If you are a geologist or a developer working with geological data, we invite you to explore geogst and contribute to its development. Your experiences and feedback are crucial to improving and growing the community that uses ggSurf.

Mojo

A new AI-oriented language is on the scene: Mojo, which aims to become a superset of Python. It follows Python semantics but, thanks to recent and advanced compilation techniques, might achieve speedups as large as 100,000 times or more compared to standard Python.

Mojo was created by Chris Lattner, the creator of LLVM, Clang, Swift, Swift for TensorFlow and MLIR.

Currently, Mojo is only supported on Linux and macOS.

The source code is hosted at GitHub, from where you can clone the repository locally. Within the downloaded repository, the 'examples/notebooks' directory contains Jupyter Notebooks that you can run using Mojo (provided that Mojo and the Jupyter plugin are installed).

A very brief introduction to the Mojo language is provided in HelloMojo.ipynb. Another notebook, Matmul.ipynb, implements a matrix multiplication example in both Python and Mojo. The example showcases how optimization in Mojo can lead to a remarkable speedup, reportedly around 455,127x compared to Python. It's important to note that achieving such optimized Mojo code requires considerable effort and a solid understanding of Mojo.

That said, I tried running the code in that notebook on my old HP laptop, which lacks a GPU, and has the following specifications:

inxi -Fxxxzr

System:
  Kernel: 5.15.0-118-generic x86_64 bits: 64 compiler: gcc v: 11.4.0
  Desktop: Xfce 4.18.1 tk: Gtk 3.24.33 info: xfce4-panel wm: xfwm 4.18.0
  vt: 7 dm: LightDM 1.30.0 Distro: Linux Mint 21.3 Virginia
  base: Ubuntu 22.04 jammy
Machine:
  Type: Laptop System: HP product: HP Laptop 15-bs0xx v: Type1ProductConfigId
  serial: <superuser required> Chassis: type: 10 serial: <superuser required>
  Mobo: HP model: 832B v: 23.37 serial: <superuser required> UEFI: Insyde
  v: F.21 date: 07/04/2017
CPU:
  Info: dual core model: Intel Core i5-7200U bits: 64 type: MT MCP
  smt: enabled arch: Amber/Kaby Lake note: check rev: 9 cache: L1: 128 KiB
  L2: 512 KiB L3: 3 MiB
  Speed (MHz): avg: 1134 high: 1245 min/max: 400/3100 cores: 1: 1040
  2: 1097 3: 1245 4: 1157 bogomips: 21599

Even though the speedup I achieved was an order of magnitude lower than the original 455,127x reported in the notebook, it was still impressive: 10,032x! It’s likely that using a newer machine, possibly equipped with a GPU, I could achieve speedups on the order of 100,000x.

As reported, tools like Cython or Numba achieve speedups typically in the range of 10-100x. 

Moreover, the Mojo code ran without any issue on my old laptop running Linux Mint.

So Mojo is a very interesting language, even in its infancy.

Other AI-oriented languages to explore, as described in the insightful post by James Thomason in VentureBeat, "Mojo Rising: The resurgence of AI-first programming languages" include Bend and JAX.

 

Note: the text was checked in ChatGPT.

 

 

SAVI index and agricultural landscapes

SAVI is the acronym for Soil Adjusted Vegetation Index that, as the name suggests, it is one of the many indices for the calculation of the vegetation index.

The general formula for the SAVI index is (indexdatabase.de):

     800 nm - 670 nm
  ________________________  ( 1 + L)

   800 nm + 670 nm + L

   
where L is between -0.9 and 1.6.

In the case of Sentinel 2 data, SAVI uses two bands, B8 and B4, that have the maximum available resolution for Sentinel 2 images, i.e., 10 meters, so the results have a quite high spatial resolution.

B8 is in the NIR (842 nm) while B4 registers the red-NIR transition (665 nm).

For Sentinel 2, the index is therefore calculated as:

     B8 - B4
  _______________   ( 1 + L)

   B8 + B4 + L
   
The default value for L in indexdatabase is 0.5, and in fact testing for a few values (-0.9, -0.5, 0.0, 0.5, 1.0, 1.5) in the prescribed range, the most contrasted results are obtained for L = 0.5.

One information that is obtained from applying such an index to intensive farming landscapes is that it allows to delineate well the agricultural field limits, to compare the texture and the "spatial styles" of farming across different areas, for instances separated by rivers that may have acted as administrative or political borders.

In the following images, I present some examples of textures and styles of limits between fields for different areas in the Piedmont region, derived from a Sentinel 2B image, acquired on 2023/04/04 at 10:25. The SAVI index was calculated using the SNAP software by ESA.

The used L ("soil brightness correction factor") value is 0.5, chosen after some experimentation.

The figures below are generated from within SNAP, using a 8-4-3 composite for infrared view (upper window), while the SAVI index band is represented in the lower window. In the 843 composite, vegetation is red while rivers and roads are black.

Strong agricultural landscapes contrast between the western side (left) and eastern  side (right) of the Sesia river (central, with a North-South orientation). The western landscape is characterized by smaller sizes than in the eastern part, and is also more regular as geometries. Possibly it is the result of different political systems between the two sides in the previous centuries.

 Strong East-West differences in field sizes, orientations and textures without a clear separating physical boundary (such as the Sesia river in the previous example). It could be related to different grown agricultural products between the two zones (?rice to the East), possibly also related to distinct administrative/political situations.  

Irregular field borders probably due to the limiting rivers and their location changes with time.

 

 

Rayleigh Correction for Sentinel data using SNAP

Removing atmospheric effects from satellite images is tricky, since the theoretical foundations are not completed understood, it is difficult to model the physical conditions, and also because we generally lack many information regarding the atmospheric condition at the time and location of satellite image acquisition.

One procedure that attempts to partially remove atmospheric effects is implemented in the SNAP software application, that is released by ESA and mainly devoted to the processing of Sentinel data, both optical and radar.

This method is named "Rayleigh Correction" and attempts to mitigate gaseous and Rayleigh effects. It does however do not consider aerosol effects.

In the screenshot below you see how to open the method window.

The GUI allows to select the bands to be corrected (B2-4 in the example below), besides other parameters (e.g., sea-level pressure). For this example the output resolution was changed to 10 meters, to match that of input bands (B2-4). The other parameters were left unchanged.

The two windows below display a an example of Sentinel 2 data (S2A MSI, 2017-05-27) pre- and post- Rayleigh correction. The displayed zone is comprised between Chivasso to the West and Verolengo to the East, and the Po River to the South and the Canale Cavour to the North (Piedmont, Italy). 

In the upper window the original image is displayed as a 4-3-2 composite. In the lower window the same zone is represented using a 4-3-2 composite of Rayleigh-corrected bands.

The original bands are darker and with a low contrast, while the corrected bands present a larger chromatic variance. Note that the lightest pixels in the original band tend to saturate to white, somewhat loosing information.

Water, vegetation and urban constructions are more vivid and looks like more "natural", with more chromatic variations than in the original bands.

Quick calculation of band indices from satellite images with GDAL

Band indices of satellite images, such as the NDVI, can be routinely calculated from within applications, such as ENVI, SNAP and so on.

Sometimes it is preferable to calculate band indices using a programming language: it could be faster than using applications, we could not dispose of the required applications or we need to automate a set of processing by using a programming language.

Having GDAL installed, it is quick and easy to calculate band indices (and a lot more) directly from the prompt/shell. I was drawn to this useful solution by the answer provided by

The used command is gdal_calc.py and we can specify one, two or more input bands to perform calculations, for instance, for calculating the Normalized Difference Water Index (NDWI) from a spatial subset of a Sentinel 2 image saved from SNAP as BEAM-Dimap, we use:

gdal_calc.py -A B3.img -B B8.img --outfile=ndwi.tif --calc="((B-A)/(B+A))"

Obviously the shell working directory is the folder containing the source image files. Morevover we do not have to define the used band within the input .img file since there is only one band.
The complete description of the command is available at the page 'gdal_calc.py'.


An example of calculation of the Normalized Difference Water Index is here presented. Three Sentinel 2 SLC image of the Po river in the Piedmont area between Chivasso and Casale Monferrato (Northern Italy) were spatially subset in the range lat 45.1 -> 45.2 and long 7.8 -> 8.6 from within SNAP (Fig. 1).

 Fig. 1. Map of the studied area in the central Piedmont region (Northern Italy). Basemap: Google Earth. Created with QGIS.

The NDWI was calculated with the command line as listed above:

gdal_calc.py -A B3.img -B B8.img --outfile=ndwi.tif --calc="((B-A)/(B+A))" 

From the analysis of the histogram of the 2017 image, a value of 10 was chosen of a threshold for masking the non-water pixels:

gdal_calc.py -A ndwi.tif --outfile=ndwi_masked.tif --calc="A>10"

The result is represented in Fig. 2A-C, created with QGIS. The result quality is more than acceptable (for the 2017 image the spotted fields to the North of the Po river could possibly correspond to rice fields, but it has to be better investigated).

A - 2017-05-27

B - 2021-03-30

C- 2023-04-04

 Fig. 2. NDWI > 10 for 2017 (a), 2021 (b) and 2023 (c). Basemap: Google Earth. Created with QGIS.

By comparing the results at different times, it is possible to study the temporal evolution of the rivers (in this case, the Po river, Orco torrent and Canale Cavour near the city of Chivasso (Fig. 3).

 Fig. 2. NDWI > 10 for 2017 (red), 2021 (green) and 2023 (blue) in correspondence of the city of Chivasso. The water bodies are the Po river, the Orco torrent and the artificial Canale Cavour. Basemap: Google Earth. Created with QGIS.

 
If we want to automate the task, when we have for instance a lot of input images or we have to perform routinely the same task on new datasets, the run command can be created on-the-fly and executed from within a script created with languages such as Go, Rust, Ocaml and so on, perhaps even concurrently.

An example of atmospheric effect in differential interferograms in the Piedmont Po Plain area (Vercelli-Trino, Northern Italy)?

With differential interferometry from Interferometric Synthetic Aperture Radar (InSAR) it is possible to monitor subsidence or uplift along line-of-sight. One popular satellite constellation for interferometry is Sentinel-1, operated by ESA. It is composed by two satellites, Sentinel-1A and Sentinel-1B, carrying a C band SAR sensor.

Sentinel data are freely available from the Copernicus Open Access Hub (https://scihub.copernicus.eu/).

I downloaded a few Sentinel-1 images of the Piedmont Po Plain (Vercelli-Trino Vercellese and surroundings) to try deriving some interferometric results to interpret.

Fig. 1: map of the studied are in the Piedmont Po Plain (Northern Italy). Google satellite basemap. Created with QGIS.

To see the temporal evolution of the results, a set of four images for the following 2020 periods were used:

• 05/01/2020: S1B_IW_SLC__1SDV_20200501T171426_20200501T171453_021391_0289B0_577E
  
• 05/13/2020: S1B_IW_SLC__1SDV_20200513T171427_20200513T171454_021566_028F1C_D01D
  
• 05/25/2020: S1B_IW_SLC__1SDV_20200525T171428_20200525T171455_021741_029439_E72E
  
• 06/06/2020: S1B_IW_SLC__1SDV_20200606T171428_20200606T171455_021916_02997C_907D
  

Three differential interferometric results were obtained, using SNAP software by ESA and following the indications presented in tutorials and videos published by ESA. The calculated LOS displacements were masked to hide areas were results have low coherence.

The displacements are presented in the following Fig. 2A-C, where Fig. 2A refers to the 05/01->13 time interval, Fig. 2B to the 05/13 -> 25 interval and  Fig. 2C refers to the final 05/25->06/06 time interval.

 (A)

(B)

(C)

Fig. 2: masked displacements along Line-Of-Sight for the three time intervals.  A) 05/01->13 time interval; B) 05/13->25 time interval; C) 05/25->06/06 time interval. The central E-W line is the profile trace used to derive Fig. 3. Google satellite basemap. Created with QGIS.


We can see that in the first time interval (05/01->13) the regional displacements are mainly negative, in the second (05/13->25) positive with a clear E-W gradient, while the situation is more nuanced in the third one (05/13->06/06), without a clear prevalence of positive or negative values.

A profile of the three displacement rasters, along a E-W direction (red line in Fig. 2), was created using the 'gst' Python module (https://gitlab.com/alberese/gst), as described in the 'Trino' Jupyter notebook (https://gitlab.com/alberese/gst/-/blob/master/docs/interferometry/Trino.ipynb).

Fig. 3: Profiles of the masked displacement rasters along the trace in Figs 2A-C. The first period (I) is in red, the second (II) in green and the third (III) in blue. Created with Python modules matplotlib and gst.

The profiles depict the same situation as visible in map, and evidence, near to rightmost profile end, in correspondence of Vigevano town, a local maximum for the first period while the second period presents a local minimum and in the third period there is a flat (Fig. 3).  

It is difficult to consider these variations, both regional and local (Figs. 2 and 3) as related to topographic surface movements, since no constant trend is observed in the three periods, but the displacements during the second period somewhat reverse those of the first periods. The observed results are therefore interpreted as the results of atmospheric effects, particularly strong in the 05/13 image, that impacted with opposite results the displacements of the first (05/01->13) and second period (05/13->25). The results of the third period are more similar to those expected for a no-movement differential interferogram and possibly reflect only local atmospheric effects for both the 05/25/2020 and 06/06/2020 images.   

Calculating lithological indices for Niger ASTER data with the QGIS plugin 'res'

From ASTER satellite data it is possible to derive indices related to surficial mineralogical assemblages at the Earth surface (Gupta, 2017 and references within). Since ASTER cell size is 30 m in SWIR and 90 m in TIR (Gupta, 2017), the provided spatial resolution is just medium, not high.

A preliminary version of 'res', a still in-progress QGIS plugin, is here presented. It allows to automatically calculate a set of ASTER ratio indices, as discussed in Gupta (2017, cf. Table 19.3), that are useful for mineralogical identification:

1. Ferric Iron Index (Rowan and Mars 2003): higher values -> greater amount of ferric iron

2. Ferrous iron index (Rowan et al., 2005)

3. Gossan index (Velosky et al., 2003): indicates the presence of highly ferruginous rocks, deriving from the alteration of sulfides (https://www.sciencedirect.com/topics/earth-and-planetary-sciences/gossan)

4. Generalized hydroxyl alteration

5. Propylitic Index (Rowan et al., 2005) -> mafic-ultramafic suite of rocks

6. Phyllic index. Three versions presented: a) Rowan and Mars (2003); b) Rowan et al. (2005); c) Ninomiya (2003a)

7. Argillic index. Two versions: Mars and Rowan (2006), Ninomiya (2003a)

8. Calcite index. Two versions: Rowan and Mars (2003), Ninomiya (2003b)

9. Quartz index. Two versions: Ninomiya et al. (2005), Rockwell and Hofstra (2008)

10. Silica index. Two versions: Ninomiya and Fu (2002), Rowan et al. (2006)

11. Carbonate index. Ninomiya and Fu (2002)

12. Mafic index and MI corrected. Ninomiya et al. (2005), Ninomiya et al. (2005)

13. Ultramafic index. Rowan et al. (2005)

Case study: data from Niger

A very basic example from Niger is presented here. Arid deserts are an ideal setting for superficial geology investigations, thanks to the aridity and the limited presence of vegetation.

Unfortunately I do not found a reliable lithological surface map, to be used to ground-control the results. If anyone is able to point me to a freely-available lithological/mineralogical map of this zone I will be grateful.

This example has no atmospheric correction and there are some clouds in the scene. The calculated indices are distorted by image features such as clouds.

The first step is the definition of the input data source, by choosing the directory where the ASTER bands are stored (Fig. 1). The available band list is automatically filled.

Fig. 1. Input definition window.

After defining the input data sources, we define the indices to be calculated (Fig. 2).

Fig. 2. Choice of indices to calculate.

We store the calculated indices, plus a natural color composite image, into a specific folder (Fig. 3).

Fig. 3. Definition of output folder. The log window is not yet populated.

Note that the log form currently does not display anything.

Now we view some results within QGIS, using as a basemap the Google Earth composite (Fig. 4). As previously said, the area is desertic, with very limited vegetation.

Fig. 4. The studied zone in Niger (Africa). Scene borders in white.

 The composite 3N-2-1 from Aster data is automatically computed by the plug-in (Fig. 5).

Fig. 5. ASTER composite 3N-2-1 with stretched pixel values.

ASTER ratio indices

Clouds are masked in the images below. All indices maps with useful results are presented applying a 2.0 standard deviation stretch.

Some bands maps are more speckled and/or with stripes (e.g., for the mafic- and ultramafic indices, Figs. 11-13, 16), due to noise possibly amplified in the index calculations, or with better spatial resolution (e.g., the gossan index, Fig. 8), due to the greater spatial resolution of involved original bands.

In general, all different versions of a specific mineralogical index produce results that are consistent with each other.

Iron-related indices

Both iron indices suggest a higher iron content in the northern region, particularly in the most northern parts (Figs. 6, 7).

 Ferrous iron index map

Fig. 6. Ferrous iron index map.

Ferric iron index map

Fig. 7. Ferric iron index map.

 

Gossan index

High Gossan index values, that suggest the presence of highly ferruginous rocks deriving from sulfide minerals alteration, are clearly correlated with red/dark red zones in the composite image (cf. Fig. 8 with 5). Generally high gossan index values are distributed in elongated strings, possibly corresponding to topographic slopes. They are mainly located in the eastern sector, particularly in its central part.

Fig. 8. Gossan index map.

Silica indices

 

The silica index is generally low in the scene region, apart from a major NW-SE elongated positive anomaly in the northern extreme of the scene, and a diffuse positive anomaly in the southermost part of the scene (Figs. 9, 10).

The NW-SE structure will also be notable in the mafic and ultramafic indices bands, with a reversed role (i.e., a negative anomaly, cf. Figs. 11-13 and 16, 17).

Silica index - Ninomiya and Fu (2002)

 

Fig. 9. Silica index (Ninomiya and Fu, 2002) map.

 Silica index - Rowan et al. (2006) 

Fig. 10. Silica index (Rowan et al., 2006) map.

Mafic indices

 The mafic indices mirror the silica indices: high silica indices zones correspond to low mafic indices zone and viceversa (Figs. 11, 12).

The northernmost NW-SE anomaly is negative, so that the outcrops could be high in silica and low in mafic minerals.  

Mafic index - Ninomiya et al. (2005)

Fig. 11. Mafic index (Ninomiya et al., 2005) map.

Corrected mafic index - Ninomiya et al. (2005)

 

Fig. 12. Corrected mafic index (Ninomiya et al., 2005) map.

Ultramafic index

The ultramafic index map presents the previously described negative anomaly in the northernmost part, as well as less intense negative anomalies in the southern part of the scene (Fig. 13).

Based on the silica, mafic and ultramafic indices, and also on the morphology in the composite ASTER image (Fig. 5), this northern structure could correspond to an outcrop of quartz sandstones, eventually forming a SE-dipping antiform, given the fan-shaped form of the anomaly.

Fig. 13. Ultramafic index map.

Filtering of results

 

The results, particularly for some indices (for instance, the ultramafic index) can be speckled. To help improve the situation, a filter tool, in particular, a median filter, has been added that can be applied to computed bands.

 

The input and output folders are defined in the 'Input/Output' tab of the 'Image Filters' plug-in sub-module (Fig. 14).

  

Fig. 14. Filtering I/O definition window.

 

 

The kernel size is defined in the 'Filters' tab (Fig. 15).

 

 

Fig. 15. Median filtering definition window.

In Figs. 16 and 17 we compare the original ultramafic index band with the filtered one: the resulting band is less speckled, even if the spatial resolution obviously do not improve.

 

 

Fig. 16. Northernmost part of the scene - original ultramafic index map.

Fig. 17. Northernmost part of the scene - ultramafic index map filtered with a 5x5 kernel.

Plug-in status

 

The 'res' plug-in is still in preparations, even if it should produce valid results in the majority of the cases. In particular, the help is still to be written. Additionally, not all calculated index bands are useful, since a few index bands consist of constant grid values. This problem has to be better investigated.

The plugin current version is 0.2. It can be downloaded as a zip file from:

https://gitlab.com/mauroalberti/res/-/releases/v0.2

 

Afterwards the zip file can be loaded into QGIS from the QGIS plugin loader.

References

Mars J.C., Rowan L.C., 2006. Regional mapping of phyllic- and argillic-altered rocks in the Zagros magmatic arc, Iran, using Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) data and logical operator algorithms. Geosphere 2 (3), 161–186.

Ninomiya Y., 2003. Rock type mapping with indices defined for multispectral thermal infrared ASTER data: case studies. Proc SPIE 4886, 123–132.

Ninomiya Y., 2003. A stabilized vegetation index and several mineralogic indices defined for ASTER VNIR and SWIR data. In: Proceedings international geoscience and remote sensing symposium (IGARSS-2003) (IEEE) 1294172. doi:10.1109/ IGARSS.2003.1294172.

Ninomiya Y., Fu B., 2002. Mapping quartz, carbonate minerals and mafic–ultramafic rocks using remotely sensed multispectral thermal infrared ASTER data. Proc SPIE 4710, 191–202.

Ninomiya Y., Fu B., Cudahy T.J., 2005. Detecting lithology with advanced spaceborne thermal emission and reflection radiometer (ASTER) multispectral thermal infrared “radiance-at-sensor” data. Remote Sens Environ 99, 127–139.

Gupta, R.P., 2017. Remote Sensing Geology. Springer, Third Edition.

Rockwell B.W., Hofstra A.H., 2008. Identification of quartz and carbonate minerals across northern Nevada using ASTER thermal infrared emissivity data, implications for geologic mapping and mineral resource investigations in well-studied and frontier areas. Geosphere 4(1), 218–246.

Rowan L.C., Mars J.C., 2003. Lithologic mapping in the Mountain Pass, California area using advanced spaceborne thermal emission and reflection radiometer (ASTER) data. Remote Sens. Environ. 84, 350–366.

Rowan L.C., Mars J.C., Simpson C.J., 2005. Lithologic mapping of the Mordor N.T., Australia ultramafic complex by using the advanced spaceborne thermal emission and reflection radiometer (ASTER). Remote Sens. Environ 99, 105–126.

Rowan L.C., Schmidt R.G., Mars J.C., 2006. Distribution of hydrothermally altered rocks in the Reko Diq, Pakistan, mineralized area based on spectral analysis of ASTER data. Remote Sens. Environ. 104, 74–87.

Velosky J.C., Stern R.J., Johnson P.R., 2003. Geological control of massive sulfide mineralization in the Neoproterozoic Wadi Bidah shear zone, southwestern Saudi Arabia, inferences from orbital remote sensing and field studies. Precambr. Res. 123 (2–4), 235–247.

SVD-based inversion of geological surfaces from topographic traces: preliminary assessment of reliability via eigenvalue log ratios

Topographic traces of geological surfaces store information about the surface attitudes. When considering the local scale, in the majority of the cases the surface can be approximated by a geometric plane.

Singular Value Decomposition (SVD) is one of the available methods for inverting the topographic traces, i.e., deriving the local, geological surface planar attitude.

SVD input and output matrices

Input data for the SVD methodology are expressed as an m x n matrix, here named A. When successful, the SVD method produces 3 matrices, named U, S and V* (* means transposed), so that A = U x S x V*. In the general case, U has dimensions m x m, S is an m x n diagonal matrix (i.e., all values are zero except for those on the main diagonal), while V* is an n x n matrix.

Input and output data for topographic point inversion

For topographic points inversions, m, i.e.. the number of matrix rows, is equal to the number of input points while n, the number of matrix columns, is equal to 3, since we are considering a 3D Cartesian space. The U matrix dimensions, i.e., m x m, are equal to the number of input points. S is an m x 3 diagonal matrix, whose main diagonal values stores the eigenvalues in descending magnitude order (s1 >= s2 >= s3).
V* is a 3 x 3 matrix that stores the eigenvectors of the A matrix.
In particular, the last row of the V* matrix, that is associated with the minimum value eigenvalue s3, contains the Cartesian components of the vector that is normal to the potential best-fit-plane.

For our purpose, only V* and S are needed for getting the estimated best-fit-plane (by using the last V* matrix row) and evaluating its reliability (by considering S matrix diagonal values).

What about the results reliability?

Some questions are: how can we be sure that the input points can be approximated by a plane? How can we evaluate how good is the fit between the computed result and the real plane? Can a quantitative metric of the result reliability therefore be defined?

To try to answer these questions, at least in a preliminary way, I've created a few simulations of points with different spatial distributions, in order to invert them and check the degree of concordance between the computed plane-normal eigenvector and the theoretical surface attitude. The eigenvalues has been used to try to infer the degree of result reliability.

To better characterize the eigenvalues, their negative log ratios were calculated:

nlr21 = - log (s2 / s1)
nlr31 = - log (s3 / s1)
nlr32 = - log (s3 / s2)

The negative signs before the log operators are used to avoid dealing with negative values. Since the numerator is always equal or lower than the denominator, the resulting negative log ratios are between 0 (for numerator equal to denominator) to infinity (for numerator equal to zero). Nan results can arise when both numerator and denominator eigenvalues are zero.

Software used for simulations and inversions

Both the simulations and the inversions were performed using the programs in the geoSurfDEM repository (currently available only in the dev branch). The SVD inversions are mainly dealt with using FORTRAN and linked Blas/Lapack libraries. C++ code is wrapped around the inversion FORTRAN routines and takes care of data input and output for both the simulations and the inversions.

Basic simulations

To start with, two very simple cases were considered, in order to get a feeling for the resulting eigenvalue distributions.

For these simple distributions, the specific points were directly written into the input files, as rows of x-y-z values (files in 'geoSurfDEM/test_data/InvertPoints' folder, named 'points_01.csv' to 'points_04.csv'). The different cases were then inverted with the code in the 'geoSurfDEM/InvertPoints' folder, specifically the compiled 'InvertPoints.cpp' program).

Co-planar points

These simulations consider perfectly co-planar points (just 4 for each plane dip subcase), lying in a horizontal (a), vertical (b) or 45°-dipping (c) plane.

Results 

The third eigenvector in the resulting V* matrix is always perfectly normal to the simulated plane. The s1 and s2 eigenvalues are equal or of the same order of magnitude (e.g., s1 = 2 and s2 = 1.41), while the s3 eigenvalue is always zero.
The nlr21 values are between 0 and 0.15 for all three cases (i.e., s1 and s2 equal or almost equal). Both nlr31 and nlr32 indices are infinity (since s3 is zero).
So for co-planar points we expect s1 = s2 >> s3, with nlr12 values around zero and very high values of both nlr31 and nlr32.

Collinear points

Three perfectly collinear, horizontal points were used as input data. Obviously no meaningful best-fit-plane exists.

Results

s1 is 1.4142 while s2 and s3 are both zero. So both nlr21 and nlr31 are infinity while nlr32 is Nan. 

Collinear points produces eigenvalues where s1 >> s2 = s3 and very large values of nlr21 and nlr31. The nlr32 index should be around 0 (apart when s2 and s3 are both zero).

Collinear and co-planar log ratios summary

The current results for the co-planar and collinear cases can be summarized as in the following diagram:

• co-planar data have very high nlr31 and nlr32 indices. 
• collinear points have very high nlr31 index and near zero nlr32 index.

 

Going uniform

When points have no preferred orientation, what are the resulting eigenvalues?

To answer this question, I've created uniform spherical distributions of point. I used a slightly modified C++ code published by Cory Simon in http://corysimon.github.io/articles/uniformdistn-on-sphere/. As before, the code is available in the geoSurfDEM repository (currently only in the dev branch and I have to check that code again since I've made modifications to the FORTRAN inversion routine...).

The simulated data to invert were created as such: 1000 random simulations were generated for 6 set of generated point number: 3, 4, 5, 10, 100 and 1000 points.
The files storing 3 points correspond in general to a unique well-defined inverted plane, apart from the cases where the points tend to collinearity (the case of three coincident points was excluded in the code). 

In the Addendum, all statistical plots for these simulations are presented. Here I show only the most significant plot, that complements the previous results.

When considering the nlr31 vs. nlr32, we see that perfectly co-planar data (n = 3) have high values of both these indices.

Spatially uniform data (simulations with n >> 3) have almost equal eigenvalues, so in the graph they are clustered towards the axes origin.

 

Conclusions

We can add together all the previous results in this summary sketch:

Where is the exact limit between almost co-planar and random points? And between collinear and co-planar points (undulating blue line)?

Unfortunately these simulations do not give us a quantitative answer. They say us when it is surely random, when surely co-planar or collinear, but further analyses where a quantitative measure of adherence of the empirical data to the theoretical (i.e., planar) model are needed. 

However, prior to applying other models, a check that our data are not random or collinear could help us to avoid producing false-positive results.

Addendum: Statistical properties of spherically uniform random points.

These are all the statistical plots for the spherical simulations.

The plots were generated with R and ggplot2. 

Simulations are subdivided by number of simulated points (from 3 to 1000).

The plots are not commented.

 

 

 

 

To invert the attitude of a geological trace

In 2016 I wrote a small utility in C++ and Fortran, geoSurfDEM, that calculates the intersections of a geological surfaces with a topographic surface (in the module "IntersectDEM", see post "geoSurfDEM: a C++ console application for determining intersections between 3D geological surfaces and topography") and that inverts the attitudes of geological traces (module "BestFitGeoplanes", see post "A Linux tool for calculating local best-fit plane attitudes from geological traces"" for the original description).

Provided a set of points of a geological trace, defined by their coordinates x-y-z, the BestFitGeoplanes module derives the local best-fit-planes, using a mathematical technique known as Singular Value Decomposition (SVD).
The used Fortran Lapack subroutine derives the eigenvectors and eigenvalues of the set of the source points expressed as a matrix (m x 3, where m is the number of points and 3 is the dimension of the embedding space, i.e., a 3D Cartesian space).

The eigenvector normal to the plane with the minimum data variance, i.e., the plane that better approximates the data assuming that they approximately subplanar, is the one in which we are interested.

The question is, how reliable is the result? If my data are collinear, no unique solution exists. The same applies when all the points are approximately coincident.

In the original version, in order to characterize the situations related to spurious results, the module calculated the spatial range of the source data along the three reference axes and also the volume of the bounding box (axes-aligned).

To verify whether they can really provided indications of errors (i.e., large differences between expected and calculated values),  a synthetic data set, made up by two geological traces, was inverted and the results analyzed.
The two synthetic geological traces derive from the analysis of the Timpa di San Lorenzo fault structure (Calabria, Southern Italy), where the northern and the southern fault structure segments attitudes can be approximated by two planes, respectively oriented 067.3°/39° and 077.3°/40° (with respect to the UTM 32632 CRS up direction)(see previous post "Along-trace profiles of the Timpa San Lorenzo fault structure (Calabria, Italy)"). The used DEM is a subset of the TINItaly DEM.

 

The intersection of two theoretical attitudes with the DEM were calculated using the 'DEM-plane intersections' of the qgSurf QGIS plugin, with the results saved as two point shapefiles. Points outside the fault spatial range were removed, the two layers were then merged into a shapefile and the data were exported as an x-y-z csv file, to be used as input for the 'BestFitGeoplanes' module.

The inversion results reproduce the theoretical attitudes very well (< 1° of difference), except for 4 values out of the total 311 (1.3 %). The spatial ranges of the source point data for each result do not appear however to be predictors of the result reliability.

To find a better quality indicator, a new code version outputs the values of the three eigenvalues for each result.
Generally the first eigenvalue (S1) has a value in the 10-50 range, the second (S2) is 1 or less, and the third eigenvalue (S3) is very low or zero.

In the analyzed data set, the "erroneous" results are characterized by low S2 values, around e-5. By considering the logs of the ratio between S2 and S1, the results with values lower than -5 are from grossly (>> 1° misfit) to mildly (1° or less) disoriented (Fig. 3). 

Indices comprised between -3.2 and -5.0 are correlated with mildly disoriented results (< 1° disorientation).

So using a threshold of -3.0 for this data set filters out all the low-quality.

Obviously, to check whether this empirical threshold can be considered general requires to consider other test cases. Possibly a collinearity index, to be defined, could be used an additional or alternative quality indicator.

References

 

geoSurfDEM - https://gitlab.com/mauroalberti/geoSurfDEM