Vizualizacija Vertikalnih i Horizontalnih Preseka Geoelektrične Otpornosti – GeoElektrična tomografija

Geoelektrična tomografija (engl. Electrical Resistivity Tomography, ERT) je geofizička metoda koja se koristi za istraživanje podzemne strukture Zemlje pomoću merenja električne otpornosti tla. Tehnika je neinvazivna, što znači da ne zahteva kopanje ili bušenje, i veoma je korisna u geologiji, inženjerskoj geofizici, hidrogeologiji, kao i u arheologiji i istraživanju zagađenja tla.

Princip rada

Geoelektrična tomografija se zasniva na Ohmovom zakonu i meri kako se električni tok (struja) ponaša kada prolazi kroz različite vrste podzemnih materijala.

Elektrode se postavljaju u liniji (ili mreži) na površini terena.
Kroz par elektroda se ubacuje električna struja, dok se sa drugih elektroda meri električni potencijal (napon).
Na osnovu razlike u naponu i poznate struje računa se tzv. prividna otpornost.
Ponavljanjem merenja za različite kombinacije elektroda dobija se dvodimenzionalna (ili trodimenzionalna) slika podzemnih struktura.

Zašto otpornost?

Otpornost zemljišta zavisi od:

Mineralnog sastava
Vlažnosti i zasićenosti vodom
Poroznosti
Temperature
Prisutnosti gline, peska, stena, pukotina itd.

Na primer:

Suva stena → visoka otpornost
Vlažna glina → niska otpornost

Rezultat

Rezultat geoelektrične tomografije je mapa podzemne otpornosti u vidu preseka (npr. kontura ili lažno obojenih mapa), gde različite boje predstavljaju različite vrednosti otpornosti.

Primena

Hidrogeologija – otkrivanje podzemnih voda, debljine vodonosnih slojeva
Inženjerska geologija – identifikacija nestabilnih zona, klizišta
Ekologija – lociranje zagađenja tla
Arheologija – otkrivanje zidova, grobnica i drugih podzemnih objekata
Rudarska istraživanja – pronalaženje rudnih tela

Tehnike elektroda

Wenner
Schlumberger
Dipol-dipol
Pole-dipol
Gradient
(Svaka konfiguracija ima svoje prednosti u dubini prodiranja i rezoluciji.)

Softver i analiza

Nakon terenskog merenja, podaci se analiziraju pomoću računarskih programa koji vrše inverziju podataka, tj. pretvaraju merenja u verovatnu raspodelu otpornosti u dubinu.

Programu za geoelektričnu tomografiju (EGT) i učitavanje .DAT fajlova — evo detaljnog opisa formata .DAT fajla, za šta se koristi, i primera sa objašnjenjem.

Šta je `.DAT` fajl u kontekstu ERT/EGT?

U geoelektričnim merenjima .DAT fajl sadrži sirove podatke merenja dobijene iz terena: raspored elektroda, rastojanja, struju, napon i izračunatu prividnu otpornost. Fajl je često u CSV ili prostorno formatiranom tekstualnom obliku, i koristi se za:

dalju analizu (npr. u Pythonu, Res2DInv, ZondRes2D, BERT, AGI EarthImager)
vizualizaciju kao konture, preseke, 2D ili 3D modele
izvođenje inverzije (rekonstrukcija stvarne strukture tla)

Tipična struktura ERT `.DAT` fajla

Evo primera iz tvoje ranije poruke:

dat,x1,y1,x2,y2,x3,y3,x4,y4,Rho
001.dat,0,0,1,0,2,0,3,0,25.4
001.dat,1,0,2,0,3,0,4,0,23.1
001.dat,2,0,3,0,4,0,5,0,18.9
...

Objašnjenje kolona:

Kolona	Opis
`dat`	Naziv fajla / identifikator profila
`x1, y1`	Pozicija prve elektrode (A – injekcija struje)
`x2, y2`	Pozicija druge elektrode (B – injekcija struje)
`x3, y3`	Pozicija treće elektrode (M – merenje napona)
`x4, y4`	Pozicija četvrte elektrode (N – merenje napona)
`Rho`	Prividna električna otpornost (u Ohm·m)

