CREATING AN ESRI SHAPEFILE WITH BATHYMETRY CONTOURS FROM POINT DATA USING QGIS
This article shows how to use QGIS operations to convert Excel bathymetry point data to depth contours. When bathymetry readings are provided from measurements at many points in a body of water, the goal is to create polygons around all points with depth classes, and remove land areas using clipping.
I was given an Excel file containing over 200,000 points of data:
This much data is too much for QGIS to handle. For a contour map, I also didn’t need that much detail. Data can be simplified by keeping (for example) only every 100th row, reducing it to 2,000 data points. In Excel, this can be done using the following formula:
Creating this formula in a new sheet, and copying it by dragging down 2,000 times produces a new worksheet with simplified data.
Adding depth categories
The contours we will be creating need discrete depth classes. In my Excel data, there are columns
depth, where depths are anywhere between 0m and 35m, with up to 6 significant digits. For contour polygons, this level of detail is not desirable and depths should be expressed in classes (0-5m, 5-10m, etc.) To do this, Excel offers the
This produces multiples of 5 (i.e. 0, 5, 10, 15, …) which will be the classes. We will place them in a new column called
class, so that we get:
The resulting file (with the additional
class column) can now be exported to CSV.
QGIS can read point data from a CSV file using
Layer -> Add Layer -> Add delimited text layer. Default options work.
The input list of points must be converted to a set of polygons, all of which must have shared edges. One way to do this is to create Voronoi polygons. This creates a polygon around each point which edges equidistant to each neighbouring point.
QGIS offers Voronoi diagram creation under
Vector -> Geometry tools -> Voronoi polygons. The Voronoi polygons tool has an option to set a buffer region. This is the amount by which the resulting polygons will extend beyond the perimeter points. Setting it to 50% will help clipping later.
For 2,000 points this works well, while 20,000 points hangs QGIS. This appears to be a limitation of QGIS’s Python implementation of the Voronoi algorithm, since algorithms exist that can handle many more points and produce results in seconds - but using them would require leaving QGIS.
The Voronoi diagram algorithm will keep the
depth fields. These can be removed in the QGIS Attribute Table. Only the
class field need be kept.
Polygons with the same depth class may be joined (dissolved). QGIS offers this under
Vector -> Geoprocessing Tools -> Dissolve. The
class field must be selected to dissolve on. Dissolve all must be unselected.
The dissolved polygons must be clipped with the area originally occupied by the input points. This can be done using
Vector -> Geometry Tools -> Delaunay Triangulation, followed by a Dissolve on all fields, which produces a single polygon.
Vector -> Geoprocessing Tools -> Clip will clip the dissolved Voronoi diagram with the clip area.
Clipping may cause the
class field to become of string type. To change it back to integer, use
Refactor fields in the QGIS toolbox. This will create a new layer.
Note that Delaunay triangulation creates a convex shape. For concave shapes, it is better to use the Concave hull algorithm found in the QGIS toolbox - though there is no true concave hull and this is merely an approximation.
The clipped result will be blocky, consisting of clearly visible Voronoi cells. This can be improved using the GRASS toolkit, with the
v.generalize.smooth tool. Various methods are available, giving different results:
The Snakes method appears to produce rounded polygons suitable for depth contours. The output shown is for default parameters.
The resulting smoothed layer can now be saved as an ESRI shapefile.