Bathymetry of southern Rio Grande do Sul continental shelf

Anchoita Cruises

The data was available by the Hydroacoustic Laboratory (Furg, Stefa & Pinho)

import os
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import griddata
import rasterio
from rasterio.plot import show
import pickle

from shapely.geometry import Point
from shapely.geometry import Polygon

mapa_base = rasterio.open(r'd:\GUTO\1_Trabs\2_Plataforma\PlaRS\BaseGeografica\recorte_plataforma_rs3.tif')
show(mapa_base)

<AxesSubplot:>

'''
Load the ASCII files with latitude, longitude and depth
'''
files = os.listdir()

crz = []
for f in files:
    if f.endswith('txt') == True:
        print(f)
        with open(f) as io:
            d = io.read().splitlines()
            
        print(len(d))
        crz.append(d)

Anchoita_III_2010.txt
358952
Anchoita_II_2010.txt
342127
Anchoita_IV_2010.txt
320374
Anchoita_I_2010.txt
228680
Anchoita_V_2010.txt
315886

def isfloat(x):
    try:
        float(x)
        return True
    except ValueError:
        return False

'''
Check consistency and convert the strings to numbers and 
'''
crz2 = []

for c in crz:
    
    nli = []
    
    for li in c:
        
        liq = li.split()
            
        if len(liq) == 3:
            
            if isfloat(liq[0]) == True:
                 
                liq = [float(x) for x in liq]
                
                nli.append(liq)
     
    crz2.append(np.array(nli))
    
    print(len(nli))

358904
341710
320373
218530
315885

'''
Merge the datasets
''' 
junta = np.zeros((1,3))
for c in crz2[:3]:
    c = c[c[:,0] < -10]
    junta = np.vstack((junta, c))
    
junta = np.delete(junta, 0, axis=0)

'''
Check the tracks... some looks weird!
'''
for c in crz2:
    c = c[c[:,0] < -10]
    plt.plot(c[:,1], c[:,0])
    plt.show()

# a glance on the bathymetric data...
plt.figure(figsize=(20,4))
plt.plot(junta[:,2])

[<matplotlib.lines.Line2D at 0x1f4557f58b0>]

# elimitating some obvious spurious depth values
junta = junta[junta[:,2] > 3]

plt.figure(figsize=(20,4))
plt.plot(junta[:,2])

[<matplotlib.lines.Line2D at 0x1f455857670>]

plt.figure(figsize=(20,4))
plt.plot(np.diff(junta[:,2]))

[<matplotlib.lines.Line2D at 0x1f4558b44c0>]

'''
Find out the steep changes in the depth, which means interruption of data
'''
breaks = np.where(
                  np.abs(
                         np.diff(junta[:,2])
                  ) > 25
                  ) 

plt.figure(figsize=(20,4))
plt.plot(np.diff(junta[:,2]))
plt.plot(breaks, np.linspace(50, 50, len(breaks)), 'go')
plt.show()

'''
insert the start and end of the series for further processing
'''
breaks = np.append(breaks, len(junta))
breaks = np.insert(breaks, 0, 0)

'''
reduces the data by averaging every 20 values but avoiding the breaks!
'''
rdz = []
for i, b in enumerate(breaks[1:]):
    ints = np.arange(breaks[i-1], b, 20)
    
    for j in range(len(ints)-1):
        media = np.mean( junta[ints[j]:ints[j+1],:], axis=0)
        rdz.append(media)
    

junta = np.array(rdz)

'''
Finding by try & error the Albardão section
'''
plt.figure(figsize=(20,4))

p1 = 40200
p2 = 42000

plt.plot(junta[:,2])
plt.axvline(p1, color='red')
plt.axvline(p2, color='green')

plt.show()
plt.plot(junta[:,1], junta[:,0])
plt.plot(junta[p1,1], junta[p1,0], 'ro')
plt.plot(junta[p2,1], junta[p2,0], 'go')

[<matplotlib.lines.Line2D at 0x1f45d6be040>]

'''
Checking the seciont and further refining the selection (try&error)
'''
secao1 = junta[p1+205:p2+130, :]

plt.plot(secao1[:,2])
plt.show()

plt.plot(secao1[:,1], secao1[:,0], '.-')

[<matplotlib.lines.Line2D at 0x1f45d8ebf40>]

'''
Interpolating
'''
xi = np.arange(np.min(junta[:,1]), np.max(junta[:,1]), 1000/111120)
yi = np.arange(np.min(junta[:,0]), np.max(junta[:,0]), 1000/111120)

xx, yy = np.meshgrid(xi, yi)

points = np.vstack((junta[:,1], junta[:,0])).T

bat = griddata(points, junta[:,2], (xx,yy), method='linear')

bat.shape

(399, 329)

conts = [10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 120, 140]
contour_plot = plt.contourf(xi, yi, bat, conts, cmap='rainbow')
plt.plot(junta[:,1], junta[:,0], '.k', markersize=.01)
plt.axis('equal')
plt.colorbar()