Za koje programe se koristi ovaj format?

Kompatibilan sa:

Res2DInv / Res3DInv (uz konverziju)
ZondRes2D
EarthImager
MATLAB obrada
Surfer (za mapiranje kontura)

Program vert.py je Python alat za vizuelizaciju geoelektričnih vertikalnih profila na osnovu merenja prividne otpornosti. Njegova osnovna namena je da automatski učita .dat ili .csv fajlove koji sadrže terenske podatke, interpolira ih u raster i nacrta kontur mapu (izolinije) otpornosti u 2D preseku dubina.

Namena programa `vert.py`

✔️ Učitavanje geoelektričnih merenja iz .dat ili .csv fajlova
✔️ Interpolacija merenih tačaka u regularnu 2D mrežu (x, z)
✔️ Prikaz i čuvanje kontur slike podzemne raspodele otpornosti
✔️ Vizualizacija vertikalnog preseka podzemne strukture

Šta program radi korak po korak

1. Učitavanje fajlova

Pretražuje sve .dat fajlove u trenutnom folderu
Za svaki fajl poziva load_data, koji razlikuje .dat i .csv:
- .dat: ručno se parsira svaka linija (formata: x, z, rho)
- .csv: koristi pandas.read_csv

2. Pretvaranje tačaka u mrežu

Pomoću funkcije interpolate_to_grid koristi se scipy.interpolate.griddata da se:
- iz nepravilno raspoređenih mernih tačaka
- napravi regularna 2D mreža (Xi, Zi) sa interpolisanim vrednostima otpornosti Ri

3. Crtanje kontura

Funkcija save_contour_as_image:
- Pravi kontur mapu sa 15 nivoa boja (plt.contourf)
- Preklapa izolinije (plt.contour)
- Dodaje bočnu skalu (colorbar), obrće osu z (dubina ide na dole)
- Čuva sliku u formatu ime_fajla_konture.jpg (300 DPI)

Upotreba

Postavi vert.py u folder sa .dat fajlovima
Pokreni komandu:

python3 vert.py

Program automatski:
- učitava sve .dat fajlove
- interpolira podatke
- generiše kontur slike npr. profil1_konture.jpg, profil2_konture.jpg, itd.

Ključne prednosti

Automatizovano za više fajlova
Radi sa nepravilnim tačkama i interpolira u raster
Konture se automatski obrću vertikalno (dubina ide nadole)
Rezultat su jasne, publikacione slike preseka

Da bi program vert.py radio ispravno, moraš imati instalirane sledeće Python biblioteke (zavisnosti):

pip3 install pandas numpy matplotlib scipy

Programski kod za vert.py

# vert.py
# The MIT License (MIT)
# Copyright (c) 2025 Aleksandar Maričić
#
# Ovim se omogućava bilo kome da koristi, kopira, menja, spaja, objavljuje,
# distribuira, daje podlicencu i/ili prodaje kopije ovog softverskog programa,
# uz uslov da u svim kopijama ili značajnim delovima softverskog programa bude
# uključena sledeća obavest:
#
# Copyright (c) 2025 Aleksandar Maričić
#
# Ovaj softverski program je pružen "takav kakav jeste", bez bilo kakvih garancija,
# izričitih ili impliciranih, uključujući, ali ne ograničavajući se na, garancije o
# prikladnosti za prodaju ili pogodnosti za određenu svrhu. U svakom slučaju, autori
# ili nosioci prava nisu odgovorni za bilo kakvu štetu ili druge obaveze koje mogu nastati
# usled upotrebe ovog softverskog programa.
# Program vert.py je Python alat za vizuelizaciju geoelektričnih vertikalnih profila na osnovu 
# merenja prividne otpornosti. Njegova osnovna namena je da automatski učita .dat ili .csv fajlove 
# koji sadrže terenske podatke, interpolira ih u raster i nacrta kontur mapu (izolinije) otpornosti u 2D preseku dubina.

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import griddata
import os
import sys
import glob

