<jupyter_start><jupyter_text>Datum transformation: Aratu -> SAD69 The first challenge to create a project with both seismic and well data is to stabilish a common datum. The seismic data in Brazil is stored in SAD69 while well data, mainly the old ones, are stored in the Aratu datum. Here I've used nine wells to find the difference between the two datums. For a more precise location it's used the Molodensky transformation. But the error between the approaches is less than 1 meter. **Important:** The codes below are valid for Aratu datum, central meridian 39. The Aratu datum is limited to some regions of Brazil. See the map below (Anciães, 2003):<jupyter_code>from IPython.display import Image
Image("imagens/aratu-brazil.png")<jupyter_output><empty_output><jupyter_text>The Aratu datum is identified by the yellow color.  The Molodensky simplified formula $\Delta \phi^\circ= \frac{1}{M_1} \{(a_1 \Delta f + f_1 \Delta a)  \ sin \ 2 \phi_1 - \Delta X \ sin \ \phi_1 \ \ cos \ \lambda_1 - \Delta Y \ sin \ \phi_1 \ \ sin \ \lambda_1 + \Delta Z \ \ cos \ \phi_1\} \frac{180}{\pi}$ $\quad \quad M = \frac{a(1 - e^2)}{(1 - e^2 \ sin^2 \ \phi)^{3/2}}$ $ \Delta \lambda^\circ = \frac{1}{N_1 \ cos \ \phi_1} \{-\Delta X \ sen \ \lambda_1 + \Delta Y cos \ \lambda_1\} \frac{180}{\pi}$ $\quad \quad N = \frac{a}{(1 - e^2 \ sin^2 \phi)^{1/2}}$ $\quad \quad \quad e^2 = f(2-f) \quad \quad e'^2 = \frac{e^2}{1-e^2}$ Where, $M$ and $N$ = meridional and prime vertical radii of curvature, respectively, at a given latitude; $\phi$ = geodetic latitude; $\lambda$ = geodetic longitude; $a$ = semi-major axis of the local geodetic datum ellipsoid; $b$ = semi-minor axis of the local geodetic datum ellipsoid; $f$ = flattening of the local geodetic datum ellipsoid; $e$, $e'$ = first and second excentricities of the ellipsoid; $\Delta a, \Delta f$ = differences between the semi-major axis and flattening of the local geodetic datum ellipsoid.  Delta values for the WGS-84 datum in reference to SAD-69 (IBGE - resolução nº 23 de 21 de fevereiro de 1989): $\Delta X = +66.87 \ m \ \pm \ 0.43 \ m$   $\Delta Y = -4.37 \ m \ \pm \ 0.44 \ m$   $\Delta Z = + 38.52 \ m \ \pm \ 0.40 \ m$ All $\Delta$ quantities are formed by subtracting local geodetic datum ellipsoid value from the reference Ellipsoid values. That values are summed (or subtracted) to obtain the transformed values as shown below: $ \phi_2^\circ = \phi_1^\circ + \Delta \phi^\circ \quad \quad $ or $\quad \quad  \Delta \phi^\circ = \phi_2^\circ - \phi_1^\circ $ $ \lambda_2^\circ = \lambda_1^\circ + \Delta \lambda^\circ \quad \quad $ or $ \quad \quad \Delta \lambda^\circ = \lambda_2^\circ - \lambda_1^\circ$ Where: $\phi_1^\circ$ = local geodetic latitude; $\phi_2^\circ$ = reference datum geodetic latitude; $\lambda_1^\circ$ = local geodetic longitude; $\lambda_2^\circ$ = reference datum geodetic longitude. **Parameters of the SAD-69 datum:**     a = 6378160 m     b = 6356774.719 m   $f = 1/298.25$    ** Parameters Aratu datum:**    a = 6378388 m  $f = 1/297$   Transformation from SAD-69 to WGS-84 (IBGE - resolução nº 23 de 21 de fevereiro de 1989): $X_2 = X_1 + 66.87$   $Y_2 = Y_1 - 4.37$   $Z_2 = Z_1 + 38.52$ Applying the logic above for nine wells with coordinates in SAD69 and in ARATU, extracted form the website [epsg.io](https://epsg.io/transforms_srs=20824&t_srs=29194&x=304413.0000000&y=7549548.0000000 "epsg.io"), I can find the $\Delta$ value for the transformation between Aratu and SAD-69 datums:<jupyter_code>import pandas as pd
d = {'X_SAD69': [304576.07, 294313.11, 311693.03, 321554.01, 294887.80, 327191.96, 323310.99, 303580.06], 'Y_SAD69': [7549440.14, 7540343.16, 7573468.07, 7564824.10, 7553656.52, 7614185.95, 7587296.03, 7565617.09], 'Y_ARATU': [7549548,7540451,7573576,7564932,7553764.4,7614294,7587404,7565725], 'X_ARATU': [304413,294150,311530,321391,294724.7,327029,323148,303417]}
df = pd.DataFrame(data=d)
df<jupyter_output><empty_output><jupyter_text>Determination of $\Delta$'s values by the mean value:<jupyter_code>import numpy as np<jupyter_output><empty_output><jupyter_text>$\Delta X = X_2 - X_1$<jupyter_code>dif_X =  df['X_SAD69'] - df['X_ARATU']
dif_X
max_X = np.amax(dif_X)
max_X
min_X = np.amin(dif_X)
min_X
Delta_X = (max_X + min_X) / 2
Delta_X
error_X = abs(Delta_X) - abs(max_X)
error_X<jupyter_output><empty_output><jupyter_text>$ \Delta Y = Y_2 - Y_1$<jupyter_code>dif_Y = df['Y_SAD69'] - df['Y_ARATU']
dif_Y
min_Y = np.amin(dif_Y)
min_Y
max_Y = np.amax(dif_Y)
max_Y
Delta_Y = (max_Y + min_Y) / 2
Delta_Y
error_Y = abs(Delta_Y) - abs(min_Y)
error_Y<jupyter_output><empty_output><jupyter_text>Parameters for transformation from Aratu (M 39 - Campos Basin) to SAD-69: $\Delta X = 163.035 m \ \pm \ 0.074 m$  $\Delta Y = -107.944 m \ \pm \ 0.104 m $  Reading the coordinates of .las files<jupyter_code># import here the .las files already integrated by the previous codes.
import lasio
las_files = ('perfis/1RJS_0051__RJ_1RJS_0051__RJ_BHC_00003.las', 'perfis/1RJS_0070__RJ_1RJS_0070__RJ_BHC_00008.las', 'perfis/1RJS_0087__RJ_1RJS_0087__RJ_BHC_00011.las', 'perfis/1RJS_0191__RJ_1RJS_0191__RJ_FDC_00013.las','perfis/1RJS_0203__RJ_1RJS_0203__RJ_FDC_00015.las', 'perfis/1RJS_0225__RJ_1RJS_0225__RJ_FDC_CNL_00018.las', 'perfis/1RJS_0351__RJ_1RJS_0351__RJ_DIT_BHC_GR_00023.las', 'perfis/1RJS_0426__RJ_1RJS_0426__DT_digitized_0001.las' )
N_files = len(las_files)
N_files
for i in range(N_files):
    globals()['log%s' % i] = lasio.read(las_files[i])
# for reading only one .las file use the following code:
# las_files = lasio.read('perfis/1RJS_0051__RJ_1RJS_0051__RJ_BHC_00003.las')
# las_files.well.fl1.value<jupyter_output><empty_output><jupyter_text>Reading X coordinates<jupyter_code>float(log0.well.fl1.value[5:15])   # identifying the location inside of the line
X_Aratu = []
for i in range(N_files):
    X = float(globals()['log%s' % i].well.fl1.value[5:15])
    X_Aratu = np.append(X_Aratu, X)
X_Aratu = np.asarray([X_Aratu])
X_Aratu.T<jupyter_output><empty_output><jupyter_text>Reading Y coordinates<jupyter_code>float(log1.well.fl1.value[25:34])    # identifying the location inside of the line
Y_Aratu = []
for i in range(N_files):
    Y = float(globals()['log%s' % i].well.fl1.value[25:34])
    Y_Aratu = np.append(Y_Aratu, Y)
Y_Aratu = np.asarray([Y_Aratu])
Y_Aratu.T
Aratu = np.concatenate((X_Aratu.T, Y_Aratu.T), axis=1)
Aratu<jupyter_output><empty_output><jupyter_text>Transforming coordinates from Aratu (M 39) to SAD69: $X_2 = X_1 + \Delta X $   $Y_2 = Y_1 + \Delta Y$<jupyter_code>X_2 = (np.array([Aratu[:, 1] + Delta_X]))
X_2.T
Y_2 = (np.array([Aratu[:, 0] + Delta_Y]))
Y_2.T<jupyter_output><empty_output><jupyter_text>Verifying the similarity between the results:<jupyter_code>np.allclose(df['X_SAD69'], X_2, atol=1e-4)
np.allclose(df['Y_SAD69'], Y_2, atol=1e-4)<jupyter_output><empty_output><jupyter_text>Exporting the wells track     col 1 - well name      col 2 - X       col 3 - Y       col 4 - KB       col 5 - Total Depth       col 6 - Ground Level elevation (GL)       col 7 - Unique Well ID  Reading additional data for creating the well track:<jupyter_code>well_name = []
KB = []
total_depth = []
ground_elevation = []
well_api = []


for i in range(N_files):
    well = str(globals()['log%s' % i].well.well.value[5:9])
    well_name = np.array([np.append(well_name, well)])
    
    K = str(globals()['log%s' % i].well.EKB.value)
    KB = np.array([np.append(KB, K)])
    
    total = str(globals()['log%s' % i].well.STOP.value)
    total_depth = np.array([np.append(total_depth, total)])
    
    ground = float(globals()['log%s' % i].well.EGL.value)
    ground_elevation = np.array([np.append(ground_elevation, ground)])
    
    well = str(globals()['log%s' % i].well.api.value)
    well_api = np.array([np.append(well_api, well)])
    
well_api.T
SAD69 = np.hstack((well_name.T, X_2.T, Y_2.T, KB.T, total_depth.T, ground_elevation.T, well_api.T))
SAD69 
## add well manually
add_well = ['1CST_0001_RJ', 292702.935, 7564517.955, 8, 2620.9, 4, 10000001]
SAD69 = np.insert(SAD69, [1], add_well, axis=0)
SAD69
np.savetxt('well_track.txt', SAD69, fmt='%s', delimiter='    ', newline='\n', header='', footer='', comments='# ')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14604957870313837, 'fraction_numerical': 0.0885333984613506, 'mean_word_length': 4.33203125, 'pattern_counts': {'":': 0, '<': 45, '<?xml version=': 0, '>': 46, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 2, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '646737', 'license': 'MIT', 'n_tokens_mistral': 3752, 'n_tokens_neox': 3302, 'n_words': 814, 'path': 'Transformacao-aratu-sad69.ipynb', 'repo_name': 'raquelsilva/Seismic-well-tie', 'size': '8189'}