Utilities to process and compute parameters. Module for Algebra calculus.

watex.utils.exmath.adaptive_moving_average(data, /, window_size_factor=0.1)

Adaptative moving average as smoothing technique.

Parameters:

• data (Arraylike) – Noise data for smoothing

• window_size_factor (float, default=0.1) – Parameter to control the adaptiveness of the moving average.

Returns:

result – Smoothed data

Return type:

Arraylike

Example

>>> import matplotlib.pyplot as plt
>>> from watex.utils.exmath import adaptive_moving_average
>>> # Sample magnetotelluric data (replace this with your own data)
>>> # Example data: a sine wave with noise
>>> time = np.linspace(0, 10, 1000)  # Replace with your actual time values
>>> mt_data = np.sin(2 * np.pi * 1 * time) + 0.2 * np.random.randn(1000)  # Example data
>>> # Function to calculate the adaptive moving average
>>> # Define the window size factor (adjust as needed)
>>> window_size_factor = 0.1  # Adjust this value based on your data characteristics
>>> # Apply adaptive moving average to the magnetotelluric data
>>> smoothed_data = adaptive_moving_average(mt_data, window_size_factor)
>>> # Plot the original and smoothed data
>>> plt.figure(figsize=(10, 6))
>>> plt.plot(time, mt_data, 'b-', label='Original Data')
>>> plt.plot(time, smoothed_data, 'r-', label='Smoothed Data (AMA)')
>>> plt.xlabel('Time')
>>> plt.ylabel('Amplitude')
>>> plt.title('Adaptive Moving Average (AMA) Smoothing')
>>> plt.legend()
>>> plt.grid(True)
>>> plt.show()

watex.utils.exmath.betaj(xj, L, W, **kws)

Weight factor function for convoluting at station/site j.

The function deals with the discrete hanning window based on ideas presented in Torres-Verdin and Bostick, 1992, https://doi.org/10.1190/1.2400625.

Parameters:

• xj – int, position of the point to compute its weight.

• W – int, window size, presumes to be the number of dipole.

• L – int : length of dipole in meters

• kws – dict , additional scipy.intergate.quad() functions.

Returns:

Weight value at the position xj, prefix-x`is used to specify the direction. Commonly the survey direction is considered as `x.

Example:

>>> from watex.exmath import betaj
>>> # compute the weight point for window-size = 5 at position j =2
>>> L= 1 ; W=5
>>> betaj (xj = 2 , L=L, W=W )
... 0.35136534572813144

watex.utils.exmath.bnn_inversion(AB, rhoa, rho0, h0, **kwargs)

Perform VES data inversion using a Bayesian Neural Network (BNN).

Parameters:

• AB (np.ndarray) – Array containing the spacing of current electrodes.

• rhoa (np.ndarray) – Array containing apparent resistivity values.

• rho0 (float) – Initial guess for resistivity.

• h0 (float) – Initial guess for thickness of the first layer.

Returns:

Inversion results including estimated resistivity and thickness.

Return type:

np.ndarray

watex.utils.exmath.butterworth_filter(data, /, freqs, fs=None, frange=None, order=5, plot=False)

Defines a bandpass filter using a Butterworth filter and then applies it to your AFMT data to remove frequencies outside the specified range.

Adjust the lowcut and highcut parameters according to the desired frequency range for your data. Removing bad frequencies from data typically involves filtering the data to eliminate unwanted noise or artifacts.

Parameters:

• data (arraylike 1D) – Noise data to filter.

• freqs (Arraylike 1d) – Array of frequencies onto apply the bandpass filter.

• fs (int,) – Sample of frequencies. If None, use the number of original frequency

• frange (list , Optional) – frequency range ( min/200., max/5) for the bandpass filter (in Hz). By default, use the minimum and maximum of original frquency array. Note that digital filter critical frequencies must be 0 < Wn < 1 i.e.

• order (int, default=5) – Order for butter bandpass.

• plot (bool, default=False) – Visualize the filtered data.

Returns:

y

Return type:

filtered data

Example

>>> import numpy as np
>>> from watex.utils.exmath import butterworth_filter
>>> time = np.linspace(0, 1, 1000)  # Replace with your actual time values
>>> freqs = np.linspace ( 1, 1000, 500)
>>> data = np.sin(2 * np.pi * 10 *freqs) + 0.5 * np.sin(2 * np.pi * 50 *freqs)
>>> _=butterworth_filter (data , freqs , fs = 1000, frange=( 5, 20), plot=True )

watex.utils.exmath.compute_anr(sfi, rhoa_array, pos_bound_indexes)

Compute the select anomaly ratio (ANR) along with the whole profile from SFI.

The standardized resistivity values`rhoa` of is averaged from \(X_{begin}\) to \(X_{end}\). The ANR is a positive value and the equation is given as:

\[ANR= sfi * \frac{1}{N} * | \sum_{i=1}^{N} \frac{ \rho_{a_i} - \bar \rho_a}{\sigma_{\rho_a}} |\]

where \(\sigma_{rho_a}\) and \(\bar \rho_a\) are the standard deviation and the mean of the resistivity values composing the selected anomaly.

Parameters:

• sfi – Is standard fracturation index. please refer to compute_sfi().

• rhoa_array (array_like) – Resistivity values of Electrical Resistivity Profiling line

• pos_bound_indexes – Select anomaly station location boundaries indexes. Refer to compute_power() of pos_bounds.

Returns:

Anomaly ratio

Return type:

float

Example:

>>> from watex.utils.exmath import compute_anr
>>> import pandas as pd
>>> anr = compute_anr(sfi=sfi,
...                  rhoa_array=data = pd.read_excel(
...                  'data/l10_gbalo.xlsx').to_numpy()[:, -1],
...              pk_bound_indexes  = [9, 13])
>>> anr

watex.utils.exmath.compute_errors(arr, /, error='std', axis=0, return_confidence=False)

Compute Errors ( Standard Deviation ) and standard errors.

Standard error and standard deviation are both measures of variability: - The standard deviation describes variability within a single sample. Its

formula is given as:

\[SD = \sqrt{ \sum |x -\mu|^2}{N}\]

where \(\sum\) means the “sum of”, \(x\) is the value in the data set,:math:mu is the mean of the data set and \(N\) is the number of the data points in the population. \(SD\) is the quantity expressing by how much the members of a group differ from the mean value for the group.

• The standard error estimates the variability across multiple samples of a population. Different formulas are used depending on whether the population standard deviation is known.

  • when the population standard deviation is known:

    \[SE =\]

rac{SD}{sqrt{N}}

• When the population parameter is unknwon

  \[SE =\]

rac{s}{sqrt{N}}

where \(SE\) is the standard error, : math:s is the sample standard deviation. When the population standard is knwon the \(SE\) is more accurate.

Note that the \(SD\) is a descriptive statistic that can be calculated from sample data. In contrast, the standard error is an inferential statistic that can only be estimated (unless the real population parameter is known).

Parameters:

arrarray_like , 1D or 2D

Array for computing the standard deviation

error: str, default=’std’

Name of error to compute. By default compute the standard deviation. Can also compute the the standard error estimation if the argument is passed to ste.

return_confidence: bool, default=False,

If True, returns the confidence interval with 95% of sample means

Returns:

err: arraylike 1D or 2D

Error array.

Examples

>>> from watex.datasets import load_huayuan
>>> from watex.utils.exmath import compute_errors
>>> emobj=load_huayuan ().emo
>>> compute_errors (emobj.freqs_ )
.. Out[104]: 14397.794665683341
>>> freq2d = emobj.make2d ('freq')
>>> compute_errors (freq2d ) [:7]
array([14397.79466568, 14397.79466568, 14397.79466568, 14397.79466568,
       14397.79466568, 14397.79466568, 14397.79466568])
>>> compute_errors (freq2d , error ='se') [:7]
array([1959.29168624, 1959.29168624, 1959.29168624, 1959.29168624,
       1959.29168624, 1959.29168624, 1959.29168624])

watex.utils.exmath.compute_lower_anomaly(erp_array, station_position=None, step=None, **kws)

Function to get the minimum value on the ERP array.

If pk is provided wil give the index of pk.

Parameters:

• erp_array (array_like) – array of apparent resistivity profile

• position (station) – array of station position (survey), if not given and step is known , set the step value and station_position will compute automatically

• step – The distance between measurement im meter. If given will recompute the station_position

Returns:

• bestSelectedDict: dict containing best anomalies

  with the anomaly resistivities range.

• anpks: Main positions of best select anomaly

• collectanlyBounds: list of arrays of select anomaly values

• min_pks: list of tuples (pk, minVal of best anomalies points.)

Return type:

tuple

Example:

>>> from watex.utils.exmath import compute_lower_anolamy
>>> import pandas as pd
>>> erp_data= 'data/l10_gbalo.xlsx'
>>> dataRes=pd.read_excel(erp_data).to_numpy()[:,-1]
>>> anomaly, *_ =  compute_lower_anomaly(erp_array=data, step =10)
>>> anomaly

watex.utils.exmath.compute_magnitude(rhoa_max=None, rhoa_min=None, rhoaMinMax=None)

Compute the magnitude Ma of selected anomaly expressed in Ω.m. ano :param rhoa_min: resistivity value of selected anomaly :type rhoa_min: float

Parameters:

rhoa_max (float) – Max boundary of the resistivity value of select anomaly.

Returns:

The absolute value between the rhoa_min and rhoa_max.

Return type:

float

Example:

>>> from watex.utils.exmath import compute_power
>>> power= compute_power(80, 130)

watex.utils.exmath.compute_power(posMinMax=None, pk_min=None, pk_max=None)

Compute the power Pa of anomaly.

Parameters:

• pk_min (float) – Min boundary value of anomaly. pk_min is min value (lower) of measurement point. It’s the position of the site in meter.

• pk_max (float) – Max boundary of the select anomaly. pk_max is the maximum value the measurement point in meter. It’s the upper boundary position of the anomaly in the site in m.

Returns:

The absolute value between the pk_min and pk_max.

Return type:

float

Example:

>>> from watex.utils.exmath import compute_power
>>> power= compute_power(80, 130)

watex.utils.exmath.compute_sfi(pk_min, pk_max, rhoa_min, rhoa_max, rhoa, pk)

SFI is introduced to evaluate the ratio of presumed existing fracture from anomaly extent. We use a similar approach as IF computation proposed by Dieng et al. (2004) to evaluate each selected anomaly extent and the normal distribution of resistivity values along the survey line. The SFI threshold is set at \(sqrt(2)\) for symmetrical anomaly characterized by a perfect distribution of resistivity in a homogenous medium.

Parameters:

• pk_min – see compute_power()

• pk_max – see compute_power()

• rhoa_max – see compute_magnitude()

• rhoa_min – see compute_magnitude()

• pk – Station position of the selected anomaly in float value.

• rhoa – Selected anomaly apparent resistivity value in ohm.m

Returns:

standard fracture index (SFI)

Return type:

float

Example:

>>> from watex.utils.exmath import compute_sfi
>>> sfi = compute_sfi(pk_min = 90,
...                      pk_max=130,
...                      rhoa_min=175,
...                      rhoa_max=170,
...                      rhoa=132,
...                      pk=110)
>>> sfi

watex.utils.exmath.convert_distance_to_m(value, converter=1000.0, unit='km')

Convert distance from km to m or vice versa even a string value is given.

Parameters:

• value – value to convert.

• unit – unit to convert to.

Paramm converter:

Equivalent if given in km rather than m.

watex.utils.exmath.d_hanning_window(x, xk, W)

Discrete hanning function.

For futher details, please refer to https://doi.org/10.1190/1.2400625

Parameters:

• x – variable point along the window width

• xk – Center of the window W. It presumes to host the most weigth.

• W – int, window-size; preferably set to odd number. It must be less than the dipole length.

Returns:

Anonymous function (x,xk, W) value

watex.utils.exmath.define_anomaly(erp_data, station_position=None, pks=None, dipole_length=10.0, **kwargs)

Function will select the different anomalies. If pk is not given, the best three anomalies on the survey lines will be computed automatically

Parameters:

• erp_data (array_like) – Electrical resistivity profiling

• pks (list or dict) – station positions anomaly boundaries (pk_begin, pk_end) If selected anomalies is more than one, set pks into dict where number of anomaly =keys and pks = values