def konvertuj_dat_u_df(ulazni_fajl):
    with open(ulazni_fajl, 'r', encoding='utf-8') as f:
        linije = f.readlines()

    podatak_pocetak = 5
    podaci = []
    for linija in linije[podatak_pocetak:]:
        linija = linija.strip()
        if linija == '':
            continue
        try:
            x, z, rho = map(lambda s: float(s.replace(',', '.')), linija.split(','))
            podaci.append([x, z, rho])
        except ValueError:
            print(f"[UPOZORENJE] Preskačem liniju: {linija}")
            continue

    df = pd.DataFrame(podaci, columns=['x', 'z', 'ρa (Ohm·m)'])
    print(f"[INFO] Učitano {len(df)} redova iz '{ulazni_fajl}'")
    return df

def load_data(filename):
    ext = os.path.splitext(filename)[1].lower()
    if ext == '.dat':
        return konvertuj_dat_u_df(filename)
    elif ext == '.csv':
        print(f"[INFO] Učitavanje CSV fajla '{filename}'...")
        return pd.read_csv(filename)
    else:
        raise ValueError(f"Nepoznat format fajla: {ext}")

def interpolate_to_grid(df, res=0.1, method='cubic'):
    x, z, rho = df['x'], df['z'], df['ρa (Ohm·m)']
    xi = np.arange(min(x), max(x), res)
    zi = np.arange(min(z), max(z), res)
    Xi, Zi = np.meshgrid(xi, zi)
    Ri = griddata((x, z), rho, (Xi, Zi), method=method)
    return Xi, Zi, Ri

def save_contour_as_image(Xi, Zi, Ri, output_filename='ert_konture.jpg'):
    plt.figure(figsize=(10, 5))
    levels = np.linspace(np.nanmin(Ri), np.nanmax(Ri), 15)
    contour = plt.contourf(Xi, Zi, Ri, levels=levels, cmap='plasma')
    plt.contour(Xi, Zi, Ri, levels=levels, colors='black', linewidths=0.5)
    plt.colorbar(contour, label='ρa (Ohm·m)')
    plt.gca().invert_yaxis()
    plt.gca().set_aspect('equal')
    plt.xlabel('Položaj (m)')
    plt.ylabel('Dubina (m)')
    plt.title('Izolinije prividne otpornosti (ERT konture)')
    plt.grid(True)
    plt.tight_layout()
    plt.savefig(output_filename, dpi=300)
    plt.close()
    print(f"[INFO] Sačuvano kao '{output_filename}'")

def main():
    folder = '.'  # trenutni folder
    dat_files = glob.glob(os.path.join(folder, '*.dat')) + glob.glob(os.path.join(folder, '*.DAT'))
    if not dat_files:
        print("[GREŠKA] Nema .DAT fajlova u folderu.")
        sys.exit(1)

    for dat_file in sorted(dat_files):
        print(f"\n[INFO] Procesiram fajl: {dat_file}")
        df = load_data(dat_file)
        Xi, Zi, Ri = interpolate_to_grid(df, res=0.1)
        output_image = os.path.splitext(dat_file)[0].lower() + '_konture.jpg'
        save_contour_as_image(Xi, Zi, Ri, output_image)

if __name__ == '__main__':
    main()

Program vert2horDAT.py automatski konvertuje više vertikalnih geoelektričnih profila (ERT preseka) u horizontalne preseke na zadatim dubinama, i rezultate snima u .dat fajlovima pogodnim za dalju vizualizaciju i analizu.

Upotreba

Stavi sve fajlove KON19A40.DAT, KON20A40.DAT, …, KON28A40.DAT u isti folder sa vert2horDAT.py
Pokreni u terminalu:

python3 vert2horDAT.py

Program će generisati .dat fajl za svaki horizontalni presek

Upotreba i značaj

Omogućava horizontalnu analizu geoelektričnih slojeva
Idealno za detekciju horizontalnih anomalija (npr. podzemne vode, objekti, barijere)
Priprema podatke za dalje prikaze (npr. konture, 3D modele)
Može se koristiti kao ulaz za surfer, matplotlib, ili druge vizualizacione alate

