0

Any way to read bathymetry surveys (*.xyz) into QGIS as a TIN or as raw points?

river

  • 2
    What is raw multibeam data? What hardware captures it? What do you usually use to read it? – HeikkiVesanto Feb 22 '17 at 9:48
  • 1
    What have you tried with QGIS so far? Where did you get stuck? You are welcome to discuss the state of GIS in the offshore oil and gas industry in the GIS Chat Room but comments about that and question processing within your question seem to be a distraction from the question you are asking, and so have been removed. – PolyGeo Feb 22 '17 at 9:49
  • 1
    What format is your data currently in? A raster could well work: do you have x, y, z data and some additional attribute(s) for each point? A multi-band raster could have the z in one band and additional attributes stored in additional bands. What are you trying to do with the data? What do you need to show? – George of all trades Feb 22 '17 at 12:49
  • This is one of the posts I mentioned in my previous comment: gis.stackexchange.com/questions/137518/… – Andre Silva Feb 22 '17 at 15:43
  • 1
    Perhaps you should cut your update and post it as an answer, so to preserve the Q/A format. – Andre Silva Mar 10 '17 at 11:23
5

Most GIS software is not really built for dealing with the amounts of data present in both raw multibeam data and raw lidar data (the two datatypes are quite similar in content and topic).

GIS software struggles when showing many points because no filtering is applied in an effort to reduce the number of individual returns shown.

All in all, you need to preprocess the data in a bit of dedicated software to reduce the point density and potentially convert it into a more suitable format, such as a raster.

2

The term "bathymetry" originally referred to the ocean's depth relative to sea level, although it has come to mean “submarine topography,” or the depths and shapes of underwater terrain.

In the same way that topographic maps represent the three-dimensional features (or relief) of overland terrain, bathymetric maps illustrate the land that lies underwater. Variations in sea-floor relief may be depicted by color and contour lines called depth contours or isobaths.

Bathymetry is the foundation of the science of hydrography, which measures the physical features of a water body. Hydrography includes not only bathymetry, but also the shape and features of the shoreline; the characteristics of tides, currents, and waves; and the physical and chemical properties of the water itself.

[https://oceanservice.noaa.gov/facts/bathymetry.html]

Solution 1: Import using GRASS

The r.in.xyz module will load and bin ungridded x,y,z ASCII data into a new raster map. The user may choose from a variety of statistical methods in creating the new raster. Gridded data provided as a stream of x,y,z points may also be imported.

r.in.xyz is designed for processing massive point cloud datasets, for example raw LIDAR or sidescan sonar swath data. It has been tested with datasets as large as tens of billion of points (705GB in a single file).

Source: https://grass.osgeo.org/grass73/manuals/r.in.xyz.html Also see this post for more info: Visualizing a LiDAR point cloud in 3D with GRASS?

Solution 2: LAStools

First convert the ascii xyz data into a las file (compressed, .laz), in ubuntu terminal. Then create a dem from the laz file, in this case stored as an asc raster. [Make] a hillshade (in this case geotiff) is also nice, you may do it for all the laz files in a directory.

Source: http://hydrogeotools.blogspot.no/2013/11/gridding-interpolate-xyz-data.html

Solution 3: Directly in QGIS

Add your xyz ascii file as a vector layer by "add delimited text layer". Then use "Interpolation" plugin to create a raster (from TIN or IDW) out of your delimited text layer (or use Raster-Analysis-Grid for more options)

Source: http://hydrogeotools.blogspot.no/2013/11/gridding-interpolate-xyz-data.html

Solution 4: Python and matplotlib

# -*- coding: utf-8 -*-
"""
__date__    2013-11-16
__author__  josef 
__web__     http://hydrogeotools.blogspot.se/
To import xyz data from an ascii file, interpolate and save as geotiff
"""
import numpy as np
import matplotlib.mlab as ml
import scipy.interpolate as il #for method2, in case the matplotlib griddata method fails
from osgeo import gdal
from osgeo import osr

fil_in = r"""/PathToFile/FileName.xyz""" #CHANGE HERE
raster_ut = r"""/PathToFile/RasterOut.tif""" #CHANGE HERE

x,y,z = np.loadtxt(fil_in, skiprows=1, delimiter=" ",unpack = True) #CHANGE HERE
xmin,xmax,ymin,ymax = [min(x),max(x),min(y),max(y)]

#size of 1 m grid
nx = (int(xmax - xmin + 1))#CHANGE HERE
ny = (int(ymax - ymin + 1))#CHANGE HERE

# Generate a regular grid to interpolate the data.
xi = np.linspace(xmin, xmax, nx)
yi = np.linspace(ymin, ymax, ny)
xi, yi = np.meshgrid(xi, yi) 

# Interpolate the values of z for all points in the rectangular grid
# Method 1 - Interpolate by matplotlib delaunay triangularizatio and nearest neigh. PLEASE NOTE! THIS FAILS QUITE OFTEN (http://matplotlib.org/api/mlab_api.html#matplotlib.mlab.griddata) But there might be a solution - install mpl_toolkits.natgrid (http://matplotlib.org/mpl_toolkits/)
zi = ml.griddata(x,y,z,xi,yi,interp='nn') #interpolation is 'nn' by default (natural neighbour based on delaunay triangulation) but 'linear' is faster (see http://matplotlib.1069221.n5.nabble.com/speeding-up-griddata-td20906.html)
# PLEASE NOTE! Method 1 fails sometimes and then using mpl_toolkits.natgrid may be a solution (http://matplotlib.org/api/mlab_api.html#matplotlib.mlab.griddata) (http://matplotlib.org/mpl_toolkits/)

# Otherwise, try Method 2 - Interpolate  using scipy interpolate griddata
#zi = il.griddata((x, y), z, (xi, yi),method='linear') #(may use 'nearest', 'linear' or 'cubic'  - although constant problems w linear) 

#---------------  Write to GeoTIFF ------------------------
nrows,ncols = np.shape(zi)
xres = (xmax-xmin)/float(ncols)
yres = (ymax-ymin)/float(nrows)
geotransform=(xmin,xres,0,ymin,0, yres) 

output_raster = gdal.GetDriverByName('GTiff').Create(raster_ut,ncols, nrows, 1 ,gdal.GDT_Float32,['TFW=YES', 'COMPRESS=PACKBITS'])  # Open the file, see here for information about compression: https://gis.stackexchange.com/questions/1104/should-gdal-be-set-to-produce-geotiff-files-with-compression-which-algorithm-sh
output_raster.SetGeoTransform(geotransform)  # Specify its coordinates
srs = osr.SpatialReference()                 # Establish its coordinate encoding
srs.ImportFromEPSG(3010)                     # This one specifies SWEREF99 16 30
output_raster.SetProjection( srs.ExportToWkt() )   # Exports the coordinate system to the file
output_raster.GetRasterBand(1).WriteArray(zi)   # Writes my array to the raster

Source: http://hydrogeotools.blogspot.no/2013/11/gridding-interpolate-xyz-data.html

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.