Bathymetry profile with QGIS and Python
In my case I have a set of locations in the Firth of Thames, an enclosed bay located in the north-east coast of the North Island of New Zealand. The goal is to generate a plot of the bathymetry profile between my locations. To do so the first thing one needs is a bathymetric dataset. The National Institute of Water and Atmospheric Research (NIWA) has several datasets available on their website.
For this example I'm using the Hauraki Gulf Bathymetry set (20 m resolution) from Mackay et al. (2012). The LINZ data service hosts a nice array of datasets as well. Once you've downloaded it, proceed to open it with QGIS.
The points of my transect stored as a CSV file:
X,Y,id 1806344.08570028,5882018.3007867,1 1812868.47817093,5888135.57116718,2 1813249.50269121,5894221.5244638,3 1812236.91697977,5906351.67494522,4
Add the points file as a "Delimited Text Layer". If you're using a basemap like Google Maps (Satellite) layer you'll see something like:
Map of the Firth of Thames and the transect locations
Before proceeding make sure you have the Profile Tool plugin installed. Click on the plugin icon (Terrain profile). A panel will appear on the bottom, just below your map. Click on "Add Layer" and select your bathymetric layer, "hg_dtm_noland" in my case. Select your points layer in Layers panel and in the Options choose "Selected layer". Magic! You should see something remotely similar to:
Screenshot of QGIS with the Profile tool panel in the bottom.
You can save your profile as is or visualize it with Matplotlib in the "Settings" tab. The plotting options are limited though so in my case I exported the profile data so I could generate a plot from scratch. Go to the "Table" tab and copy the data to clipboard, with coordinates. Paste it in a file and save it as "transect_bathymetry.csv" replacing blank spaces with commas. Each record contains the X (distance in m), the Y (elevation in m) and the coordinates.
The following code reads the freshly exported data and saves the plot as a PNG file. It's also available as a Github gist. The code:
import numpy as np import matplotlib.pyplot as plt import pandas as pd import seaborn as sns data_path = './transect_bathymetry.csv' df = pd.read_csv(data_path) df = df.apply(pd.to_numeric) # sites distance in m sites = [ 0.0, 3630.692532275101, 8490.190136457848, 25767.366984381923 ] dfsites = df.loc[df['x'].isin(sites)] sns.set(rc={'figure.figsize': (11, 4)}) fig, ax = plt.subplots(nrows=1, ncols=1) plt.yticks(np.arange(-20, 7, 4)) ax.plot((df.x/1000), df.y) plt.stem(dfsites.x/1000, dfsites.y, bottom=-20, linefmt='C0--', basefmt=' ') plt.xlabel("Distance (km)") plt.ylabel("Elevation (m +MSL)") ax.axis([-0.5, 27, -20, 6]) # x and y limits ax.set_xticks([s/1000 for s in sites], minor=True) # Sites labels y = 0 for i, r in dfsites.iterrows(): plt.annotate("S%i" % (y + 1), xy=(r["x"]/1000, r["y"]), xytext=(r["x"]/1000 + 0.1, r["y"] + 0.5), size="15") y += 1 fig.tight_layout() plt.savefig("./profile.png", bbox_inches='tight', dpi=200) plt.show() plt.close()
And the result:
Bathymetry profile of sites 1 to 5
The more resolution the bathymetry/elevation dataset you use the more accurate the final plot will be. The dataset I used has a resolution of 20m which is nice but could be better.
References
Mackay, K.A., Mackay, E.J., Neil, H.L., Mitchell, J.S., Bardsley, S.A. (2012). Hauraki Gulf NIWA Chart, Miscellaneous Series 91. ftp://ftp.niwa.co.nz/bathymetry/HGBathy_DTM_binary_grid.zip