Zavisnosti

Program koristi sledeće Python biblioteke:

pip3 install numpy pandas scipy

Programski kod za vert2horDAT.py

# The MIT License (MIT)
# Copyright (c) 2025 Aleksandar Maričić
#
# Ovim se omogućava bilo kome da koristi, kopira, menja, spaja, objavljuje,
# distribuira, daje podlicencu i/ili prodaje kopije ovog softverskog programa,
# uz uslov da u svim kopijama ili značajnim delovima softverskog programa bude
# uključena sledeća obavest:
#
# Copyright (c) 2025 Aleksandar Maričić
#
# Ovaj softverski program je pružen "takav kakav jeste", bez bilo kakvih garancija,
# izričitih ili impliciranih, uključujući, ali ne ograničavajući se na, garancije o
# prikladnosti za prodaju ili pogodnosti za određenu svrhu. U svakom slučaju, autori
# ili nosioci prava nisu odgovorni za bilo kakvu štetu ili druge obaveze koje mogu nastati
# usled upotrebe ovog softverskog programa.
# Program vert2horDAT.py automatski konvertuje više vertikalnih geoelektričnih profila (ERT preseka)
# u horizontalne preseke na zadatim dubinama, i rezultate snima u .dat 
# fajlovima pogodnim za dalju vizualizaciju i analizu.


import os
import numpy as np
import pandas as pd
from scipy.interpolate import interp1d

folder = '.'  # folder sa KONxxA40.DAT fajlovima
brojevi_profila = range(19, 29)  # KON19 do KON28
razmak_profila = 1.0  # udaljenost profila u metrima (po y)

# Generišemo dubine u koracima od 0.5m od 1m do 6m uključivo
#dubine_za_prekseke = np.arange(0.5, 5.5 + 0.5, 0.5)
dubine_za_prekseke = np.arange(0.5, 5.5 + 0.25, 0.25)
def konvertuj_dat_u_df(putanja):
    with open(putanja, 'r', encoding='utf-8') as f:
        linije = f.readlines()
    podaci = []
    start = 5
    for linija in linije[start:]:
        linija = linija.strip()
        if linija == '' or linija.startswith('0'):
            continue
        try:
            x, z, rho = map(lambda s: float(s.replace(',', '.')), linija.split(','))
            podaci.append([x, z, rho])
        except Exception:
            continue
    df = pd.DataFrame(podaci, columns=['x', 'z', 'ρa'])
    return df

def ucitaj_sve_profili():
    profili = []
    for i, br in enumerate(brojevi_profila):
        ime = f'KON{br}A40.DAT'
        putanja = os.path.join(folder, ime)
        if os.path.exists(putanja):
            df = konvertuj_dat_u_df(putanja)
            df['y'] = i * razmak_profila  # širina preseka (udaljenost profila)
            profili.append(df)
        else:
            print(f"[UPOZORENJE] Fajl ne postoji: {ime}")
    return profili

def napravi_horizontalni_presek(profili, dubina, korak_x=0.8):
    svi_podaci = pd.concat(profili, ignore_index=True)
    x_min = svi_podaci['x'].min()
    x_max = svi_podaci['x'].max()
    y_min = svi_podaci['y'].min()
    y_max = svi_podaci['y'].max()

    x_koord = np.arange(x_min, x_max + korak_x, korak_x)
    y_koord = np.arange(y_min, y_max + razmak_profila, razmak_profila)

    otpornost = np.full((len(y_koord), len(x_koord)), np.nan)

    for i, yv in enumerate(y_koord):
        df_prof = svi_podaci[svi_podaci['y'] == yv]
        if df_prof.empty:
            continue
        for j, xv in enumerate(x_koord):
            df_x = df_prof[np.isclose(df_prof['x'], xv, atol=1e-3)]
            if df_x.empty:
                df_prof_x = df_prof[(df_prof['x'] >= xv - korak_x/2) & (df_prof['x'] <= xv + korak_x/2)]
                if df_prof_x.empty:
                    continue
                najblizi_x = df_prof_x.iloc[(df_prof_x['x'] - xv).abs().argmin()]['x']
                podaci_za_interp = df_prof[df_prof['x'] == najblizi_x]
            else:
                podaci_za_interp = df_x

            z_vals = podaci_za_interp['z'].values
            rho_vals = podaci_za_interp['ρa'].values

            if len(z_vals) < 2 or dubina < z_vals.min() or dubina > z_vals.max():
                continue

            f_interp = interp1d(z_vals, rho_vals, bounds_error=False, fill_value=np.nan)
            otpornost[i, j] = f_interp(dubina)

    return x_koord, y_koord, otpornost