<matplotlib.colorbar.Colorbar at 0x1f45f1375e0>

# Albardão lighthouse coordinates... (from Google Earth)

longitude = -52.706006
latitude = -33.202726

'''
load the mask digitized... see cells at the end of the notebook!
'''
with open('bathymetry_mask.pkl', 'rb') as io:
    mask = pickle.load(io)
    
# to close the polygon!
mask = np.vstack((mask, mask[0,:]))

fig, ax = plt.subplots(figsize=(13,11))

conts = [10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 120, 140]
contour_plot = ax.contourf(xi, yi, bat, conts, cmap='rainbow')
ax.plot(junta[:,1], junta[:,0], '.k', markersize=.01)
ax.axis('equal')
plt.colorbar(contour_plot)

ax.plot(longitude, latitude, 'ro')
ax.set_ylim(-35.5, -31)
ax.set_xlim(-54, -49)
show(mapa_base, ax=ax)

ax.plot(mask[:,0], mask[:,1], 'r')

[<matplotlib.lines.Line2D at 0x1f460145e80>]

Erasing the extrapolated points outside the sampled area!

# https://stackoverflow.com/questions/43892459/check-if-geo-point-is-inside-or-outside-of-polygon

# create the object Polygon from the 2D-Array
polygon = Polygon(tuple(map(tuple, mask)))

# create the object Point (the Albardão lighthouse)
point = Point(longitude, latitude)

print(polygon.contains(point)) # check if polygon contains point
print(point.within(polygon)) # check if a point is in the polygon 

False
False

'''
'Erasing' any data outside the boundaries!
'''
li, co = xx.shape

for l in range(li):
    for c in range(co):
        point = Point(xx[l,c], yy[l,c])
        
        if point.within(polygon) == False:
            bat[l, c] = np.nan

plt.rcParams.update({'font.size':15})
fig, ax = plt.subplots(figsize=(13,11))

conts = [10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 120, 140]
contour_plot = ax.contourf(xi, yi, bat, conts, cmap='rainbow')
ax.plot(junta[:,1], junta[:,0], '.k', markersize=.01)
ax.axis('equal')
plt.colorbar(contour_plot)

ax.plot(longitude, latitude, 'ro')
ax.set_ylim(-35.5, -31)
ax.set_xlim(-54, -49)
show(mapa_base, ax=ax)

ax.plot(secao1[:,1], secao1[:,0], '.', color='tab:red')
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')
ax.set_title('Bathymety (m)')

Text(0.5, 1.0, 'Bathymety (m)')

'''
Distrance from the Albardao Lighthouse
'''
dx = secao1[:,1] - longitude
dy = secao1[:,0] - latitude

dist = (dx**2 + dy**2)**.5 * 111.120

fig, ax = plt.subplots(figsize=(15,5))

ax.plot(dist, -secao1[:,2])
ax.set_ylabel('Profundidade (m)')
ax.set_xlabel('Distância (km)')
ax.set_ylim(-125, 0)

(-125.0, 0.0)

Digitizing the boundaries! The polygon that contains the data!

1) We need to create a polygon to delimit the area! We can use the 'ginput'... it's a little bit tricky!

https://stackoverflow.com/questions/41403406/matplotlib-use-of-ginput-on-jupyter-matplotlib-notebook

a) activate the external plotting with '%matplotlib qt'
b) enable the 'plt.ginput' to capture points! setup it accordingly!
c) the points will be stored as tuples!
d) if you succeed, save the captured points (pickle) de deactivate the 'qt' to back to work 'on line'
e) this is better made in separated notebook to avoid confusion... if you forget to de-activate the qt and ginput, the notebook will freeze until you use the ginput!

%matplotlib inline

'''
Enable to allow to use 'qt'
'''
# import matplotlib
# matplotlib.use('TkAgg')

"\nEnable to allow to use 'qt'\n"

'''
Enable to digitize the points to create the mask
'''
# %matplotlib qt
# fig, ax = plt.subplots(figsize=(13,11))

# conts = [10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 120, 140]
# contour_plot = ax.contourf(xi, yi, bat, conts, cmap='rainbow')
# ax.plot(junta[:,1], junta[:,0], '.k', markersize=.01)
# ax.axis('equal')
# plt.colorbar(contour_plot)

# ax.plot(longitude, latitude, 'ro')
# ax.set_ylim(-35.5, -31)
# ax.set_xlim(-54, -49)
# show(mapa_base, ax=ax)

# pts = plt.ginput(n=-1, timeout=-1)

'\nEnable to digitize the points to create the mask\n'

'''
Enable to save the points 
'''
# x = np.array(pts)
# with open('bathymetry_mask.pkl', 'wb') as io:
#     pickle.dump(x, io)

'\nEnable to save the points \n'