Lectura del dato¶

Utilizamos el archivo datos2D.dat de anomalías de Bouguer dispuestos en una cuadrícula. Para comenzar, podemos emplear algunos comandos GNU/Linux para observar la estructura del dato en la terminal.

$> head "datos2D.dat"

# x[km]      y[km]    AB[mGal]

Viendo el encabezado, notamos que el dato está presentado en la forma

$x$ (km)	$y$ (km)	$g_z$ (mGal)
...	...	...
...	...	...
...	...	...

Ahora procedemos a leer el archivo para procesar. Para ello podemos escribir:

data     = np.loadtxt('datos2D.dat')       # las líneas que comienzan con '#' en el archivo son ignoradas
x, y, gz = data[:,0], data[:,1], data[:,2] # columnas del dato en cada vector

N  = len(gz)           # número de observaciones de AAB [mgal]
nx = len(np.unique(x)) # número de coordenadas únicas en x
ny = len(np.unique(y)) # número de coordenadas únicas en y

print(f"Contamos con {N} observaciones.")
print(f"nx={nx}, ny={ny}, (nx * ny = {nx*ny}).")

GZ = gz.reshape(ny,nx) # genero matriz con el vector de mediciones de gravedad (sólo para graficar)

Contamos con 176 observaciones.
nx=16, ny=11, (nx * ny = 176).

Ajuste de polinomios de tendencia¶

Ajustamos al dato $g_z$ polinomios de grado $M$,

$$ \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} $$

$\newcommand{\B}[1]{\mathbf{#1}}$ Consideramos $M=1$ y calculamos la matriz de diseño $\B{A}= [\B{x}\: \B{y}\: \B{1}]$:

x    = x.reshape(N,1)
y    = y.reshape(N,1)
ones = np.ones((N,1))
A    = np.hstack((x, y, ones))  # A = [ x y 1 ]

print(f"Tamaño de A: {np.shape(A)}")

Tamaño de A: (176, 3)

Resolvemos el sistema $\mathbf{A}\mathbf{m}=\mathbf{g^o}$, para $\mathbf{m}$ por cuadrados mínimos. Aquí $\mathbf{g^o}$ es el dato observado dispuesto en un vector columna. La solución dada por cuadrados mínimos la denotamos $\mathbf{m}^*$ y en código m_estimated:

from scipy.linalg import lstsq
m_estimated,_,_,_ = lstsq(A,gz)

Grado del polinomio de ajuste, M = 1
Vector de coeficientes, m*       = [-0.022066 -0.209915 8.734554]

Como valor agregado, podemos verificar el resultado si utilizamos un algoritmo de gradientes conjugados (CG) programado en una clase preliminar:

def CG(A,b):
    ...
    return m

m_estimated = CG(A.T @ A,A.T @ gz)

CG: Según la tolerancia 1e-10, CG converge en la iteración 4
vector de coeficientes, m* =  [-0.022066 -0.209915 8.734554]

Calculo de la superficie de tendencia y el mapa residual¶

Con el vector de coeficientes $\mathbf{m^*}$, obtenemos por medio del problema directo la superficie de tendencia en los puntos observados:

$$\mathbf{\hat{g}} = \mathbf{A}\mathbf{m^*},$$

y el campo residual,

$$\Delta\mathbf{g} = \mathbf{g^o}-\mathbf{\hat{g}}.$$

En python podemos generar la superficie de tendencia y el campo residual así:

# Superficie de tendencia:
GZ_fit = (A @ m) # en python 3, '@' es producto matricial
GZ_fit = GZ_fit.reshape(ny,nx) # construyo cuadrícula del ajuste 
# Campo residual (resta de cuadrículas):
diff = GZ - GZ_fit

Resultados¶

Visualizamos el dato original $g_z$, la superficie ajustada según mínimos cuadrados $\hat{g}_z$ y el residuo $\Delta g_z=g_z-\hat{g}_z$. En los gráficos, utilizamos la misma escala vertical para el dato original y la superficie de tendencia. Para visualizar el campo residual, utilizamos una escala diferente centrada en 0 mGal. Observamos que las líneas de contorno no son suaves debido a que los valores no han sido interpolados en una cuadrícula más densa.

Eso es todo por hoy.

Extra¶

Lectura del dato cuadriculado e interpolado¶

Leemos el archivo datos2D_interpolado.dat de anomalías de Bouguer dispuestas en una cuadrícula e interpoladas con un método específico para los métodos potenciales. Procedemos de la misma manera, realizando el ajuste para $M=1$ y graficando los resultados.

Interpolando el archivo original utilizando GMT¶

Podemos interpolar el archivo original utilizando GMT (Generic Mapping tools) y escribiendo unas pocas líneas en la terminal. El dato original tiene $n_x=16$ y $n_y=11$. El método de interpolación utilizado es el continuous curvature splines in tension. Generamos el dato cuadrículado 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.

$> which GMT # si esta instalado, arroja el path al ejecutable