def snimi_dat_format(x_koord, y_koord, otpornost, dubina, ime_fajla):
    with open(ime_fajla, 'w', encoding='utf-8') as f:
        f.write(f"hor.pres.prof {','.join(str(b) for b in brojevi_profila)} nivo 1\n")
        f.write(f"{1.00:.2f}\n")
        f.write(f"{1}\n")
        f.write(f"{int(dubina*100)}\n")
        f.write(f"{1}\n")
        f.write(f"0\n")

        for i, yv in enumerate(y_koord):
            for j, xv in enumerate(x_koord):
                if np.isnan(otpornost[i,j]):
                    continue
                f.write(f"{xv:.3f},{yv:.2f},{otpornost[i,j]:.2f}\n")

        for _ in range(6):
            f.write("0\n")

    print(f"[INFO] Snimljen DAT fajl: {ime_fajla}")

if __name__ == '__main__':
    profili = ucitaj_sve_profili()
    for dubina in dubine_za_prekseke:
        x_koord, y_koord, otpornost = napravi_horizontalni_presek(profili, dubina)
        ime_dat = f'horizontalni_presek_{int(dubina*100)}cm.dat'
        snimi_dat_format(x_koord, y_koord, otpornost, dubina, ime_dat)

Evo kratkog uputstva za korišćenje programa hor.py, koji vizualizuje horizontalne preseke geoelektrične otpornosti na osnovu .DAT fajlova generisanih programom vert2horDAT.py.

Naziv programa: `hor.py`

Namena:

Prikaz i automatsko čuvanje kontur slika horizontalnih preseka prividne geoelektrične otpornosti (ρa) na zadatim dubinama. Program koristi .DAT fajlove koje je prethodno napravio vert2horDAT.py.

Ulazni fajlovi:

Tekstualni .DAT fajlovi u formatu: x, y, ρa počevši od šestog reda (prvih 5 se ignorišu kao zaglavlje).
Primer imena fajlova: horizontalni_presek_50cm.dat horizontalni_presek_75cm.dat ...

Kako koristiti program:

Postavi hor.py u isti folder sa .DAT fajlovima koje je napravio vert2horDAT.py.
U terminalu pokreni: bashCopyEditpython3 hor.py
Program će:
- automatski pronaći sve .dat fajlove
- interpolirati podatke u raster
- nacrtati i sačuvati kontur slike u .jpg formatu
  (npr. horizontalni_presek_50cm_konture.jpg, horizontalni_presek_75cm_konture.jpg)

Rezultat:

Svaka slika prikazuje horizontalni sloj prividne otpornosti u ravni x-y na fiksnoj dubini.
Boje predstavljaju vrednosti otpornosti (Ohm·m), sa izolinijama.
Osa X = dužina preseka, osa Y = širina (udaljenost između profila)

Zavisnosti koje treba imati:

Instaliraj ih ako već nisu instalirane:

pip3 install pandas numpy matplotlib scipy

Napomena:

Naslov slike automatski prikazuje ime fajla
Dubine nisu prikazane direktno na slici — ali su sadržane u imenu fajla (npr. 50cm = 0.5m)

Programski kod za hor.py

