Métodos potenciales de prospección, FCAG, 2024.
Leemos el archivo cca-salado.dat, y procesamos para quedarnos con las columnas de interés. La estructura del dato podemos conocerla mirando el encabezado del archivo por medio de comandos GNU/Linux:
$> head -15 cca-salado.dat
#RELEVAMIENTO GRAVIMETRICO EXTRAIDO DEL:
#CATALOGO GENERAL DE ESTACIONES GRAVIMETRICAS
#OBSERVATORIO ASTRONOMICO UNLP
#SERIE GEODESICA - TOMO IX.
#
#MATEO J. LEVIN E et al; 1976
#
#
# LONGITUD LATITUD COTA GOBSERV GBOUGUER
# (ø) (ø) (m) (mGal) (mGal)
-58.56335 -34.63831 23.18 979705.0 -5.600
-59.42500 -34.64497 40.69 707.6 -0.100
-58.61838 -34.64832 26.20 704.8 -6.000
-59.36666 -34.67001 40.71 708.6 -1.200
Observamos que al final del archivo también hay comentarios:
$> tail -8 cca-salado.dat
-57.81001 -38.09832 55.29 036.9 34.900
-57.83996 -38.11501 58.64 033.3 30.400
-57.89501 -38.11997 55.24 030.0 26.100
-57.98334 -38.13837 56.60 032.3 27.000
#COMENTARIO: GOBSERV REPRESENTA LA GRAVEDAD OBSERVADA
# YA CORREGIDA POR LATITUD
#REFERENCIA: 1930
Leemos las coordenadas de los puntos de observación y la anomalía de Bouguer del archivo.
import numpy as np # para calcular, NUMerical PYthon.
import matplotlib.pyplot as plt # para graficar, matplotlib.
d = np.loadtxt('cca-salado.dat') # las lineas que comienzan con '#' son ignoradas
x, y, gz = d[:,0],d[:,1],d[:,4] # columnas del dato en cada vector
Describimos el dato calculando la cantidad de estaciones y los rangos de longitud y latitud abarcados en la campaña:
Número de estaciones: 161 Longitud [grados] : -60.99998 , -56.67497 Latitud [grados] : -38.13837 , -34.63831
Visualizamos la ubicación de las observaciones:
# Gráfico (simplificado)
plt.figure()
plt.scatter(x, y, c = gz)
plt.show()
Para interpretar los datos, generamos una cuadrícula a los fines de graficar el campo observado, la superficie de tendencia y el campo residual.
# Construcción de la cuadrícula
nx,ny = 50, 50
grid_x = np.linspace(min(x), max(x), nx)
grid_y = np.linspace(min(y), max(y), ny)
X,Y = np.meshgrid(grid_x, grid_y)
# Gráfico (simplificado)
plt.figure()
plt.scatter(X, Y)
plt.scatter(x, y)
plt.show()
La cuadrícula rectangular es 50 x 50
Interpolamos el dato observado en la cuadrícula para graficar. Para ello, utilizamos la subrutina SciPy scipy.interpolate.griddata:
from scipy.interpolate import griddata
GZ = griddata(points=(x,y), values=gz, xi=(X,Y), method='cubic', fill_value=np.nan)
# Gráfico (simplificado)
plt.figure()
plt.contour(X, Y, GZ, 10, colors="black")
plt.imshow(GZ , origin="lower", interpolation=None,extent=[np.min(x),np.max(x),np.min(y),np.max(y)])
plt.scatter(x, y)
plt.show()
Como podemos observar, la interpolación de los valores de $g_z$ permite facilitar la interpretación de las observaciones. En estos gráficos, los puntos que caen fuera del polígono que encierra las observaciones originales (convex hull) no son gráficados, ya que optamos por la opción fill_value=np.nan.
Ajustamos al dato $g_z$ polinomios de grado $M$,
\begin{equation} \hat{g}_z(x,y;M) = \begin{cases} c, & \; M=0;\\ ax + by + c, & \; M=1;\\ ax^2 + by^2 + cxy + dx + ey + f, & \; M=2. \end{cases} \end{equation}Analizamos para $M=0,1,2$.
En notación matricial, el vector de coeficientes es $\mathbf{m} = (a,b,c,...)^T$ con dimensiones $m \times 1$. El problema directo está dado por el sistema lineal
$$\mathbf{A}\mathbf{m} = \mathbf{g},$$donde $\mathbf{A}$ $n\times m$ es la matriz de diseño y $\mathbf{g}$ los valores de la superficie de tendencia agrupados en un vector columna $n \times 1$.
El número de filas de $\mathbf{A}$ viene dado por el número de observaciones $n$ y el número de columnas $m$ por el número de coeficientes a ajustar. Por ejemplo:
Para $M=0$, el vector de coeficientes es $\mathbf{m}=(c)$ y la matriz de diseño $\mathbf{A}$ tiene dimensiones $n \times 1$: $\mathbf{A}=\big[ \mathbf{1} \big]$,
Para $M=1$, $\mathbf{m}=(a,b,c)^T$ y $\mathbf{A}$ es $n \times 3$: $\mathbf{A}=\big[\mathbf{x}| \mathbf{y} | \mathbf{1}\big]$;
donde $\mathbf{x}$ es el vector columna $n \times 1$ de coordenadas en x, $\mathbf{y}$ es el vector columna $n \times 1$ de coordenadas en y, mientras que $\mathbf{1}$ es un vector $n \times 1$ con entradas iguales a $1$.
Calculamos la matriz de diseño $\mathbf{A}$ para un cierto valor del grado del polinomio $M$:
# Caso M = 1: m = (a,b,c)^T es 3x1
m = np.zeros(shape=(3,1))
A = np.zeros(shape=(len(gz),len(m)))
A[:,0] = x
A[:,1] = y
A[:,2] = np.ones(len(gz))
Caso M = 1 : La matriz A tiene dimensiones (161, 3)
Realizamos el ajuste de los coeficientes que definen la superficie de tendencia por medio del problema inverso, obteniendo $\mathbf{m}^*$. El vector de parámetros $\mathbf{m}^*$ es aquel que produce la mejor aproximación (en el sentido de los cuadrados mínimos) a los datos observados según el modelo dado en $\mathbf{A}$:
$$\mathbf{A}\mathbf{m}^* \approx \mathbf{g^o},$$siendo $\mathbf{g^o}$ los datos observados agrupados en un vector columna $n \times 1$.
Obtenemos $\mathbf{m}^*$ utilizando la subrutina SciPy scipy.linalg.lstsq:
from scipy.linalg import lstsq
m,_,_,_ = lstsq(A,gz)
m^* = [2.327000 -4.305000 -2.794000]
Una vez realizado el ajuste sobre lo datos originales, procedemos a interpolar para obtener la superficie de tendencia sobre la cuadrícula generada en un apartado anterior. La superficie de tendencia calculada sobre la cuadrícula facilitará la visualización y la interpretación.
Para realizar la interpolación, calculamos la matriz de diseño A_grid en los puntos de la cuadrícula:
nx = len(grid_x)
ny = len(grid_y)
A_grid = np.zeros(shape=(nx*ny,len(m)))
# Caso M = 1
A_grid[:,0] = X.reshape(nx*ny)
A_grid[:,1] = Y.reshape(nx*ny)
A_grid[:,2] = np.ones(nx*ny)
Caso M = 1 : La matriz [A_grid] para interpolar tiene dimensiones (2500, 3)
Luego, por medio del problema directo, obtenemos la superficie de tendencia, $\mathbf{\hat{g}}=\mathbf{A}\mathbf{m}^*$, y el campo residual, $\mathbf{\Delta g}=\mathbf{\hat{g}^o}-\mathbf{\hat{g}}$, en los puntos de la cuadrícula. Aquí $\mathbf{\hat{g}^o}$ es el campo observado que ha sido interpolado sobre la cuadrícula.
GZ_fit = A_grid @ m # En Python>3, '@' es el producto entre una matriz y un vector
GZ_res = GZ - GZ_fit.reshape(nx,ny) # campo residual interpolado
Notemos que no es necesario interpolar para realizar el ajuste de tendencias, ya que estamos utilizando cuadrados mínimos. Interpolar antes de realizar la inversión representaria: a) hacer más cálculos (habría más muestras, lo que conduce a una matriz de diseño de mayor dimensión) y b) ajustar sobre dato interpolado y no sobre el dato observado.
Visualizamos el dato original $g_z$ interpolado, el campo modelado según el ajuste de la superficie de tendencia $\hat{g}_z$ y el campo residual $\Delta g_z$.
# Gráficos (simplificado):
for CUADRICULA in [GZ, GZ_fit.reshape(nx,ny), GZ_res]:
plt.figure()
plt.contour(X, Y, CUADRICULA, colors="black")
plt.imshow(CUADRICULA, origin="lower", interpolation=None, extent=[np.min(x),np.max(x),np.min(y),np.max(y)])
plt.show()
Observamos que las líneas de contorno son suaves debido a que se han interpolado los resultados con el objeto de graficar e interpretar. En estos gráficos, los puntos que caen fuera de la región de interés (convex hull) no son graficados.
En otra práctica, trabajaremos con el método de la capa equivalente para realizar la interpolación.
Eso es todo por hoy.
Repasamos lo analizado en la cátedra de análisis de señales y estadística. Luego aplicamos el método geoestadístico de Kriging ordinario para interpolar el dato de anomalía observado. En este ejemplo, utilizamos un variograma esférico y la libería pykrige.
from pykrige.ok import OrdinaryKriging
variogram_model = "spherical"
coordinates_type = "geographic"
OK = OrdinaryKriging(x, y, gz, variogram_model=variogram_model, coordinates_type=coordinates_type)
gz_int,_ = OK.execute('grid', grid_x, grid_y)
gz_int = gz_int.data
...
Podemos interpolar el archivo original utilizando el software GMT (Generic Mapping Tools) escribiendo unas pocas líneas en la terminal. El dato original va de longitudes $-61^\circ$ a $-56.6^\circ$ y latitudes $-38^\circ$ a $-34.5^\circ$. El método de interpolación utilizado es el continuous curvature splines in tension. Generamos el dato en una cuadrícula e interpolado en binario (formato NetCDF) y luego lo convertimos a formato ASCII para utilizarlo en esta clase. Utilizamos el parámetro $T=0.25$ recomendado para datos de métodos potenciales.
Nota: Si no tenemos GMT instalado, las siguientes lineas no correrán.
Filtramos el dato con herramientas superpoderosas de GNU/Linux:
$> grep -v '#' cca-salado.dat > cca-salado_clean.dat # quito comentarios
$> awk '{print $1,"",$2,"",$5}' "cca-salado_clean.dat" > AB-cca-salado.dat
Y utilizamos el maravilloso GMT para interpolar:
$> gmt surface AB-cca-salado.dat -R-61/-56.6/-38.2/-34.6 -I5m -T0.25 -Gtmp.grd # interpolamos
Cambiamos el binario (en formato NetCDF) a formato ASCII para leer fácilmente en la notebook y graficar:
$> gmt grd2xyz tmp.grd > AB-cca-salado-interpolado.dat # cambiamos el formato a ASCII
