Grid Averaging and Plotting Bathymetry Data with Python¶
Ok, so you built your own autonomous surface vehicle (ASV), integrated all the components, got all your sensors working, and you’ve finally collected some bathymetry data. Now what? ReefMaster offers a free trial, but what if you don’t have $200? Fear not, we are going to grid your bathymetry data here in python.
Python Libraries¶
Here, we are using Python 3. The code in this notebook requires the following Python libraries.
Example Data¶
The example data used here were provided by …
import numpy as np
from numpy import genfromtxt
First, import numpy, as we will use it later. If you’re new to Python or numpy, check out the documentation
We will use numpy’s genfromtxt to read in the csv data. In our case, our data consists of four columns with values separated by commas. The first line of the data file is shown below. The columns are latitude, longitude, depth, and some unknown value
22.22917557,113.8987656,45.834,1.59299E+12
…
Remember to change the data_dir to point to the file on your local computer.
data_dir = '/Users/dfc/Documents/Other/MakerBoat/python/depth/'
file_name = data_dir + 'example_survey_data.csv'
bath_data = genfromtxt(file_name, delimiter=',')
Now that we have used genfromtxt to read in the data from the csv file, let’s check the type of the variable bath_data
type(bath_data)
numpy.ndarray
print(bath_data.ndim)
2
print(bath_data.shape)
(3998, 4)
So we can see that our variable type is numpy’s ndarray, it has two dimensions, and the shape consists of nearly 4000 rows and 4 columns.
Now let’s plot the lat/lon positions to make sure that they make sense.
lat = bath_data[:,0]
lon = bath_data[:,1]
depth = bath_data[:,2]
import matplotlib.pyplot as plt
plt.title("Bathy Measurement Locations")
plt.xlabel("Longitude")
plt.ylabel("Latitude")
plt.plot(lon,lat,'.')
dL = 0.0005
plt.xlim([lon.min()-dL, lon.max()+dL])
plt.ylim([lat.min()-dL, lat.max()+dL])
plt.show()
Ok, we can see what appears to be a survey grid, but the x-ticks and y-ticks look a bit funky. Let’s clean that up using matplotlib‘s ticker
import matplotlib.ticker as mtick
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(lon,lat,'.')
ax.yaxis.set_major_formatter(mtick.FormatStrFormatter('%.4f'))
ax.xaxis.set_major_formatter(mtick.FormatStrFormatter('%.4f'))
plt.xlabel("Longitude")
plt.ylabel("Latitude")
plt.show()
Now, a quick scatter plot of the depth value at each location
from matplotlib import cm
fig = plt.figure(figsize=(20,6))
ax = fig.add_subplot(131)
sc = ax.scatter(lon,lat,s=20,c=depth, marker = 'o', cmap = cm.jet );
ax.yaxis.set_major_formatter(mtick.FormatStrFormatter('%.4f'))
ax.xaxis.set_major_formatter(mtick.FormatStrFormatter('%.4f'))
plt.xlabel("Longitude")
plt.ylabel("Latitude")
plt.colorbar(sc)
plt.show()
from pyproj import Proj
myProj = Proj("+proj=utm +zone=49N, +north +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
UTMx, UTMy = myProj(lon, lat)
UTMx.max() - UTMx.min()
231.75115337385796
UTMy.max() - UTMy.min()
319.2370445705019
dX = 10
X = np.arange(UTMx.min(), UTMx.max(), dX)
Y = np.arange(UTMy.min(), UTMy.max(), dX)
grid_x, grid_y = np.meshgrid(X, Y, sparse=False, indexing='xy')
from scipy.interpolate import griddata
zi = griddata((UTMx, UTMy), depth, (grid_x, grid_y),method='linear')
print(grid_x.shape)
(32, 24)
print(zi.shape)
(32, 24)
fig = plt.figure(figsize=(20,6))
ax = fig.add_subplot(121)
cp = plt.contourf(grid_x, grid_y, zi)
plt.show()
lon_spc = np.linspace(lon.min(), lon.max(), 40)
lat_spc = np.linspace(lat.min(), lat.max(), 40)
lon_grid, lat_grid = np.meshgrid(lon_spc, lat_spc, sparse=False, indexing='xy')
zL = griddata((lon, lat), depth, (lon_grid, lat_grid),method='cubic')
fig = plt.figure(figsize=(30,6))
ax = fig.add_subplot(141)
cp = plt.contourf(lon_grid, lat_grid, zL, cmap = cm.jet)
ax.yaxis.set_major_formatter(mtick.FormatStrFormatter('%.4f'))
ax.xaxis.set_major_formatter(mtick.FormatStrFormatter('%.4f'))
plt.xlabel("Longitude")
plt.ylabel("Latitude")
plt.colorbar(cp)
plt.show()
