estoy buscando plotear una sísmica (en gris) con matplotlib y encima de la sísmica una correlación (que seria de varios colores), la correlación hace alusión a los horizontes rastreados previamente, horizon_number es el número de horizontes que se quiere ver (también es opción visualizar todos con "all"), view es para diferenciar entre inline y xline (tipos de lineas en la sísmica), line_number es el número de línea que se va a mostrar y selem es un parámetro para la dilation, el código es sencillo pero la gráfica sale incompleta y no sé por qué. También estoy abierta a opciones si a alguien se le ocurre algo con xarray. Acepto sugerencias y recomendaciones por favor. Adjunto código e imagen.
def ShowHorizon (seismic_data, horizons, horizon_number, view, line_number, selem):
correlated_cube= np.zeros(seismic_data['data'].shape, dtype=np.float32)
if horizon_number == 'all':
num=0
correlated_cube= np.zeros(seismic_data['data'].shape, dtype=np.float32)
for i in range (0, len(horizons)):
reflector=horizons[i]
num=num+1
for p in reflector:
x=p[0]
y=p[1]
z=p[2]
correlated_cube[x,y,z]=num
else:
reflector=horizons[horizon_number]
for p in reflector:
x=p[0]
y=p[1]
z=p[2]
correlated_cube[x,y,z]=correlated_cube[x,y,z]+1 #+1 FOR Å FÅ EN IDE OM HITRATE
fig, ax = plt.subplots()
if view == 'inline':
A=correlated_cube[:,line_number,:]
I=seismic_data.isel(iline=line_number)
variables1 = I.data_vars
arrays1 = [ I[i].values.flatten() for i in variables1 ]
arrays1 = [i for j in arrays1 for i in j ]
array1 = np.array(arrays1).reshape(-1,1)
B= array1.copy()
if view == 'crossline':
A=correlated_cube[:,:,line_number]
X=seismic_data.isel(xline=line_number)
variables2 = X.data_vars
arrays2 = [ X[i].values.flatten() for i in variables2 ]
arrays2 = [i for j in arrays2 for i in j ]
array2 = np.array(arrays2).reshape(-1,1)
B= array2.copy()
A = dilation(A, selem)
masked_data = np.ma.masked_where(A < 1, A)
plt.imshow(B, 'gray', aspect='auto')
plt.imshow(masked_data, 'jet', alpha=0.9, aspect='auto')
if horizon_number == 'all':
cbar = plt.colorbar(ax=None)
cbar.set_label('Horizon number')
else:
cbar = plt.colorbar(ax=None)
cbar.set_label('Autotracking accuracy (Hit rate)', rotation=270)
cbar.ax.set_yticklabels([]) #no labels on colorbar
plt.xlabel('intersecting lines')
plt.title('Autotracked horizons superimposed onto seismic')
plt.ylabel('Samples')
plt.tight_layout()
plt.show()
return