• dipole_length (float) – Distance between two measurements in meters Change the `dipole lengh

• station_position – station position array

Returns:

list of anomalies bounds

watex.utils.exmath.define_conductive_zone(erp, stn=None, sres=None, *, distance=None, dipole_length=None, extent=7)

Detect the conductive zone from `s`ves point.

Parameters:

• erp – Resistivity values of electrical resistivity profiling(ERP).

• stn – Station number expected for VES and/or drilling location.

• sres – Resistivity value at station number stn. If sres is given, the auto-search will be triggered to find the station number that fits the resistivity value.

• distance – Distance from the first station to stn. If given, be sure to provide the dipole_length

• dipole_length – Length of the dipole. Comonly the distante between two close stations. Since we use config AB/2

• extent – Is the width to depict the anomaly. If provide, need to be consistent along all ERP line. Should keep unchanged for other parameters definitions. Default is 7.

Returns:

• CZ:Conductive zone including the station position

• sres: Resistivity value of the station number

• ix_stn: Station position in the CZ

Note

If many stations got the same sres value, the first station is flagged. This may not correspond to the station number that is searching. Use sres only if you are sure that the resistivity value is unique on the whole ERP. Otherwise it’s not recommended.

Example:

>>> import numpy as np
>>> from watex.utils.exmath import define_conductive_zone
>>> sample = np.random.randn(9)
>>> cz, stn_res = define_conductive_zone(sample, 4, extent = 7)
... (array([ 0.32208638,  1.48349508,  0.6871188 , -0.96007639,
            -1.08735204,0.79811492, -0.31216716]),
     -0.9600763919368086,
     3)

watex.utils.exmath.detect_station_position(s, p)

Detect station position and return the index in positions

Parameters:

• s – str, int - Station location in the position array. It should be the positionning of the drilling location. If the value given is type string. It should be match the exact position to locate the drilling. Otherwise, if the value given is in float or integer type, it should be match the index of the position array.

• p – Array-like - Should be the conductive zone as array of station location values.

Returns:

• s_index- the position index location in the conductive zone.

• s- the station position in distance.

Example:

>>> import numpy as np
>>> from watex.utils.exmath import detect_station_position
>>> pos = np.arange(0 , 50 , 10 )
>>> detect_station_position (s ='S30', p = pos)
... (3, 30.0)
>>> detect_station_position (s ='40', p = pos)
... (4, 40.0)
>>> detect_station_position (s =2, p = pos)
... (2, 20)
>>> detect_station_position (s ='sta200', p = pos)
... WATexError_station: Station sta200             is out of the range; max position = 40

watex.utils.exmath.dummy_basement_curve(func, ks, slope=45)

Compute the pseudodepth from the search zone.

Parameters:

• f – callable - Polyfit1D function

• mz – array-zone - Expected Zone for groundwater search

• ks – float - The depth from which the expected fracture zone must starting looking for groundwater.

• slope – float - Degree angle for slope in linear function of the dummy curve

Returns:

• lambda function of basement curve func45

• beta is intercept value compute for keysearch ks

watex.utils.exmath.exportEDIs(ediObjs, new_Z, savepath=None, **kws)

Export EDI files from multiples EDI or z objects

Export new EDI files from the former object with a given new impedance tensors. The export is assumed a new output EDI resulting from multiples corrections applications.

Parameters:

• ediObjs (list of string watex.edi.Edi) – Full path to Edi file/object or object from class:EM objects.

• new_Z (list of ndarray (nfreq, 2, 2)) – A collection of Ndarray of impedance tensors Z. The tensor Z is 3D array composed of number of frequency nfreq`and four components (``xx`, xy, yx, and yy) in 2X2 matrices. The tensor Z is a complex number.

• savepath (str, Optional) – Path to save a new EDI file. If None, outputs to the current directory.

Return type:

ediObj from watex.edi.Edi

See also

exportedi

Export single EDI from

watex.utils.exmath.find_bound_for_integration(ix_arr, b0=[])

Recursive function to find the roots between f curve and basement curves so to detect the integration bounds.

The function use entirely numpy for seaching integration bound. Since it is much faster than find_limit_for_integration() although both did the same tasks.

Parameters:

• ix_arr – array-like - Indexes array from masked array where the value are true i.e. where \(b-f > 0 \Rightarrow b > f\) .

• b0 – list - Empy list to hold the limit during entire loop

Returns:

list - integration bounds

Note

\(b > f \Longrightarrow\) Curve b (basement) is above the fitting curve \(f\) . \(b < f\) otherwise. The pseudoarea is the area where \(b > f\) .

watex.utils.exmath.find_closest(arr, /, values)

Get the closest value in array from given values.

Parameters:

• arr (Arraylike) – Array to find the values

• values (float, arraylike) –

Return type:

closest values in float or array containing in the given array.

Examples

>>> import numpy as np
>>> from watex.utils.exmath import find_closest
>>> find_closest (  [ 2 , 3, 4, 5] , ( 2.6 , 5.6 )  )
array([3., 5.])
>>> find_closest (  np.array ([[2 , 3], [ 4, 5]]), ( 2.6 , 5.6 ) )
array([3., 5.])
array([3., 5.])

watex.utils.exmath.find_limit_for_integration(ix_arr, b0=[])

Use the roots between f curve and basement curves to detect the limit of integration.

Parameters:

• ix_arr – array-like - Indexes array from masked array where the value are true i.e. where :math:` b-f > 0 Rightarrow b> f` .

• b0 – list - Empy list to hold the limit during entire loop

Note

\(b > f \Longrightarrow\) Curve b (basement) is above the fitting curve \(f\) . \(b < f\) otherwise. The pseudoarea is the area where :math:` b > f` .

Returns:

list - integration bounds

watex.utils.exmath.fitfunc(x, y, deg=None, sample=1000)

Create polyfit function from a specifc sample data points.

Parameters:

• x – array-like of x-axis.

• y – array-like of y-axis.

• deg – polynomial degree. If None should compute using the length of extrema (local + global).

• sample – int - Number of data points should use for fitting function. Default is 1000.

Returns:

• Polynomial function f

• new axis x_new generated from the samples.

• projected sample values got from f.

watex.utils.exmath.fittensor(refreq, compfreq, z, fill_value=nan)

Fit each tensor component to the complete frequency range.

The complete frequency is the frequency with clean data. It contain all the frequency range on the site. During the survey, the missing frequencies lead to missing tensor data. So the function will indicate where the tensor data is missing and fit to the prior frequencies.

Parameters:

• refreq (ArrayLike) – Reference frequency - Should be the complete frequency collected in the field.

• comfreq (array-like,) – The specific frequency collect in the site. Sometimes due to the interferences, the frequency at individual site could be different from the complete. However, the frequency values at the individual site must be included in the complete frequency refreq.

• z (array-like,) – should be the tensor value (real or imaginary part ) at the component xx, xy, yx, yy.

• fill_value (float . default='NaN') – Value to replace the missing data in tensors.

Returns:

Z – new Z filled by invalid value NaN where the frequency is missing in the data.

Return type:

Arraylike

Examples

>>> import numpy as np
>>> from watex.utils.exmath import fittensor
>>> refreq = np.linspace(7e7, 1e0, 20) # 20 frequencies as reference
>>> freq_ = np.hstack ((refreq.copy()[:7], refreq.copy()[12:] ))
>>> z = np.random.randn(len(freq_)) *10 # assume length of  freq as
...                 # the same like the tensor Z value
>>> zn  = fittensor (refreq, freq_, z)
>>> z # some frequency values are missing but not visible.
...array([-23.23448367,   2.93185982,  10.81194723, -12.46326732,
         1.57312908,   7.23926576, -14.65645799,   9.85956253,
         3.96269863, -10.38325124,  -4.29739755,  -8.2591703 ,
        21.7930423 ,   0.21709129,   4.07815217])
>>> # zn show where the frequencies are missing
>>> # the NaN value means in a missing value in  tensor Z at specific frequency
>>> zn
... array([-23.23448367,   2.93185982,  10.81194723, -12.46326732,
         1.57312908,   7.23926576, -14.65645799,          nan,
                nan,          nan,          nan,          nan,
         9.85956253,   3.96269863, -10.38325124,  -4.29739755,
        -8.2591703 ,  21.7930423 ,   0.21709129,   4.07815217])
>>> # let visualize where the missing frequency value in tensor Z
>>> refreq
... array([7.00000000e+07, 6.63157895e+07, 6.26315791e+07, 5.89473686e+07,
       5.52631581e+07, 5.15789476e+07, 4.78947372e+07, 4.42105267e+07*,
       4.05263162e+07*, 3.68421057e+07*, 3.31578953e+07*, 2.94736848e+07*,
       2.57894743e+07, 2.21052638e+07, 1.84210534e+07, 1.47368429e+07,
       1.10526324e+07, 7.36842195e+06, 3.68421147e+06, 1.00000000e+00])
>>> refreq[np.isnan(zn)] #we can see the missing value between [7:12](*) in refreq
... array([44210526.68421052, 40526316.21052632, 36842105.73684211,
       33157895.2631579 , 29473684.78947368])

watex.utils.exmath.forward_model(AB, resistivity, thickness)

Placeholder forward model function. This function should compute the theoretical apparent resistivity values based on the resistivity and thickness of subsurface layers and compare them with observed values.

Parameters:

• AB (np.ndarray) – Array containing the spacing of current electrodes.

• resistivity (float or np.ndarray) – Resistivity of subsurface layers.

• thickness (float or np.ndarray) – Thickness of subsurface layers.

Returns:

Predicted apparent resistivity values.

Return type:

np.ndarray

watex.utils.exmath.get2dtensor(z_or_edis_obj_list, /, tensor='z', component='xy', kind='modulus', return_freqs=False, **kws)

Make tensor into two dimensional array from a collection of Impedance tensors Z.

Out 2D resistivity, phase-error and tensor matrix from a collection of EDI-objects.

Matrix depends of the number of frequency times number of sites. The function asserts whether all data from all frequencies are available. The missing values should be filled by NaN. Note that each element of z is (nfreq, 2, 2) dimension for:

xx ( 0, 0) ------- xy ( 0, 1)
yx ( 1, 0) ------- yy ( 1, 1)

Parameters:

• z_or_edis_obj_list (list of watex.edi.Edi or watex.externals.z.Z) – A collection of EDI- or Impedances tensors objects.

• tensor (str, default='z') – Tensor name. Can be [ resistivity|phase|z|frequency]

• component (str, default='xy' (TE mode)) – EM mode. Can be [‘xx’, ‘xy’, ‘yx’, ‘yy’]

• out (str) – kind of data to output. Be sure to provide the component to retrieve the attribute from the collection object. Except the error and frequency attribute, the missing component to the attribute will raise an error. for instance resxy for xy component. Default is resxy.

• kind (str , default='modulus') – focuses on the tensor output. Note that the tensor is a complex number of ndarray (nfreq, 2,2 ). If set to``modulus`, the modulus of the complex tensor should be outputted. If real or``imag``, it returns only the specific one. Default is complex.

• return_freqs (Arraylike ,) – If True , returns also the full frequency ranges.

• kws (dict) – Additional keywords arguments from :meth:`~EM.getfullfrequency `.

Returns:

mat2d – the matrix of number of frequency and number of Edi-collectes which correspond to the number of the stations/sites.

Return type:

arraylike2d

Examples

>>> from watex.datasets import load_huayuan
>>> from watex.methods import get2dtensor
>>> box= load_huayuan ( key ='raw', clear_cache = True, samples =7)
>>> data = box.data
>>> phase_yx = get2dtensor ( data, tensor ='phase', component ='yx')
>>> phase_yx.shape
(56, 7)
>>> phase_yx [0, :]
array([        nan,         nan,         nan,         nan, 18.73244951,
       35.00516522, 59.91093054])

watex.utils.exmath.get_anomaly_ratio(erp, czposix=None, cz=None, cz_sfi=None, raise_exception=True, p=None, e_spacing=None, **sfi_kws)

Computes the selected anomaly ratio (ANR) from the whole ERP line.

The standardized resistivity values`rhoa` of is averaged from \(S_{begin}\) to \(S_{end}\). The ANR is a positive value and the equation is given as:

\[ANR= sfi * \frac{1}{N} * | \sum_{i=1}^{N} \frac{\rho_{a_i} - \bar \rho_a}{\sigma_{\rho_a}} |\]

where \(\sigma_{rho_a}\) and \(\bar \rho_a\) are the standard deviation and the mean of the resistivity values composing the selected anomaly.

Parameters:

• erp (array_like 1d) – The ERP survey line. The line is an array of resistivity values. Note that if a dataframe is passed, be sure that the frame matches the DC resistivity data (ERP), otherwise an error occurs. At least, the frame columns includes the resistivity and stations.

• cz_sfi (float,) – The pseudo-fracturing index value. It can be computed from sfi(). It is given , p and e_spacing are not needed.

• czposix (list of int,) – The selected anomaly station location boundaries indexes. The indexes might correspond to the first and last stations indexes that represents the selected conductive zone. For instance, czposix=[2, 13] means the third ( second+ 1) stations to the 14 (13+1) th stations. Note that the index counts is Python indexes so it starts by 0.

• p (arraylike,) – is the station positions of the whole ERP line. It must be consistent with the ERP line.

• e_spacing (float, int,) – The electrode spacing. It is needed in complience with p especially when the czposix is not supplied.

• raise_exception (bool, default= True,) – Raise exception when the czposix is not given. However another alternative way when the p is not given too, is to use the cz resistivity values from detecting the czposix, however this has a risk of biais the position then raises an exception for user to be aware. Note that user can force this approach to take effect by setting raise_exception to False.

cz: array_like 1d

the selected conductive zone. If None, only the erp should be displayed. Note that cz is an subset of erp array.

Examples

>>> from watex.datasets import make_erp
>>> from watex.utils.exmath import get_anomaly_ratio
>>>  # for data reproducibility seed to 123
>>> d = make_erp (n_stations =70, min_rho =10 , max_rho =1e4 , seed =123 )
>>> selected_cz,* _= defineConductiveZone (d.resistivity , auto=True)
>>> ANR = get_anomaly_ratio (d.resistivity, cz= selected_cz , e_spacing =10)
>>> print(ANR)
... 0.06 #6%

watex.utils.exmath.get_azimuth(xlon, ylat, *, data=None, utm_zone=None, projection='ll', isdeg=True, mode='soft', extrapolate=Ellipsis, view=Ellipsis)

Compute azimuth from coordinate locations ( latitude, longitude).

If easting and northing are given rather than longitude and latitude, the projection should explicitely set to UTM to perform the ideal conversion. However if mode is set to soft (default), the type of projection is automatically detected . Note that when UTM coordinates are provided, xlon and ylat fit easting and northing respectively.

Parameters:

• xlon (Arraylike 1d or str, str) – ArrayLike of easting/longitude and arraylike of nothing/latitude. They should be one dimensional. In principle if data is supplied, they must be series. If xlon and ylat are given as string values, the data must be supplied. xlon and ylat names must be included in the dataframe otherwise an error raises.

• ylat (Arraylike 1d or str, str) – ArrayLike of easting/longitude and arraylike of nothing/latitude. They should be one dimensional. In principle if data is supplied, they must be series. If xlon and ylat are given as string values, the data must be supplied. xlon and ylat names must be included in the dataframe otherwise an error raises.

• data (pd.DataFrame,) – Data containing x and y names. Need to be supplied when x and y are given as string names.

• utm_zone (Optional, string) – zone number and ‘S’ or ‘N’ e.g. ‘55S’. Default to the centre point of coordinates points in the survey area. It should be a string (##N or ##S) in the form of number and North or South hemisphere, 10S or 03N

• projection (str, ['utm'|'ll']) – The coordinate system in which the data points for the profile is collected. when mode=’soft’, the auto-detection will be triggered and find the suitable coordinate system. However, it is recommended to explicitly provide projection when data is in UTM coordinates. Note that if x and y are composed of value greater than 180 degrees for longitude and 90 degrees for latitude, and method is still in the soft` mode, it should be considered as  longitude-latitude ``UTM coordinates system.

• isdeg (bool, default=True) – By default xlon and xlat are in degree coordinates. If both arguments are given in radians, set to False instead.

• mode (str , ['soft'|'strict']) – strict mode does not convert any coordinates system to other at least it is explicitly set to projection whereas the soft does.

• extrapolate (bool, default=False) – In principle, the azimuth is compute between two points. Thus, the number of values computed for \(N\) stations should be \(N-1\). To fit values to match the number of size of the array, extrapolate should be True. In that case, the first station holds a <<fake>> azimuth as the closer value computed from interpolation of all azimuths.

• view (bool, default=False,) – Quick view of the azimuth. It is usefull especially when extrapolate is set to True.

Returns:

azim – Azimuth computed from locations.

Return type:

ArrayLike

Examples

>>> import watex as wx
>>> from watex.utils.exmath import get_azimuth
>>> # generate a data from ERP
>>> data = wx.make_erp (n_stations =7 ).frame
>>> get_azimuth ( data.longitude, data.latitude)
array([54.575, 54.575, 54.575, 54.575, 54.575, 54.575])
>>> get_azimuth ( data.longitude, data.latitude, view =True, extrapolate=True)
array([54.57500007, 54.575     , 54.575     , 54.575     , 54.575     ,
       54.575     , 54.575     ])

watex.utils.exmath.get_bearing(latlon1, latlon2, to_deg=True)

Calculate the bearing between two points.

A bearing can be defined as a direction of one point relative to another point, usually given as an angle measured clockwise from north. The formula of the bearing \(eta\) between two points 1(lat1 , lon1) and 2(lat2, lon2) is expressed as below:

\[eta = atan2(sin(y_2-y_1)*cos(x_2), cos(x_1)*sin(x_2) – sin(x_1)*cos(x_2)*cos(y_2-y_1))\]

where:

• :math:`x_1`(lat1): the latitude of the first coordinate

• :math:`y_1`(lon1): the longitude of the first coordinate

• :math:`x_2`(lat2) : the latitude of the second coordinate

• :math:`y_2`(lon2): the longitude of the second coordinate

Parameters:

• latlon (Tuple ( latitude, longitude)) – A latitude and longitude coordinates of the first point in degree.

• latlon2 (Tuple ( latitude, longitude)) – A latitude and longitude of coordinates of the second point in degree.

• to_deg (bool, default=True) – Convert the bearing from radians to degree.

Returns:

• b (Value of bearing in degree ( default).)

• See More

• ———-

• See more details by clicking in the link below – https://mapscaping.com/how-to-calculate-bearing-between-two-coordinates/

Examples

>>> from watex.utils import get_bearing
>>> latlon1 = (28.41196763902007, 109.3328724432221) # (lat, lon) point 1
>>> latlon2= (28.38756530909265, 109.36931920880758) # (lat, lon) point 2
>>> get_bearing (latlon1, latlon2 )
127.26739270447973 # in degree

watex.utils.exmath.get_distance(x, y, *, return_mean_dist=False, is_latlon=False, **kws)

Compute distance between points

Parameters:

• x (ArrayLike 1d,) – One dimensional arrays. x can be consider as the abscissa of the landmark and y as ordinates array.

• y (ArrayLike 1d,) – One dimensional arrays. x can be consider as the abscissa of the landmark and y as ordinates array.

• return_mean_dist (bool, default =False,) – Returns the average value of the distance between different points.

• is_latlon (bool, default=False,) – Convert x and y latitude and longitude coordinates values into UTM before computing the distance. x, y should be considered as easting and northing respectively.

• kws (dict,) – Keyword arguments passed to watex.site.Location.to_utm_in()

Returns:

d – Is the distance between points.

Return type:

Arraylike of shape (N-1)

Examples

>>> import numpy as np
>>> from watex.utils.exmath import get_distance
>>> x = np.random.rand (7) *10
>>> y = np.abs ( np.random.randn (7) * 12 )
>>> get_distance (x, y)
array([ 8.7665511 , 12.47545656,  8.53730212, 13.54998351, 14.0419387 ,
       20.12086781])
>>> get_distance (x, y, return_mean_dist= True)
12.91534996818084

watex.utils.exmath.get_full_frequency(z_or_edis_obj_list, /, to_log10=False)

Get the frequency with clean data.

The full or plain frequency is array frequency with no missing frequency during the data collection. Note that when using Natural Source Audio-Magnetotellurics, some data are missing due to the weak of missing frequency at certain band especially in the attenuation band.

Parameters:

• z_or_edis_obj_list (list of watex.edi.Edi or watex.externals.z.Z) – A collection of EDI- or Impedances tensors objects.

• to_log10 (bool, default=False) – Export frequency to base 10 logarithm

Returns:

f – frequency with clean data. Out of attenuation band if survey is completed with Natural Source Audio-Magnetotellurics.

Return type:

Arraylike of shape(N, )

Examples

>>> from watex.datasets import load_huayuan
>>> from watex.methods.em import get_full_frequency
>>> box= load_huayuan ( key ='raw', clear_cache = True, samples =7)
>>> edi_data = box.data
>>> f = get_full_frequency (edi_data )
>>> f
array([8.19200e+04, 7.00000e+04, 5.88000e+04, 4.95000e+04, 4.16000e+04,
       3.50000e+04, 2.94000e+04, 2.47000e+04, 2.08000e+04, 1.75000e+04,
       ...
       3.25000e+01, 2.75000e+01, 2.25000e+01, 1.87500e+01, 1.62500e+01,
       1.37500e+01, 1.12500e+01, 9.37500e+00, 8.12500e+00, 6.87500e+00,
       5.62500e+00])
>>> len(f)
56
>>> # Get only the z component objects
>>> zobjs = [ box.emo.ediObjs_[i].Z for i in  range (len(box.emo.ediObjs_))]
>>> len(zobjs)
56

watex.utils.exmath.get_minVal(array)

Function to find the three minimum values on array and their corresponding indexes.

Parameters:

array (array_like) – array of values

Returns:

Three minimum values of rho, index in rho_array

Return type:

tuple

watex.utils.exmath.get_profile_angle(easting=None, northing=None, msg='ignore')

compute geoprofile angle. :param * easting: easting coordiantes values :type * easting: array_like :param * northing: northing coordinates values :type * northing: array_like :param * msg: :type * msg: output a little message if msg is set to “raises”

Returns:

• float – profile_angle

• float – geo_electric_strike

watex.utils.exmath.get_shape(rhoa_range)

Find anomaly shape from apparent resistivity values framed to the best points using the mean line.

Parameters:

rhoa_range (array_like or list) – The apparent resistivity from selected anomaly bounds anom_boundaries

Returns:

• V

• W

• K

• C

• M

• U

Example:

>>> from watex.utils.exmath import get_shape
>>> x = [60, 70, 65, 40, 30, 31, 34, 40, 38, 50, 61, 90]
>>> shape = get_shape (rhoa_range= np.array(x))
... U

watex.utils.exmath.get_station_number(dipole, distance, from0=False, **kws)

Get the station number from dipole length and the distance to the station.

Parameters:

• distance – Is the distance from the first station to s in meter (m). If value is given, please specify the dipole length in the same unit as distance.

• dipole – Is the distance of the dipole measurement. By default the dipole length is in meter.

• kws – convert_distance_to_m() additional arguments

watex.utils.exmath.get_strike(profile_angle=None, easting=None, northing=None, gstrike=None, msg='ignore')

Compute geoelectric strike from profile angle, easting and northing.

Parameters:

• profile_angle (*) – If not provided , will comput with easting and northing coordinates

• easting (*) – Easting coordiantes values

• northing (*) – Northing coordinates values

• gstrike (*) – strike value , if provided, will recomputed geo_electric strike . * msg: output a little message if msg is set to “raises”

Returns:

• float – profile_angle in degree E of N

• float – geo_electric_strike in degrees E of N

watex.utils.exmath.get_type(erp_array, posMinMax, pk, pos_array, dl)

Find anomaly type from app. resistivity values and positions locations

Parameters:

• erp_array (array_like) – App.resistivty values of all erp lines

• posMinMax (list or tuple or nd.array(1,2)) – Selected anomaly positions from startpoint and endpoint

• pk (float or int) – Position of selected anomaly in meters

• pos_array (array_like) – Stations locations or measurements positions

• dl – Distance between two receiver electrodes measurement. The same as dipole length in meters.

Returns:

• EC for Extensive conductive.

• NC for narrow conductive.

• CP for conductive plane

• CB2P for contact between two planes.

Example:

>>> from watex.utils.exmath import get_type
>>> x = [60, 61, 62, 63, 68, 65, 80,  90, 100, 80, 100, 80]
>>> pos= np.arange(0, len(x)*10, 10)
>>> ano_type= get_type(erp_array= np.array(x),
...            posMinMax=(10,90), pk=50, pos_array=pos, dl=10)
>>> ano_type
...CB2P

watex.utils.exmath.get_type2(erp_array, posMinMax, pk, pos_array, dl=None)

Find anomaly type from app. resistivity values and positions locations

Parameters:

• erp_array (array_like) – App.resistivty values of all erp lines

• posMinMax (list or tuple or nd.array(1,2)) – Selected anomaly positions from startpoint and endpoint

• pk (float or int) – Position of selected anomaly in meters

• pos_array (array_like) – Stations locations or measurements positions

• dl – Distance between two receiver electrodes measurement. The same as dipole length in meters.

Returns:

• EC for Extensive conductive.

• NC for narrow conductive.

• CP for conductive plane

• CB2P for contact between two planes.

Example:

>>> from watex.core.erp import get_type
>>> x = [60, 61, 62, 63, 68, 65, 80,  90, 100, 80, 100, 80]
>>> pos= np.arange(0, len(x)*10, 10)
>>> ano_type= get_type(erp_array= np.array(x),
...            posMinMax=(10,90), pk=50, pos_array=pos, dl=10)
>>> ano_type
...CB2P

watex.utils.exmath.get_z_from(edi_obj_list, /)

Extract z object from Edi object.

Parameters:

z_or_edis_obj_list (list of watex.edi.Edi or watex.externals.z.Z) – A collection of EDI- or Impedances tensors objects.

Returns:

Z – List of impedance tensor Objects.

Return type:

list of watex.externals.z.Z

watex.utils.exmath.getshape(rhoa_range)

Find anomaly shape from apparent resistivity values framed to the best points.

Parameters:

rhoa_range (array_like or list) – The apparent resistivity from selected anomaly bounds anom_boundaries

Returns:

• V

• W

• K

• C

• M

• U

Example:

>>> from watex.core.erp import get_shape
>>> x = [60, 70, 65, 40, 30, 31, 34, 40, 38, 50, 61, 90]
>>> shape = get_shape (rhoa_range= np.array(x))
...U

watex.utils.exmath.gettype(erp_array, posMinMax, pk, pos_array, dl)

Find anomaly type from app. resistivity values and positions locations

Parameters:

• erp_array (array_like) – App.resistivty values of all erp lines

• posMinMax (list or tuple or nd.array(1,2)) – Selected anomaly positions from startpoint and endpoint

• pk (float or int) – Position of selected anomaly in meters

• pos_array (array_like) – Stations locations or measurements positions

• dl – Distance between two receiver electrodes measurement. The same as dipole length in meters.

Returns:

• EC for Extensive conductive.

• NC for narrow conductive.

• CP for conductive plane

• CB2P for contact between two planes.

Example:

>>> from watex.methods.erp import get_type
>>> x = [60, 61, 62, 63, 68, 65, 80,  90, 100, 80, 100, 80]
>>> pos= np.arange(0, len(x)*10, 10)
>>> ano_type= get_type(erp_array= np.array(x),
...            posMinMax=(10,90), pk=50, pos_array=pos, dl=10)
>>> ano_type
...CB2P

watex.utils.exmath.gradient_descent(z, s, alpha=0.01, n_epochs=100, kind='linear', degree=1, raise_warn=False)

Gradient descent algorithm to fit the best model parameter.

Model can be changed to polynomial if degree is greater than 1.

Parameters:

• z (arraylike,) – vertical nodes containing the values of depth V

• s (Arraylike,) – vertical vector containin the resistivity values

• alpha (float,) – step descent parameter or learning rate. Default is ``0.01`

• n_epochs (int,) – number of iterations. Default is 100. Can be changed to other values

• kind (str, {"linear", "poly"}, default= 'linear') – Type of model to fit. Linear model is selected as the default.

• degree (int, default=1) – As the linear model is selected as the default since the degree is set to 1

Returns:

• - `F` (New model values with the best W parameters found.)

• - `W` (vector containing the parameters fits)

• - `cost_history` (Containing the error at each Itiretaions.)

Examples

>>> import numpy as np
>>> from watex.utils.exmath import gradient_descent
>>> z= np.array([0, 6, 13, 20, 29 ,39, 49, 59, 69, 89, 109, 129,
                 149, 179])
>>> res= np.array( [1.59268,1.59268,2.64917,3.30592,3.76168,
                    4.09031,4.33606, 4.53951,4.71819,4.90838,
      5.01096,5.0536,5.0655,5.06767])
>>> fz, weights, cost_history = gradient_descent(
    z=z, s=res,n_epochs=10,alpha=1e-8,degree=2)
>>> import matplotlib.pyplot as plt
>>> plt.scatter (z, res)
>>> plt.plot(z, fz)

watex.utils.exmath.hmcmc_inversion(AB, rhoa, rho0, h0, **kwargs)

Perform VES data inversion using Hybrid Monte Carlo (HMCMC).

Parameters:

• AB (np.ndarray) – Array containing the spacing of current electrodes.

• rhoa (np.ndarray) – Array containing apparent resistivity values.

• rho0 (float) – Initial guess for resistivity.

• h0 (float) – Initial guess for thickness of the first layer.

Returns:

Inversion results including estimated resistivity and thickness.

Return type:

np.ndarray

watex.utils.exmath.interpolate1d(arr, kind='slinear', method=None, order=None, fill_value='extrapolate', limit=None, **kws)

Interpolate array containing invalid values NaN

Usefull function to interpolate the missing frequency values in the tensor components.

Parameters:

• arr (array_like) – Array to interpolate containg invalid values. The invalid value here is NaN.

• kind (str or int, optional) – Specifies the kind of interpolation as a string or as an integer specifying the order of the spline interpolator to use. The string has to be one of linear, nearest, nearest-up, zero, slinear,``quadratic``, cubic, previous, or next. zero, slinear, quadratic``and ``cubic refer to a spline interpolation of zeroth, first, second or third order; previous and next simply return the previous or next value of the point; nearest-up and nearest differ when interpolating half-integers (e.g. 0.5, 1.5) in that nearest-up rounds up and nearest rounds down. If method param is set to pd which refers to pd.interpolate method , kind can be set to polynomial or pad interpolation. Note that the polynomial requires you to specify an order while pad requires to specify the limit. Default is slinear.

• method (str, optional, default='mean') – Method of interpolation. Can be base for scipy.interpolate.interp1d mean or bff for scaling methods and pd``for pandas interpolation methods. Note that the first method is fast and efficient when the number of NaN in the array if relatively few. It is less accurate to use the `base` interpolation when the data is composed of many missing values. Alternatively, the scaled method(the  second one) is proposed to be the alternative way more efficient. Indeed, when ``mean argument is set, function replaces the NaN values by the nonzeros in the raw array and then uses the mean to fit the data. The result of fitting creates a smooth curve where the index of each NaN in the raw array is replaced by its corresponding values in the fit results. The same approach is used for bff method. Conversely, rather than averaging the nonzeros values, it uses the backward and forward strategy to fill the NaN before scaling. mean and bff are more efficient when the data are composed of lot of missing values. When the interpolation method is set to pd, function uses the pandas interpolation but ended the interpolation with forward/backward NaN filling since the interpolation with pandas does not deal with all NaN at the begining or at the end of the array. Default is base.

• fill_value (array-like or (array-like, array_like) or extrapolate, optional) – If a ndarray (or float), this value will be used to fill in for requested points outside of the data range. If not provided, then the default is NaN. The array-like must broadcast properly to the dimensions of the non-interpolation axes. If a two-element tuple, then the first element is used as a fill value for x_new < x[0] and the second element is used for x_new > x[-1]. Anything that is not a 2-element tuple (e.g., list or ndarray, regardless of shape) is taken to be a single array-like argument meant to be used for both bounds as below, above = fill_value, fill_value. Using a two-element tuple or ndarray requires bounds_error=False. Default is extrapolate.

• kws (dict) – Additional keyword arguments from spi.interp1d.

Return type:

array like - New interpoolated array. NaN values are interpolated.

Notes

When interpolated thoughout the complete frequencies i.e all the frequency values using the base method, the missing data in arr can be out of the arr range. So, for consistency and keep all values into the range of frequency, the better idea is to set the param fill_value in kws argument of spi.interp1d to extrapolate. This will avoid an error to raise when the value to interpolated is extra-bound of arr.

References

https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html https://www.askpython.com/python/examples/interpolation-to-fill-missing-entries

Examples

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from watex.utils.exmath  import interpolate1d,
>>> z = np.random.randn(17) *10 # assume 17 freq for 17 values of tensor Z
>>> z [[7, 10, 16]] =np.nan # replace some indexes by NaN values
>>> zit = interpolate1d (z, kind ='linear')
>>> z
... array([ -1.97732415, -16.5883156 ,   8.44484348,   0.24032979,
          8.30863276,   4.76437029, -15.45780568,          nan,
         -4.11301794, -10.94003412,          nan,   9.22228383,
        -15.40298253,  -7.24575491,  -7.15149205, -20.9592011 ,
                 nan]),
>>> zn
...array([ -1.97732415, -16.5883156 ,   8.44484348,   0.24032979,
         8.30863276,   4.76437029, -15.45780568,  -4.11301794,
       -10.94003412,   9.22228383, -15.40298253,  -7.24575491,
        -7.15149205, -20.9592011 , -34.76691014, -48.57461918,
       -62.38232823])
>>> zmean = interpolate1d (z,  method ='mean')
>>> zbff = interpolate1d (z, method ='bff')
>>> zpd = interpolate1d (z,  method ='pd')
>>> plt.plot( np.arange (len(z)),  zit, 'v--',
          np.arange (len(z)), zmean, 'ok-',
          np.arange (len(z)), zbff, '^g:',
          np.arange (len(z)), zpd,'<b:',
          np.arange (len(z)), z,'o',
          )
>>> plt.legend(['interp1d', 'mean strategy', 'bff strategy',
                'pandas strategy', 'data'], loc='best')

watex.utils.exmath.interpolate2d(arr2d, method='slinear', **kws)

Interpolate the data in 2D dimensional array.

If the data contains some missing values. It should be replaced by the interpolated values.

Parameters:

• arr2d (np.ndarray, shape (N, M)) – 2D dimensional data

• method (str, default linear) – Interpolation technique to use. Can be nearest``or ``pad.

• kws (dict) – Additional keywords. Refer to interpolate1d().

Returns:

arr2d – 2D dimensional data interpolated

Return type:

np.ndarray, shape (N, M)

Examples

>>> from watex.methods.em import EM
>>> from watex.utils.exmath import interpolate2d
>>> # make 2d matrix of frequency
>>> emObj = EM().fit(r'data/edis')
>>> freq2d = emObj.make2d (out = 'freq')
>>> freq2d_i = interpolate2d(freq2d )
>>> freq2d.shape
...(55, 3)
>>> freq2d
... array([[7.00000e+04, 7.00000e+04, 7.00000e+04],
       [5.88000e+04, 5.88000e+04, 5.88000e+04],
       ...
        [6.87500e+00, 6.87500e+00, 6.87500e+00],
        [        nan,         nan, 5.62500e+00]])
>>> freq2d_i
... array([[7.000000e+04, 7.000000e+04, 7.000000e+04],
       [5.880000e+04, 5.880000e+04, 5.880000e+04],
       ...
       [6.875000e+00, 6.875000e+00, 6.875000e+00],
       [5.625000e+00, 5.625000e+00, 5.625000e+00]])

References

https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.interpolate.html https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.interpolate.interp2d.html

watex.utils.exmath.invertVES(data=None, rho0=None, h0=None, typeof='HMCMC', **kwargs)

Invert the VES data collected in the exploration area.

Parameters:

• data (DataFrame) – DataFrame containing the depth measurement AB from current electrodes, the potential electrodes MN, and the collected apparent resistivities.

• rho0 (float, optional) – Value of the starting resistivity model. If None, rho0 is set to half the minimum value of the collected apparent resistivity. Units are in Ω.m.

• h0 (float, optional) – Thickness in meters of the first layer. If None, it defaults to 1 meter.

• typeof (str, default 'HMCMC') – Type of inversion scheme. Default is Hybrid Monte Carlo (HMCMC). Another option is Bayesian neural network approach (BNN).

• **kwargs (dict) – Additional keyword arguments for VES data operations. See watex.utils.exmath.vesDataOperator for further details.

Returns:

Inversion results, specifics depend on the inversion scheme used.

Return type:

Tuple[ArrayLike]

Notes

This function is a placeholder for the inversion process. The actual implementation would require detailed algorithms for HMCMC or BNN approaches.

watex.utils.exmath.linkage_matrix(df, columns=None, kind='design', metric='euclidean', method='complete', as_frame=False, optimal_ordering=False)

Compute the distance matrix from the hierachical clustering algorithm

Parameters:

• df (dataframe or NDArray of (n_samples, n_features)) – dataframe of Ndarray. If array is given , must specify the column names to much the array shape 1

• columns (list) – list of labels to name each columns of arrays of (n_samples, n_features) If dataframe is given, don’t need to specify the columns.

• kind (str, ['squareform'|'condense'|'design'], default is {'design'}) – kind of approach to summing up the linkage matrix. Indeed, a condensed distance matrix is a flat array containing the upper triangular of the distance matrix. This is the form that pdist returns. Alternatively, a collection of \(m\) observation vectors in \(n\) dimensions may be passed as an \(m\) by \(n\) array. All elements of the condensed distance matrix must be finite, i.e., no NaNs or infs. Alternatively, we could used the squareform distance matrix to yield different distance values than expected. the design approach uses the complete inpout example matrix also called ‘design matrix’ to lead correct linkage matrix similar to squareform and condense`.

• metric (str or callable, default is {'euclidean'}) – The metric to use when calculating distance between instances in a feature array. If metric is a string, it must be one of the options allowed by sklearn.metrics.pairwise.pairwise_distances(). If X is the distance array itself, use “precomputed” as the metric. Precomputed distance matrices must have 0 along the diagonal.

• method (str, optional, default is {'complete'}) – The linkage algorithm to use. See the Linkage Methods section below for full descriptions.

• optimal_ordering (bool, optional) – If True, the linkage matrix will be reordered so that the distance between successive leaves is minimal. This results in a more intuitive tree structure when the data are visualized. defaults to False, because this algorithm can be slow, particularly on large datasets. See also scipy.cluster.hierarchy.linkage().

Returns:

row_clusters – consist of several rows where each rw represents one merge. The first and second columns denotes the most dissimilar members of each cluster and the third columns reports the distance between those members

Return type:

linkage matrix

watex.utils.exmath.magnitude(cz)

Compute the magnitude of selected conductive zone.

The magnitude parameter is the absolute resistivity value between the minimum \(\min \rho_a\) and maximum \(\max \rho_a\) value of selected anomaly:

\[magnitude=|\min\rho_a -\max\rho_a|\]

Parameters:

cz – array-like. Array of apparent resistivity values composing the conductive zone.

Returns:

Absolute value of anomaly magnitude in ohm.meters.

watex.utils.exmath.moving_average(y, *, window_size=3, method='sma', mode='same', alpha=0.5)

A moving average is used with time series data to smooth out short-term fluctuations and highlight longer-term trends or cycles.

Funtion analyzes data points by creating a series of averages of different subsets of the full data set.

Parameters:

• y (array_like, shape (N,)) – the values of the time history of the signal.

• window_size (int) – the length of the window. Must be greater than 1 and preferably an odd integer number.Default is 3

• method (str) – variant of moving-average. Can be sma, cma, wma and ema for simple, cummulative, weight and exponential moving average. Default is sma.

• mode (str) – returns the convolution at each point of overlap, with an output shape of (N+M-1,). At the end-points of the convolution, the signals do not overlap completely, and boundary effects may be seen. Can be full, same and valid. See ~np.convole for more details. Default is same.

• alpha (float,) – smoothing factor. Only uses in exponential moving-average. Default is .5.

Returns:

ya – Averaged time history of the signal

Return type:

array like, shape (N,)

Notes

The first element of the moving average is obtained by taking the average of the initial fixed subset of the number series. Then the subset is modified by “shifting forward”; that is, excluding the first number of the series and including the next value in the subset.

Examples

>>> import numpy as np ; import matplotlib.pyplot as plt
>>> from watex.utils.exmath  import moving_average
>>> data = np.random.randn (37)
>>> # add gaussion noise to the data
>>> data = 2 * np.sin( data)  + np.random.normal (0, 1 , len(data))
>>> window = 5  # fixed size to 5
>>> sma = moving_average(data, window)
>>> cma = moving_average(data, window, method ='cma' )
>>> wma = moving_average(data, window, method ='wma' )
>>> ema = moving_average(data, window, method ='ema' , alpha =0.6)
>>> x = np.arange(len(data))
>>> plt.plot (x, data, 'o', x, sma , 'ok--', x, cma, 'g-.', x, wma, 'b:')
>>> plt.legend (['data', 'sma', 'cma', 'wma'])

References

watex.utils.exmath.ohmicArea(data=None, search=45.0, sum=False, objective='ohmS', **kws)

Compute the ohmic-area from the Vertical Electrical Sounding data collected in exploration area.

Parameters:

* data: Dataframe pandas - contains the depth measurement AB from current

electrodes, the potentials electrodes MN and the collected apparents resistivities.

* search: float - The depth in meters from which one expects to find a

fracture zone outside of pollutions. Indeed, the search parameter is used to speculate about the expected groundwater in the fractured rocks under the average level of water inrush in a specific area. For instance in Bagoue region , the average depth of water inrush is around 45m. So the search can be specified via the water inrush average value.

* objective: str - Type operation to outputs. By default, the function

outputs the value of pseudo-area in \(\Omega .m^2\). However, for plotting purpose by setting the argument to view, its gives an alternatively outputs of X and Y, recomputed and projected as weel as the X and Y values of the expected fractured zone. Where X is the AB dipole spacing when imaging to the depth and Y is the apparent resistivity computed

kws: dict - Additionnal keywords arguments from |VES| data operations.

See watex.utils.exmath.vesDataOperator() for futher details.

Returns:

List of twice tuples:

• Tuple(ohmS, error, roots):

  • `ohmS`is the pseudo-area computed expected to be a fractured zone

  • error is the integration error

  • roots is the integration boundaries of the expected fractured

    zone where the basement rocks is located above the resistivity transform function. At these points both curves values equal to null.

• Tuple (XY, fit XY,XYohmSarea):

  • XY is the ndarray(nvalues, 2) of the operated of AB dipole

    spacing and resistivity rhoa values.

  • fit XY is the fitting ndarray(nvalues, 2) uses to redraw the

    dummy resistivity transform function.

  • XYohmSarea is ndarray(nvalues, 2) of the dipole spacing and

    resistiviy values of the expected fracture zone.

Raises:

VESError

If the search is greater or equal to the maximum investigation depth in meters.

Notes

The ohmS value calculated from pseudo-area is a fully data-driven parameter and is used to evaluate a pseudo-area of the fracture zone from the depth where the basement rock is supposed to start. Usually, when exploring deeper using the Vertical Electrical Sounding, we are looking for groundwater in thefractured rock that is outside the anthropic pollution (Biemi, 1992). Since the VES is an indirect method, we cannot ascertain whether the presumed fractured rock contains water inside. However, we assume that the fracture zone could exist and should contain groundwater. Mathematically, based on the VES1D model proposed by `Koefoed, O. (1976)`_ , we consider a function \(\rho_T(l)\), a set of reducing resistivity transform function to lower the boundary plane at half the current electrode spacing \((l)\). From the sounding curve \(\rho_T(l)\), curve an imaginary basement rock \(b_r (l)\) of slope equal to 45° with the horizontal \(h(l)\) was created. A pseudo-area \(S(l)\) should be defined by extending from \(h(l)\) the \(b_r (l)\) curve when the sounding curve \(\rho_T(l)\) is below \(b_r(l)\), otherwise \(S(l)\) is equal to null. The computed area is called the ohmic-area \(ohmS\) expressed in \(\Omega .m^2\) and constitutes the expected fractured zone. Thus \(ohmS\) ≠ \(0\) confirms the existence of the fracture zone while of \(Ohms=0\) raises doubts. The equation to determine the parameter is given as:

\[ \begin{align}\begin{aligned}ohmS & = &\int_{ l_i}^{l_{i+1}} S(l)dl \quad {s.t.}\\\begin{split}S(l) & = & b_r (l) - \rho_T (l) \quad \text{if} \quad b_r (l) > \rho_T (l) \\ & = & 0. \quad \text{if} \quad b_r (l) \leq \rho_T (l)\end{split}\\b_r(l) & = & l + h(l) \quad ; \quad h(l) = \beta\\\rho_T(l) & = & l^2 \int_{0}^{\infty} T_i( \lambda ) h_1( \lambda l) \lambda d\lambda\end{aligned}\end{align} \]

where \(l_i \quad \text{and} \quad l_{i+1}\) solve the equation \(S(l=0)\); \(l\) is half the current electrode spacing \(AB/2\), and \(h_1\) denotes the first-order of the Bessel function of the first kind, \(\beta\) is the coordinate value on y-axis direction of the intercept term of the \(b_r(l)\) and \(h(l)\), \(T_i(\lambda )\) resistivity transform function, \(lamda\) denotes the integral variable, where n denotes the number of layers, \(rho_i\) and \(h_i\) are the resistivity and thickness of the \(i-th\) layer, respectively. Get more explanations and cleareance of formula in the paper of `Kouadio et al 2022`_.

. _Cote d’Ivoire: https://en.wikipedia.org/wiki/Ivory_Coast

Examples

>>> from watex.utils.exmath import ohmicArea
>>> from watex.utils.coreutils import vesSelector
>>> data = vesSelector (f= 'data/ves/ves_gbalo.xlsx')
>>> (ohmS, err, roots), *_ = ohmicArea(data = data, search =45, sum =True )
... (13.46012197818152, array([5.8131967e-12]), array([45.        , 98.07307307]))
# pseudo-area is computed between the spacing point AB =[45, 98] depth.
>>> _, (XY.shape, XYfit.shape, XYohms_area.shape) = ohmicArea(
                AB= data.AB, rhoa =data.resistivity, search =45,
                objective ='plot')
... ((26, 2), (1000, 2), (8, 2))

watex.utils.exmath.plotOhmicArea(data=None, search=45.0, pre_computed=False, xy=None, xyf=None, xyarea=None, colors=None, fbtw=False, **plot_kws)

Plot the Vertical Electrical Sounding data ohmic -area

Parameters:

• data (*) – contains the depth measurement AB from current electrodes, the potentials electrodes MN and the collected apparent resistivities.

• search (*) – The depth in meters from which one expects to find a fracture zone outside of pollutions. Indeed, the search parameter is used to speculate about the expected groundwater in the fractured rocks under the average level of water inrush in a specific area. For instance in Bagoue region , the average depth of water inrush is around 45m. So the search can be specified via the water inrush average value.

• pre_computed (bool, default=False,) – If True computed the ohmic_area parameters. If False, the ohmic area arguments must be passed to xy, xyf and xyarea, otherwise an errors will raise.

• xy (array-like of shape (n_AB, 2)) – Arraylike of the sanitized depth measurement AB from current. electrodes n_AB. See vesDataOperator().

• xyf (array-like of shape (n_fit_samples, 2)) – Array-like of the fitted samples i.e the number of points for fitting the sounding resistivity values from the surface thin the total depth. The fitted rhoa showns a smooth curves. The default point is 1000.

• xyarea (array-like of shape (n_area, 2)) – Arraylike of the resistivity positions of the depth measurment AB where the fractured zone is found.

• fbtw (bool, default=False,) – If True, filled the computed fractured zone using the parameters computed from xyf and xyarea.

• kws (dict - Additionnal keywords arguments from Vertical Electrical Sounding data operations.) – See watex.utils.exmath.vesDataOperator() for futher details.

Notes

The first and second columns of xy, xyfit and xyarea are the position AB/2 and their corresponding resistivity values.

Examples

>>> from watex.datasets import load_semien
>>> from watex.utils.exmath import plotOhmicArea
>>> ves_data = load_semien ()
>>> plotOhmicArea (ves_data)

watex.utils.exmath.plot_(*args, fig_size=None, raw=False, style='seaborn', dtype='erp', kind=None, fig_title_kws=None, fbtw=False, fig=None, ax=None, **kws)

Quick visualization for fitting model, Electrical Resistivity Profiling and Vertical Electrical Sounding curves.

param x:

array-like - array of data for x-axis representation

param y:

array-like - array of data for plot y-axis representation

param fig_size:

tuple - Matplotlib (MPL) figure size; should be a tuple value of integers e.g. figsize =(10, 5).

param raw:

bool- Originally the plot_ function is intended for the fitting Electrical Resistivity Profiling model i.e. the correct value of Electrical Resistivity Profiling data. However, when the raw is set to True, it plots the both curves: The fitting model as well as the uncorrected model. So both curves are overlaining or supperposed.

param style:

str - Pyplot style. Default is seaborn

param dtype:

str - Kind of data provided. Can be Electrical Resistivity Profiling data or Vertical Electrical Sounding data. When the Electrical Resistivity Profiling data are provided, the common plot is sufficient to visualize all the data insight i.e. the default value of kind is kept to None. However, when the data collected is Vertical Electrical Sounding data, the convenient plot for visualization is the loglog for parameter kind` while the dtype can be set to VES to specify the labels into the x-axis. The default value of dtype is erp for common visualization.

param kind:

str - Use to specify the kind of data provided. See the explanation of dtype parameters. By default kind is set to None i.e. its keep the normal plots. It can be loglog, semilogx and semilogy.

param fbtw:

bool, default=False, Mostly used for Vertical Electrical Sounding data. If True, filled the computed fractured zone using the parameters computed from ohmicArea().

param fig_title_kws:

dict, Additional keywords argument passed in dictionnary to customize the figure title.

param fig:

Matplotlib.pyplot.figure add plot on the same figure.

param kws:

dict - Additional Matplotlib plot keyword arguments. To cus- tomize the plot, one can provide a dictionnary of MPL keyword additional arguments like the example below.

Example:

>>> import numpy as np
>>> from watex.utils.exmath import plot_
>>> x, y = np.arange(0 , 60, 10) ,np.abs( np.random.randn (6))
>>> KWS = dict (xlabel ='Stations positions', ylabel= 'resistivity(ohm.m)',
            rlabel = 'raw cuve', rotate = 45 )
>>> plot_(x, y, '-ok', raw = True , style = 'seaborn-whitegrid',
          figsize = (7, 7) ,**KWS )

. _Cote d’Ivoire: https://en.wikipedia.org/wiki/Ivory_Coast

watex.utils.exmath.plot_confidence_in(z_or_edis_obj_list, /, tensor='res', view='1d', drop_outliers=True, distance=None, c_line=False, view_ci=True, figsize=(6, 2), fontsize=4.0, dpi=300.0, top_label='Stations', rotate_xlabel=90.0, fbtw=True, savefig=None, **plot_kws)

Plot data confidency from tensor errors.

The default tensor for evaluating the data confidence is the resistivity at TE mode (‘xy’).

Check confidence in the data before starting the concrete processing seems meaningful. In the area with complex terrain, with high topography addition to interference noises, signals are weals or missing especially when using AMT survey. The most common technique to do this is to eliminate the bad frequency and interpolate the remains one. However, the tricks for eliminating frequency differ from one author to another. Here, the tip using the data confidence seems meaningful to indicate which frequencies to eliminate (at which stations/sites) and which ones are still recoverable using the tensor recovering strategy.

The plot implements three levels of confidence:

• High confidence: \(conf. \geq 0.95\) values greater than 95%

• Soft confidence: \(0.5 \leq conf. < 0.95\). The data in this confidence range can be beneficial for tensor recovery to restore the weak and missing signals.

• bad confidence: \(conf. <0.5\). Data in this interval must be deleted.

Parameters:

• z_or_edis_obj_list (list of watex.edi.Edi or watex.externals.z.Z) – A collection of EDI- or Impedances tensors objects.

• tensor (str, default='res') – Tensor name. Can be [ ‘resistivity’|’phase’|’z’|’frequency’]

• view (str, default='1d') – Type of plot. Can be [‘1D’|’2D’]

• drop_outliers (bool, default=True) – Suppress the ouliers in the data if True.

• distance (float, optional) – Distance between stations/sites

• fontsize (float, default=3.) – label font size.

• figsize (Tuple, default=(6, 2)) – Figure size.

• c_line (bool, default=True,) – Display the confidence line in two dimensinal view.

• dpi (int, default=300) – Image resolution in dot-per-inch

• rotate_xlabel (float, default=90.) – Angle to rotate the stations/sites labels

• top_label (str,default='Stations') – Labels the sites either using the survey name.

• view_ci (bool,default=True,) – Show the marker of confidence interval.

• fbtw (bool, default=True,) – Fill between confidence interval.

• plot_kws (dict,) – Additional keywords pass to the plot()

See also

watex.methods.Processing.zrestore

For more details about the function for tensor recovering technique.

Examples

>>> from watex.utils.exmath import plot_confidence_in
>>> from watex.datasets import fetch_data
>>> emobj  = fetch_data ( 'huayuan', samples = 25, clear_cache =True,
                         key='raw').emo
>>> plot_confidence_in (emobj.ediObjs_ ,
                        distance =20 ,
                        view ='2d',
                        figsize =(6, 2)
                        )
>>> plot_confidence_in (emobj.ediObjs_ , distance =20 ,
                        view ='1d', figsize =(6, 3), fontsize =5,
                        )

watex.utils.exmath.plot_sfi(cz, p=None, s=None, dipolelength=None, fig_size=(10, 4), style='classic', **plotkws)

Plot sfi parameter components.

Parameters:

• cz (array-like 1d,) – Selected conductive zone

• p (array-like 1d,) – Station positions of the conductive zone.

• dipolelength (float. If p is not given, it will be set) – automatically using the default value to match the cz size. The default value is 10.

• fig_size (tuple, default=(10, 4)) – Matplotlib (MPL) figure size; should be a tuple value of integers

See also

watex.utils.exmath.sfi

for more details about the sfi parameter computation.

Examples

>>> import numpy as np
>>> from watex.utils.exmath import plot_sfi
>>> rang = np.random.RandomState (42)
>>> condzone = np.abs(rang.randn (7))*1e2
>>> plotkws  = dict (rlabel = 'Selected conductive zone (cz)',
                     color=f'{P().frcolortags.get("fr3")}',
                     )
>>> plot_sfi (condzone, **plotkws)

watex.utils.exmath.power(p)

Compute the power of the selected conductive zone. Anomaly power is closely referred to the width of the conductive zone.

The power parameter implicitly defines the width of the conductive zone and is evaluated from the difference between the abscissa \(X_{LB}\) and the end \(X_{UB}\) points of the selected anomaly:

\[power=|X_{LB} - X_{UB} |\]

Parameters:

p – array-like. Station position of conductive zone.

Returns:

Absolute value of the width of conductive zone in meters.

watex.utils.exmath.qc(z_or_edis_obj_list, /, tol=0.5, *, interpolate_freq=False, return_freq=False, tensor='res', return_data=False, to_log10=False, return_qco=False)

Check the quality control in the collection of Z or EDI objects.

Analyse the data in the EDI collection and return the quality control value. It indicates how percentage are the data to be representative.

Parameters:

• z_or_edis_obj_list (list of watex.edi.Edi or watex.externals.z.Z) – A collection of EDI- or Impedances tensors objects.

• tol (float, default=.5) – the tolerance parameter. The value indicates the rate from which the data can be consider as meaningful. Preferably it should be less than 1 and greater than 0. Default is .5 means 50 %. Analysis becomes soft with higher tol values and severe otherwise.

• interpolate_freq (bool,) – interpolate the valid frequency after removing the frequency which data threshold is under the ``1-tol``% goodness

• return_freq (bool, default=False) – returns the interpolated frequency.

• return_data (bool, default= False,) – returns the valid data from up to 1-tol% goodness.

• tensor (str, default='z') – Tensor name. Can be [ resistivity|phase|z|frequency]. Impedance is used for data quality assessment.

• to_log10 (bool, default=True) – convert the frequency value to log10.

• qco (return) –

  retuns quality control object that wraps all usefull informations after control. The following attributes can be fetched as:

  • rate_: the rate of the quality of the data

  • component_: The selected component where data is selected for analysis By default used either xy or yx.

  • mode_: The EM mode. Either the [‘TE’|’TM’] modes

  • freqs_: The valid frequency in the data selected according to the tol parameters. Note that if interpolate_freq is True, it is used instead.

  • invalid_freqs_: Useless frequency dropped in the data during control

  • data_: Valid tensor data either in TE or TM mode.

Returns:

• return the quality control value and interpolated frequency if

return_freq is set to True otherwise return the only the quality control ratio.

• return the the quality control object.

Return type:

Tuple (float ) or (float, array-like, shape (N, )) or QCo

Examples

>>> import watex as wx
>>> data = wx.fetch_data ('huayuan', samples =20, return_data =True ,
                          key='raw')
>>> r,= wx.qc (data)
r
Out[61]: 0.75
>>> r, = wx.qc (data, tol=.2 )
0.75
>>> r, = wx.qc (data, tol=.1 )

watex.utils.exmath.quickplot(arr, dl=10)

Quick plot to see the anomaly

watex.utils.exmath.rhoa2z(rhoa, phs, freq)

Convert apparent resistivity to impendance tensor z

Parameters:

• rhoa (ndarray, shape (N, M)) – Apparent resistivity in \(\Omega.m\)

• phs (ndarray, shape (N, M)) – Phase in degrees

• freq (array-like , shape (N, )) – Frequency in Hertz

: :return: Impendance tensor; Tensor is a complex number in \(\Omega\). :rtype: ndarray, shape (N, M), dtype = ‘complex’

Example:

>>> import numpy as np
>>> rhoa = np.array([1623.73691735])
>>> phz = np.array([45.])
>>> f = np.array ([1014])
>>> rhoa2z(rhoa, phz, f)
... array([[2.54950976+2.54950976j]])

watex.utils.exmath.rhophi2z(rho, phi, freq)

Convert impedance-style information given in Rho/Phi format into complex valued Z.

Parameters:

• rho (ArrayLike 1D/2D) – Resistivity array in \(\Omega.m\). If array is two-dimensional, it should be 2x2 array (real).

• phi (ArrayLike 1D/2D) – Phase array in degree (\(\degree\)). If array is two-dimensional, it should be 2x2 array (real).

• freq (float, arraylike 1d) – Frequency in Hz

Returns:

Z – Z dimension depends to the inputs array rho and phi.

Return type:

Arraylike 1d or 2d , complex

Examples

>>> import numpy as np
>>> from watex.utils.exmath import rhophi2z
>>> rhophi2z (823 , 25 , 500 )
array([1300.00682824+606.20313966j])
>>> rho = np.array ([[823, 700], [723, 526]] )
>>> phi = np.array ([[45, 50], [90, 180]])
>>> rhophi2z (rho, phi , freq= 500  )
array([[ 1.01427314e+03+1.01427314e+03j,  8.50328081e+02+1.01338154e+03j],
       [ 8.23227764e-14+1.34443297e+03j, -1.14673449e+03+1.40434473e-13j]])
>>> rhophi2z (np.array ( [ 823, 700])  , np.array ([45, 50 ])  , [500, 700] )
array([1014.27313876+1014.27313876j, 1006.12175325+1199.04921402j])
>>> rho  = np.abs (np.random.randn (7, 3 ) * 100 )
>>> phi = np.abs ( np.random.randn (7, 3 ) *180 % 90 )
>>> freq = np.abs ( np.random.randn (7) * 100 )
>>> rhophi2z (rho   , phi  , freq )

watex.utils.exmath.savgol_coeffs(window_length, polyorder, deriv=0, delta=1.0, pos=None, use='conv')

Compute the coefficients for a 1-D Savitzky-Golay FIR filter.

Parameters:

• window_length (int) – The length of the filter window (i.e., the number of coefficients). window_length must be an odd positive integer.

• polyorder (int) – The order of the polynomial used to fit the samples. polyorder must be less than window_length.

• deriv (int, optional) – The order of the derivative to compute. This must be a nonnegative integer. The default is 0, which means to filter the data without differentiating.

• delta (float, optional) – The spacing of the samples to which the filter will be applied. This is only used if deriv > 0.

• pos (int or None, optional) – If pos is not None, it specifies evaluation position within the window. The default is the middle of the window.

• use (str, optional) – Either ‘conv’ or ‘dot’. This argument chooses the order of the coefficients. The default is ‘conv’, which means that the coefficients are ordered to be used in a convolution. With use=’dot’, the order is reversed, so the filter is applied by dotting the coefficients with the data set.

Returns:

coeffs – The filter coefficients.

Return type:

1-D ndarray

References

A. Savitzky, M. J. E. Golay, Smoothing and Differentiation of Data by Simplified Least Squares Procedures. Analytical Chemistry, 1964, 36 (8), pp 1627-1639.

See also

savgol_filter

Examples

>>> from watex.exmath.signal import savgol_coeffs
>>> savgol_coeffs(5, 2)
array([-0.08571429,  0.34285714,  0.48571429,  0.34285714, -0.08571429])
>>> savgol_coeffs(5, 2, deriv=1)
array([ 2.00000000e-01,  1.00000000e-01,  2.07548111e-16, -1.00000000e-01,
       -2.00000000e-01])

Note that use=’dot’ simply reverses the coefficients.

>>> savgol_coeffs(5, 2, pos=3)
array([ 0.25714286,  0.37142857,  0.34285714,  0.17142857, -0.14285714])
>>> savgol_coeffs(5, 2, pos=3, use='dot')
array([-0.14285714,  0.17142857,  0.34285714,  0.37142857,  0.25714286])

x contains data from the parabola x = t**2, sampled at t = -1, 0, 1, 2, 3. c holds the coefficients that will compute the derivative at the last position. When dotted with x the result should be 6.

>>> x = np.array([1, 0, 1, 4, 9])
>>> c = savgol_coeffs(5, 2, pos=4, deriv=1, use='dot')
>>> c.dot(x)
6.0

watex.utils.exmath.savgol_filter(x, window_length, polyorder, deriv=0, delta=1.0, axis=-1, mode='interp', cval=0.0)

Apply a Savitzky-Golay filter to an array.

This is a 1-D filter. If x has dimension greater than 1, axis determines the axis along which the filter is applied.

Parameters:

• x (array_like) – The data to be filtered. If x is not a single or double precision floating point array, it will be converted to type numpy.float64 before filtering.

• window_length (int) – The length of the filter window (i.e., the number of coefficients). window_length must be a positive odd integer. If mode is ‘interp’, window_length must be less than or equal to the size of x.

• polyorder (int) – The order of the polynomial used to fit the samples. polyorder must be less than window_length.

• deriv (int, optional) – The order of the derivative to compute. This must be a nonnegative integer. The default is 0, which means to filter the data without differentiating.

• delta (float, optional) – The spacing of the samples to which the filter will be applied. This is only used if deriv > 0. Default is 1.0.

• axis (int, optional) – The axis of the array x along which the filter is to be applied. Default is -1.

• mode (str, optional) – Must be ‘mirror’, ‘constant’, ‘nearest’, ‘wrap’ or ‘interp’. This determines the type of extension to use for the padded signal to which the filter is applied. When mode is ‘constant’, the padding value is given by cval. See the Notes for more details on ‘mirror’, ‘constant’, ‘wrap’, and ‘nearest’. When the ‘interp’ mode is selected (the default), no extension is used. Instead, a degree polyorder polynomial is fit to the last window_length values of the edges, and this polynomial is used to evaluate the last window_length // 2 output values.

• cval (scalar, optional) – Value to fill past the edges of the input if mode is ‘constant’. Default is 0.0.

Returns:

y – The filtered data.

Return type:

ndarray, same shape as x

See also

savgol_coeffs

Notes

Details on the mode options:

‘mirror’:

Repeats the values at the edges in reverse order. The value closest to the edge is not included.

‘nearest’:

The extension contains the nearest input value.

‘constant’:

The extension contains the value given by the cval argument.

‘wrap’:

The extension contains the values from the other end of the array.

For example, if the input is [1, 2, 3, 4, 5, 6, 7, 8], and window_length is 7, the following shows the extended data for the various mode options (assuming cval is 0):

mode       |   Ext   |         Input          |   Ext
-----------+---------+------------------------+---------
'mirror'   | 4  3  2 | 1  2  3  4  5  6  7  8 | 7  6  5
'nearest'  | 1  1  1 | 1  2  3  4  5  6  7  8 | 8  8  8
'constant' | 0  0  0 | 1  2  3  4  5  6  7  8 | 0  0  0
'wrap'     | 6  7  8 | 1  2  3  4  5  6  7  8 | 1  2  3

New in version 0.14.0.

Examples

>>> from watex.utils.exmath import savgol_filter
>>> np.set_printoptions(precision=2)  # For compact display.
>>> x = np.array([2, 2, 5, 2, 1, 0, 1, 4, 9])

Filter with a window length of 5 and a degree 2 polynomial. Use the defaults for all other parameters.

>>> savgol_filter(x, 5, 2)
array([1.66, 3.17, 3.54, 2.86, 0.66, 0.17, 1.  , 4.  , 9.  ])

Note that the last five values in x are samples of a parabola, so when mode=’interp’ (the default) is used with polyorder=2, the last three values are unchanged. Compare that to, for example, mode=’nearest’:

>>> savgol_filter(x, 5, 2, mode='nearest')
array([1.74, 3.03, 3.54, 2.86, 0.66, 0.17, 1.  , 4.6 , 7.97])

watex.utils.exmath.savitzky_golay1d(y, window_size, order, deriv=0, rate=1, mode='same')

Smooth (and optionally differentiate) data with a Savitzky-Golay filter.

The Savitzky-Golay filter removes high frequency noise from data. It has the advantage of preserving the original shape and features of the signal better than other types of filtering approaches, such as moving averages techniques.

Parameters:

• y (array_like, shape (N,)) – the values of the time history of the signal.

• window_size (int) – the length of the window. Must be an odd integer number.

• order (int) – the order of the polynomial used in the filtering. Must be less then window_size - 1.

• deriv (int) – the order of the derivative to compute (default = 0 means only smoothing)

• mode (str) – mode of the border prepending. Should be valid or same. same is used for prepending or appending the first value of array for smoothing.Default is same.

Returns:

ys – the smoothed signal (or it’s n-th derivative).

Return type:

ndarray, shape (N)

Notes

The Savitzky-Golay is a type of low-pass filter, particularly suited for smoothing noisy data. The main idea behind this approach is to make for each point a least-square fit with a polynomial of high order over a odd-sized window centered at the point.

Examples

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from watex.utils.exmath import savitzky_golay1d
>>> t = np.linspace(-4, 4, 500)
>>> y = np.exp( -t**2 ) + np.random.normal(0, 0.05, t.shape)
>>> ysg = savitzky_golay1d(y, window_size=31, order=4)
>>> plt.plot(t, y, label='Noisy signal')
>>> plt.plot(t, np.exp(-t**2), 'k', lw=1.5, label='Original signal')
>>> plt.plot(t, ysg, 'r', label='Filtered signal')
>>> plt.legend()
>>> plt.show()

References

[1]

A. Savitzky, M. J. E. Golay, Smoothing and Differentiation of Data by Simplified Least Squares Procedures. Analytical Chemistry, 1964, 36 (8), pp 1627-1639.

[2]

Numerical Recipes 3rd Edition: The Art of Scientific Computing W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery Cambridge University Press ISBN-13: 9780521880688

[3]

https://en.wikipedia.org/wiki/Savitzky%E2%80%93Golay_filter#Moving_average

watex.utils.exmath.scalePosition(ydata, xdata=None, func=None, c_order=0, show=False, **kws)

Correct data location or position and return new corrected location

Parameters:

• ydata (array_like, series or dataframe) – The dependent data, a length M array - nominally f(xdata, ...).

• xdata (array_like or object) – The independent variable where the data is measured. Should usually be an M-length sequence or an (k,M)-shaped array for functions with k predictors, but can actually be any object. If None, xdata is generated by default using the length of the given ydata.

• func (callable) – The model function, f(x, ...). It must take the independent variable as the first argument and the parameters to fit as separate remaining arguments. The default func is linear function i.e for f(x)= ax +b. where a is slope and b is the intercept value. Setting your own function for better fitting is recommended.

• c_order (int or str) – The index or the column name if ydata is given as a dataframe to select the right column for scaling.

• show (bool) – Quick visualization of data distribution.

• kws (dict) – Additional keyword argument from scipy.optimize_curvefit parameters. Refer to scipy.optimize.curve_fit.

Returns:

• - ydata - array -like - Data scaled

• - popt - array-like Optimal values for the parameters so that the sum of

• the squared residuals of f(xdata, \*popt) - ydata is minimized.

• - pcov - array like The estimated covariance of popt. The diagonals provide

• the variance of the parameter estimate. To compute one standard deviation

• errors on the parameters use perr = np.sqrt(np.diag(pcov)). How the

• sigma parameter affects the estimated covariance depends on absolute_sigma

• argument, as described above. If the Jacobian matrix at the solution

• doesn’t have a full rank, then ‘lm’ method returns a matrix filled with

• np.inf, on the other hand ‘trf’ and ‘dogbox’ methods use Moore-Penrose

• pseudoinverse to compute the covariance matrix.

Examples

>>> from watex.utils import erpSelector, scalePosition
>>> df = erpSelector('data/erp/l10_gbalo.xlsx')
>>> df.columns
... Index(['station', 'resistivity', 'longitude', 'latitude', 'easting',
'northing'],
dtype='object')
>>> # correcting northing coordinates from easting data
>>> northing_corrected, popt, pcov = scalePosition(ydata =df.northing ,
xdata = df.easting, show=True)
>>> len(df.northing.values) , len(northing_corrected)
... (20, 20)
>>> popt  # by default popt =(slope:a ,intercept: b)
...  array([1.01151734e+00, 2.93731377e+05])
>>> # corrected easting coordinates using the default x.
>>> easting_corrected, *_= scalePosition(ydata =df.easting , show=True)
>>> df.easting.values
... array([790284, 790281, 790277, 790270, 790265, 790260, 790254, 790248,
...       790243, 790237, 790231, 790224, 790218, 790211, 790206, 790200,
...       790194, 790187, 790181, 790175], dtype=int64)
>>> easting_corrected
... array([790288.18571705, 790282.30300999, 790276.42030293, 790270.53759587,
...       790264.6548888 , 790258.77218174, 790252.88947468, 790247.00676762,
...       790241.12406056, 790235.2413535 , 790229.35864644, 790223.47593938,
...       790217.59323232, 790211.71052526, 790205.8278182 , 790199.94511114,
...       790194.06240407, 790188.17969701, 790182.29698995, 790176.41428289])
.. _Bagoue region: https://en.wikipedia.org/wiki/Bagou%C3%A9

. _Cote d’Ivoire: https://en.wikipedia.org/wiki/Ivory_Coast

watex.utils.exmath.scale_positions(x, y, *, is_latlon=False, step=None, use_average_dist=False, utm_zone=None, shift=True, view=False, **kws)

Correct the position coordinates.

By default, it consider x and y as easting/latitude and northing/longitude coordinates respectively. It latitude and longitude are given, specify the parameter is_latlon to True.

Parameters:

• x (ArrayLike 1d,) – One dimensional arrays. x can be consider as the abscissa of the landmark and y as ordinates array.

• y (ArrayLike 1d,) – One dimensional arrays. x can be consider as the abscissa of the landmark and y as ordinates array.

• is_latlon (bool, default=False,) – Convert x and y latitude and longitude coordinates values into UTM before computing the distance. x, y should be considered as easting and northing respectively.

• step (float, Optional) – The positions separation. If not given, the average distance between all positions should be used instead.

• use_average_dist (bool, default=False,) – Use the distance computed between positions for the correction.

• utm_zone (str, Optional (##N or ##S)) – UTM zone in the form of number and North or South hemisphere. For instance ‘10S’ or ‘03N’. Note that if x and y are UTM coordinates, the utm_zone must be provide to accurately correct the positions, otherwise the default value 49R should be used which may lead to less accuracy.

• shift (bool, default=True,) – Shift the coordinates from the units of step. This is the default behavor. If False, the positions are just scaled.

• view (bool, default=True) – Visualize the scaled positions

• kws (dict,) – Keyword arguments passed to get_distance()

Returns:

xx, yy – The arrays of position correction from x and y using the bearing.

Return type:

Arraylike 1d,

See also

watex.utils.exmath.get_bearing

Compute the direction of one point relative to another point.

Examples

>>> from watex.utils.exmath import scale_positions
>>> east = [336698.731, 336714.574, 336730.305]
>>> north = [3143970.128, 3143957.934, 3143945.76]
>>> east_c , north_c= scale_positions (east, north, step =20, view =True  )
>>> east_c , north_c
(array([336686.69198337, 336702.53498337, 336718.26598337]),
 array([3143986.09866306, 3143973.90466306, 3143961.73066306]))

watex.utils.exmath.scaley(y, x=None, deg=None, func=None)

Scaling value using a fitting curve.

Create polyfit function from a specifc data points x to correct y values.

Parameters:

• y – array-like of y-axis. Is the array of value to be scaled.

• x – array-like of x-axis. If x is given, it should be the same length as y, otherwise and error will occurs. Default is None.

• func – callable - The model function, f(x, ...). It must take the independent variable as the first argument and the parameters to fit as separate remaining arguments. func can be a linear function i.e for f(x)= ax +b where a is slope and b is the intercept value. It is recommended according to the y value distribution to set up a custom function for better fitting. If func is given, the deg is not needed.

• deg – polynomial degree. If value is None, it should be computed using the length of extrema (local and/or global) values.

Returns:

• y: array scaled - projected sample values got from f.

• x: new x-axis - new axis x_new generated from the samples.

• linear of polynomial function f

References:

Wikipedia, Curve fitting, https://en.wikipedia.org/wiki/Curve_fitting Wikipedia, Polynomial interpolation, https://en.wikipedia.org/wiki/Polynomial_interpolation

Example:

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from watex.exmath import scale_values
>>> rdn = np.random.RandomState(42)
>>> x0 =10 * rdn.rand(50)
>>> y = 2 * x0  +  rnd.randn(50) -1
>>> plt.scatter(x0, y)
>>> yc, x , f = scale_values(y)
>>> plt.plot(x, y, x, yc)

watex.utils.exmath.select_anomaly(rhoa_array, pos_array=None, auto=True, dipole_length=10.0, **kws)

Select the anomaly value from rhoa_array and find its boundaries if auto is set to True. If auto is False, it’s usefull to provide the anomaly boundaries from station position. Change the argument dipole_length i.e. the distance between measurement electrode is not equal to 10 m else give the pos_array. If the pos_array is given, the dipole_length will be recomputed.

Parameters:

• rhoa_array – The apparent resistivity value of Electrical Resistivity Profiling.

• pos_array (array_like) – The array of station position in meters

• auto –

  bool

  Automaticaly of manual computation to select the best anomaly point. Be sure if auto is set to False to provide the anomaly boundary by setting pos_bounds:

  \[pos_bounds=(90, 130)\]

  where 90 is the pk_min and 130 is the pk_max If pos_bounds is not given an station error will probably occurs from StationError.

• dipole_length (float) – Is the distance between two closest measurement. If the value is known it’s better to provide it and don’t need to provied a pos_array value.

• pos_bounds – Is the tuple value of anomaly boundaries composed of pk_min and pk_max. Please refer to compute_power(). When provided the pos_bounds value, please set the dipole_length to accurate the computation of compute_power().

:type pos_bounds:tuple

Returns:

• rhoa : The app. resistivity value of the selected anomaly

• pk_min and the pk_max: refer to compute_power().

• rhoa_max and rhoa_min: refer to compute_magnitude()

Note:

If the auto param is True, the automatic computation will give at most three best animalies ranking according to the resitivity value.

watex.utils.exmath.sfi(cz, p=None, s=None, dipolelength=None, view=False, raw=False, return_components=False, **plotkws)

Compute the pseudo-fracturing index known as sfi.

The sfi parameter does not indicate the rock fracturing degree in the underground but it is used to speculate about the apparent resistivity dispersion ratio around the cumulated sum of the resistivity values of the selected anomaly. It uses a similar approach of IF parameter proposed by Dieng et al (2004). Furthermore, its threshold is set to \(sqrt{2}\) for symmetrical anomaly characterized by a perfect distribution of resistivity in a homogenous medium. The formula is given by:

\[sfi=\sqrt{(P_a^{*}/P_a )^2+(M_a^{*}/M_a )^2}\]

where \(P_a\) and \(M_a\) are the anomaly power and the magnitude respectively. \(P_a^{*}\) is and \(M_a^{*}\) are the projected power and magnitude of the lower point of the selected anomaly.

Parameters:

• cz (array-like,) – Selected conductive zone

• p (array-like,) – Station positions of the conductive zone.

• dipolelength (float. If p is not given, it will be set) – automatically using the default value to match the cz size. The default value is 10..

• view (bool, default=False,) – Visualize the fitting curve. Default is False.

• raw (bool,default=False,) – Overlaining the fitting curve with the raw curve from cz.

• return_components (bool, default=False,) – If True, it returns the different components used for compute sfi especially for external visualization.

• plotkws (dict) – Matplotlib plot keyword arguments.

Returns:

sfi – value computed for pseudo-fracturing index

Return type:

float

Examples

>>> import numpy as np
>>> from watex.property import P
>>> from watex.utils.exmath import sfi
>>> rang = np.random.RandomState (42)
>>> condzone = np.abs(rang.randn (7))
>>> # no visualization and default value `s` with global minimal rho
>>> pfi = sfi (condzone)
... 3.35110143
>>> # visualize fitting curve
>>> plotkws  = dict (rlabel = 'Conductive zone (cz)',
                     label = 'fitting model',
                     color=f'{P().frcolortags.get("fr3")}',
                     )
>>> sfi (condzone, view= True , s= 5, figsize =(7, 7),
          **plotkws )
Out[598]: (array([ 0., 10., 20., 30.]), 1)

References

• See Numpy Polyfit

• See Stackoverflow

  the answer of AkaRem edited by Tobu and Migilson.

• See Numpy Errorstate and

  how to implement the context manager.

watex.utils.exmath.shape(cz, s=Ellipsis, p=Ellipsis)

Compute the shape of anomaly.

The shape parameter is mostly used in the basement medium to depict the better conductive zone for the drilling location. According to Sombo et al. (2011; 2012), various shapes of anomalies can be described such as:

“V”, “U”, “W”, “M”, “K”, “C”, and “H”

The shape consists to feed the algorithm with the Electrical Resistivity Profiling resistivity values by specifying the station \((S_{VES})\). Indeed, mostly, \(S_{VES}\) is the station with a very low resistivity value expected to be the drilling location.

Parameters:

• cz – array-like - Conductive zone resistivity values

• s – int, str - Station position index or name.

• p – Array-like - Should be the position of the conductive zone.

Note

If s is given, p should be provided. If p is missing an error will raises.

Returns:

str - the shape of anomaly.

Example:

>>> import numpy as np
>>> rang = np.random.RandomState(42)
>>> from watex.utils.exmath import shape
>>> test_array1 = np.arange(10)
>>> shape (test_array1)
...  'C'
>>> test_array2 = rang.randn (7)
>>> shape(test_array2)
... 'H'
>>> test_array3 = np.power(10, test_array2 , dtype =np.float32)
>>> shape (test_array3)
... 'H'   # does not change whatever the resistivity values.

References

Sombo, P. A., Williams, F., Loukou, K. N., & Kouassi, E. G. (2011).

Contribution de la Prospection Électrique à L’identification et à la Caractérisation des Aquifères de Socle du Département de Sikensi (Sud de la Côte d’Ivoire). European Journal of Scientific Research, 64(2), 206–219.

Sombo, P. A. (2012). Application des methodes de resistivites electriques

dans la determination et la caracterisation des aquiferes de socle en Cote d’Ivoire. Cas des departements de Sikensi et de Tiassale (Sud de la Cote d’Ivoire). Universite Felix Houphouet Boigny.

watex.utils.exmath.shortPlot(erp, cz=None)

Quick plot to visualize the sample of ERP data overlained to the selected conductive zone if given.

Parameters:

• erp – array_like, the electrical profiling array

• cz – array_like, the selected conductive zone. If None, cz should be plotted.

Example:

>>> import numpy as np
>>> from watex.utils.exmath import shortPlot, define_conductive_zone
>>> test_array = np.random.randn (10)
>>> selected_cz ,*_ = define_conductive_zone(test_array, 7)
>>> shortPlot(test_array, selected_cz )

watex.utils.exmath.smooth1d(ar, /, drop_outliers=True, ma=True, absolute=False, interpolate=False, view=False, x=None, xlabel=None, ylabel=None, fig_size=(10, 5))

Smooth one-dimensional array.

Parameters:

• ar (ArrayLike 1d) – Array of one-dimensional

• drop_outliers (bool, default=True) – Remove the outliers in the data before smoothing

• ma (bool, default=True,) – Use the moving average for smoothing array value. This seems more realistic.

• interpolate (bool, default=False) –

  Interpolate value to fit the original data size after NaN filling.

  New in version 0.2.8.

• absolute (bool, default=False,) – keep postive the extrapolated scaled values. Indeed, when scaling data, negative value can be appear due to the polyfit function. to absolute this value, set absolute=True. Note that converting to values to positive must be considered as the last option when values in the array must be positive.

• view (bool, default =False) – Display curves

• x (ArrayLike, optional) – Abscissa array for visualization. If given, it must be consistent with the given array ar. Raises error otherwise.

• xlabel (str, optional) – Label of x

• ylabel (str, optional) – label of y

• fig_size (tuple , default=(10, 5)) – Matplotlib figure size

Returns:

yc – Smoothed array value.

Return type:

ArrayLike

Examples

>>> import numpy as np
>>> from watex.utils.exmath import smooth1d
>>> # add Guassian Noise
>>> np.random.seed (42)
>>> ar = np.random.randn (20 ) * 20 + np.random.normal ( 20 )
>>> ar [:7 ]
array([6.42891445e+00, 3.75072493e-02, 1.82905357e+01, 2.92957265e+01,
       6.20589038e+01, 2.26399535e+01, 1.12596434e+01])
>>> arc = smooth1d (ar, view =True , ma =False )
>>> arc [:7 ]
array([12.08603102, 15.29819907, 18.017749  , 20.27968322, 22.11900412,
       23.5707141 , 24.66981557])
>>> arc = smooth1d (ar, view =True )# ma=True by default
array([ 5.0071604 ,  5.90839339,  9.6264018 , 13.94679804, 17.67369252,
       20.34922943, 22.00836725])

watex.utils.exmath.smoothing(ar, /, drop_outliers=True, ma=True, absolute=False, interpolate=False, axis=0, view=False, fig_size=(7, 7), xlabel=None, ylabel=None, cmap='binary')

Smooth data along axis.

Parameters:

• ar (ArrayLike 1d or 2d) – One dimensional or two dimensional array.

• drop_outliers (bool, default=True) – Remove the outliers in the data before smoothing along the given axis

• ma (bool, default=True,) – Use the moving average for smoothing array value along axis. This seems more realistic rather than using only the scaling method.

• absolute (bool, default=False,) – keep positive the extrapolated scaled values. Indeed, when scaling data, negative value can be appear due to the polyfit function. to absolute this value, set absolute=True. Note that converting to values to positive must be considered as the last option when values in the array must be positive.

• axis (int, default=0) – Axis along with the data must be smoothed. The default is the along the row.

• view (bool, default =False) – Visualize the two dimensional raw and smoothing grid.

• xlabel (str, optional) – Label of x

• ylabel (str, optional) – label of y

• fig_size (tuple , default=(7, 5)) – Matplotlib figure size

• cmap (str, default='binary') – Matplotlib.colormap to manage the view color

Returns:

arr0 – Smoothed array value.

Return type:

ArrayLike

Examples

>>> import numpy as np
>>> from watex.utils.exmath import smoothing
>>> # add Guassian Noises
>>> np.random.seed (42)
>>> ar = np.random.randn (20, 7 ) * 20 + np.random.normal ( 20, 7 )
>>> ar [:3, :3 ]
array([[ 31.5265026 ,  18.82693352,  34.5459903 ],
       [ 36.94091413,  12.20273182,  32.44342041],
       [-12.90613711,  10.34646896,   1.33559714]])
>>> arc = smoothing (ar, view =True , ma =False )
>>> arc [:3, :3 ]
array([[32.20356863, 17.18624398, 41.22258603],
       [33.46353806, 15.56839464, 19.20963317],
       [23.22466498, 13.8985316 ,  5.04748584]])
>>> arcma = smoothing (ar, view =True )# ma=True by default
>>> arcma [:3, :3 ]
array([[23.96547827,  8.48064226, 31.81490918],
       [26.21374675, 13.33233065, 12.29345026],
       [22.60143346, 16.77242118,  2.07931194]])
>>> arcma_1 = smoothing (ar, view =True, axis =1 )
>>> arcma_1 [:3, :3 ]
array([[18.74017857, 26.91532187, 32.02914421],
       [18.4056216 , 21.81293014, 21.98535213],
       [-1.44359989,  3.49228057,  7.51734762]])

watex.utils.exmath.torres_verdin_filter(arr, /, weight_factor=0.1, beta=1.0, logify=False, axis=Ellipsis)

Calculates the adaptive moving average of a given data array from Torres and Verdin algorithm [1].

Parameters:

• arr (Arraylike 1d) – List or array-like of data points. If two-dimensional array is passed, axis must be specified to apply the filter onto.

• weight_factor (float, default=.1) – Base smoothing factor for window size which gets adjusted by a factor dependent on the rate of change in the data.

• beta (float, default =1.) – Scaling factor to adjust weight_factor during high volatility. It controls how much the weight_factor is adjusted during periods of high volatility.

• logify (bool, default=False,) – By default , Torres uses exponential moving average. So if the values can be logarithmized to ensure the weight be ranged between 0 and 1. This is important when data are resistivity or phase.

• axis (int, default=0) – Axis along which to apply the AMA filter.

Returns:

ama

Return type:

Adaptive moving average

References

[1]

Torres-Verdin and Bostick, 1992, Principles of spatial surface electric field filtering in magnetotellurics: electromagnetic array profiling (EMAP), Geophysics, v57, p603-622.https://doi.org/10.1190/1.2400625

Example

>>> import matplotlib.pyplot as plt
>>> from watex.utils.exmath import torres_verdin_filter
>>> data = np.random.randn(100)
>>> ama = torres_verdin_filter(data)
>>> plt.plot (range (len(data)), data, 'k', range(len(data)), ama, '-or')
>>> # apply on two dimensional array
>>> data2d = np.random.randn(7, 10)
>>> ama2d = torres_verdin_filter ( data2d, axis =0)
>>> fig, ax  = plt.subplots (nrows = 1, ncols = 2 , sharey= True,
                         figsize = (7,7) )
>>> ax[0].imshow(data2d , label ='Raw data', cmap = 'binary' )
>>> ax[1].imshow (ama2d,  label = 'AMA data', cmap ='binary' )
>>> ax[0].set_title ('Raw data')
>>> ax[1].set_title ('AMA data')
>>> plt.legend
>>> plt.show ()

watex.utils.exmath.type_(erp)

Compute the type of anomaly.

The type parameter is defined by the African Hydraulic Study Committee report (CIEH, 2001). Later it was implemented by authors such as (Adam et al., 2020; Michel et al., 2013; Nikiema, 2012). Type comes to help the differenciation of two or several anomalies with the same shape. For instance, two anomalies with the same shape W will differ from the order of priority of their types. The type depends on the lateral resistivity distribution of underground (resulting from the pace of the apparent resistivity curve) along with the whole Electrical Resistivity Profiling survey line. Indeed, four types of anomalies were emphasized:

“EC”, “CB2P”, “NC” and “CP”.

For more details refers to references.

Parameters:

erp – array-like - Array of Electrical Resistivity Profiling line composed of apparent resistivity values.

Returns:

str -The type of anomaly.

Example:

>>> import numpy as np
>>> from watex.utils.exmath import type_
>>> rang = np.random.RandomState(42)
>>> test_array2 = rang.randn (7)
>>> type_(np.abs(test_array2))
... 'EC'
>>> long_array = np.abs (rang.randn(71))
>>> type(long_array)
... 'PC'

References

Adam, B. M., Abubakar, A. H., Dalibi, J. H., Khalil Mustapha,M., & Abubakar,

A. H. (2020). Assessment of Gaseous Emissions and Socio-Economic Impacts From Diesel Generators used in GSM BTS in Kano Metropolis. African Journal of Earth and Environmental Sciences, 2(1),517–523. https://doi.org/10.11113/ajees.v3.n1.104

CIEH. (2001). L’utilisation des méthodes géophysiques pour la recherche

d’eaux dans les aquifères discontinus. Série Hydrogéologie, 169.

Michel, K. A., Drissa, C., Blaise, K. Y., & Jean, B. (2013). Application

de méthodes géophysiques à l ’ étude de la productivité des forages d ’eau en milieu cristallin : cas de la région de Toumodi ( Centre de la Côte d ’Ivoire). International Journal of Innovation and Applied Studies, 2(3), 324–334.

Nikiema, D. G. C. (2012). Essai d‘optimisation de l’implantation géophysique

des forages en zone de socle : Cas de la province de Séno, Nord Est du Burkina Faso (IRD). (I. / I. Ile-de-France, Ed.). IST / IRD Ile-de-France, Ouagadougou, Burkina Faso, West-africa. Retrieved from http://documentation.2ie-edu.org/cdi2ie/opac_css/doc_num.php?explnum_id=148

watex.utils.exmath.vesDataOperator(AB=None, rhoa=None, data=None, typeofop='mean', outdf=False)

Process VES data to handle duplicated spacing distances (AB) by applying specified operations to the corresponding resistivity values (rhoa).

In VES measurements, it’s common to encounter duplicated AB values with different resistivity readings due to minor adjustments in the MN distance. This function consolidates such duplicates by either averaging, taking the median, or randomly selecting one of the resistivity readings.

Parameters:

• AB (Optional[ArrayLike], optional) – 1D array of AB spacings from current electrodes, representing the exploration depth measurements in meters. If data is provided, this parameter is ignored.

• rhoa (Optional[ArrayLike], optional) – 1D array of apparent resistivity values corresponding to AB spacings, measured in ohm-meters (Ω·m). If data is provided, this parameter is ignored.

• data (Optional[DataFrame], optional) – DataFrame containing both AB spacings and corresponding rhoa values. Overrides AB and rhoa parameters when provided.

• typeofop (str, optional) –

  Specifies the operation to apply to rhoa values for duplicated AB spacings: - ‘mean’: Calculates the mean of rhoa values for each unique AB spacing. - ‘median’: Determines the median of rhoa values for each unique AB spacing. - ‘leaveoneout’: Randomly selects one rhoa value from the duplicates for each

  unique AB spacing. This approach is useful for experiments with significant measurement variance at the same AB spacing.

  Default is ‘mean’.

• outdf (bool, optional) – Determines the format of the output. If True, returns a DataFrame with processed AB and rhoa values; otherwise, returns a tuple (AB, rhoa). Default is False.

Returns:

Processed AB and rhoa values. The format is dictated by the outdf parameter: a DataFrame if True, or a tuple (AB, rhoa) if False.

Return type:

Union[Tuple[ArrayLike, ArrayLike], DataFrame]

Examples

>>> from watex.utils.exmath import vesDataOperator
>>> AB = np.array([10, 10, 20, 30, 30])
>>> rhoa = np.array([100, 105, 150, 200, 195])
>>> # Processing with mean operation
>>> AB_proc, rhoa_proc = vesDataOperator(AB, rhoa, typeofop='mean')
>>> print(AB_proc)
[10 20 30]
>>> print(rhoa_proc)
[102.5 150 197.5]

>>> # Using DataFrame input and median operation
>>> data = pd.DataFrame({'AB': AB, 'rhoa': rhoa})
>>> df_proc = vesDataOperator(data=data, typeofop='median', outdf=True)
>>> print(df_proc)
     AB  rhoa
0  10  102.5
1  20  150.0
2  30  197.5

watex.utils.exmath.vesDataOperator2(AB=None, rhoa=None, data=None, typeofop=None, outdf=False)

Check the data in the given deep measurement and set the suitable operations for duplicated spacing distance of current electrodes AB.

Sometimes at the potential electrodes (MN), the measurement of AB are collected twice after modifying the distance of MN a bit. At this point, two or many resistivity values are targetted to the same distance AB (AB still remains unchangeable while while MN is changed). So the operation consists whether to average (mean) the resistiviy values or to take the median values or to leaveOneOut (i.e. keep one value of resistivity among the different values collected at the same point`AB`) at the same spacing AB. Note that for the LeaveOneOut`, the selected resistivity value is randomly chosen.

Parameters:

• AB (array-like 1d,) – Spacing of the current electrodes when exploring in deeper. Is the depth measurement (AB/2) using the current electrodes AB. Units are in meters.

• rhoa (array-like 1d) – Apparent resistivity values collected by imaging in depth. Units are in \(\Omega {.m}\) not \(log10(\Omega {.m})\)

• data (DataFrame,) – It is composed of spacing values AB and the apparent resistivity values rhoa. If data is given, params AB and rhoa should be kept to None.

• typeofop (str,['mean'|'median'|'leaveoneout'], default='mean') – Type of operation to apply to the resistivity values rhoa of the duplicated spacing points AB. The default operation is mean.

• outdf (bool , default=False,) – Outpout a new dataframe composed of AB and rhoa; data renewed.

Returns:

• - Tuple of (AB, rhoa) (New values computed from typeofop)

• - DataFrame (New dataframe outputed only if outdf is True.)

Notes

By convention AB and MN are half-space dipole length which correspond to AB/2 and MN/2 respectively.

Examples

>>> from watex.utils.exmath import vesDataOperator
>>> from watex.utils.coreutils import vesSelector
>>> data = vesSelector ('data/ves/ves_gbalo.xlsx')
>>> len(data)
... (32, 3) # include the potentiel electrode values `MN`
>>> df= vesDataOperator(data.AB, data.resistivity,
                        typeofop='leaveOneOut', outdf =True)
>>> df.shape
... (26, 2) # exclude `MN` values and reduce(-6) the duplicated values.

watex.utils.exmath.z2rhoa(z, freq)

Convert impendance tensor z to apparent resistivity

Parameters:

• z (ndarray, shape (N, M)) – Impedance tensor in \(\Omega\)

• freq (array-like , shape (N, )) – Frequency in Hertz

: :return: Apparent resistivity in \(\Omega.m\) :rtype: ndarray, shape (N, M)

Example:

>>> import numpy as np
>>> z = np.array([2 + 1j *3 ])
>>> f = np.array ([1014])
>>> z2rhoa(z, f)
... array([[1623.73691735]])