# hor.py
# The MIT License (MIT)
# Copyright (c) 2025 Aleksandar Maričić
#
# Ovim se omogućava bilo kome da koristi, kopira, menja, spaja, objavljuje,
# distribuira, daje podlicencu i/ili prodaje kopije ovog softverskog programa,
# uz uslov da u svim kopijama ili značajnim delovima softverskog programa bude
# uključena sledeća obavest:
#
# Copyright (c) 2025 Aleksandar Maričić
#
# Ovaj softverski program je pružen "takav kakav jeste", bez bilo kakvih garancija,
# izričitih ili impliciranih, uključujući, ali ne ograničavajući se na, garancije o
# prikladnosti za prodaju ili pogodnosti za određenu svrhu. U svakom slučaju, autori
# ili nosioci prava nisu odgovorni za bilo kakvu štetu ili druge obaveze koje mogu nastati
# usled upotrebe ovog softverskog programa.
# Prikaz i automatsko čuvanje kontur slika horizontalnih preseka prividne geoelektrične otpornosti (ρa)
# na zadatim dubinama. Program koristi .DAT fajlove koje je prethodno napravio vert2horDAT.py.



import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import griddata
import os
import sys
import glob

def konvertuj_dat_u_df(ulazni_fajl):
    with open(ulazni_fajl, 'r', encoding='utf-8') as f:
        linije = f.readlines()

    podatak_pocetak = 5
    podaci = []
    for linija in linije[podatak_pocetak:]:
        linija = linija.strip()
        if linija == '':
            continue
        try:
            x, z, rho = map(lambda s: float(s.replace(',', '.')), linija.split(','))
            podaci.append([x, z, rho])
        except ValueError:
            print(f"[UPOZORENJE] Preskačem liniju: {linija}")
            continue

    df = pd.DataFrame(podaci, columns=['x', 'z', 'ρa (Ohm·m)'])
    print(f"[INFO] Učitano {len(df)} redova iz '{ulazni_fajl}'")
    return df

def load_data(filename):
    ext = os.path.splitext(filename)[1].lower()
    if ext == '.dat':
        return konvertuj_dat_u_df(filename)
    elif ext == '.csv':
        print(f"[INFO] Učitavanje CSV fajla '{filename}'...")
        return pd.read_csv(filename)
    else:
        raise ValueError(f"Nepoznat format fajla: {ext}")

def interpolate_to_grid(df, res=0.1, method='cubic'):
    x, z, rho = df['x'], df['z'], df['ρa (Ohm·m)']
    xi = np.arange(min(x), max(x), res)
    zi = np.arange(min(z), max(z), res)
    Xi, Zi = np.meshgrid(xi, zi)
    Ri = griddata((x, z), rho, (Xi, Zi), method=method)
    return Xi, Zi, Ri

def save_contour_as_image(Xi, Zi, Ri, output_filename='ert_konture.jpg'):
    plt.figure(figsize=(10, 5))
    levels = np.linspace(np.nanmin(Ri), np.nanmax(Ri), 15)
    contour = plt.contourf(Xi, Zi, Ri, levels=levels, cmap='plasma')
    plt.contour(Xi, Zi, Ri, levels=levels, colors='black', linewidths=0.5)
    plt.colorbar(contour, label='ρa (Ohm·m)')
    plt.gca().invert_yaxis()
    plt.gca().set_aspect('equal')
    plt.xlabel('Dužina (m)')
    plt.ylabel('Širina (m)')
    plt.title(f" '{output_filename}' i prividne otpornosti (ERT konture)")
    plt.grid(True)
    plt.tight_layout()
    plt.savefig(output_filename, dpi=300)
    plt.close()
    print(f"[INFO] Sačuvano kao '{output_filename}'")

def main():
    folder = '.'  # trenutni folder
    dat_files = glob.glob(os.path.join(folder, '*.dat')) + glob.glob(os.path.join(folder, '*.DAT'))
    if not dat_files:
        print("[GREŠKA] Nema .DAT fajlova u folderu.")
        sys.exit(1)

    for dat_file in sorted(dat_files):
        print(f"\n[INFO] Procesiram fajl: {dat_file}")
        df = load_data(dat_file)
        Xi, Zi, Ri = interpolate_to_grid(df, res=0.1)
        output_image = os.path.splitext(dat_file)[0].lower() + '_konture.jpg'
        save_contour_as_image(Xi, Zi, Ri, output_image)

if __name__ == '__main__':
    main()

Vertikalni preseci