Marine Science
this page is a work in progress
Tools for marine scientists
Bathymetry processing
Please expand
- Basic DEM creation: v.in.ascii + v.surf.rst
- r.in.xyz for processing multibeam sonar swaths
- r.surf.nnbathy add-on script for creating bathymetic DEMs from input x,y,z data
Bathymetric data
- See the Global datasets wiki page for more datasets (ETOPO, GEBCO, etc.)
- Smith and Sandwell 1-minute global elevation v10.1, May 13, 2008
global_topo_1min/README_V10.1.txt file:
Version 9.1 has a very different FORMAT than V8.2 The main differences are that the grid spacing in longitude is now 1 minute rather than 2 minutes. In addition, the latitude range is increased to +/- 80.738. Like the old versions, the elevation(+) and depth(-) are stored as 2-byte integers to the nearest meter. An odd depth of say -2001m signifies that this pixel was constrained by a real depth sounding while an even depth of say -2000m is a predicted depth. Here are the parameters for the old and new versions: param V8.2 V9.2 ___________________________ nlon 10800 21600 nlat 12672 17280 rlt0 -72.006 -80.738 rltf 72.006 80.738 ___________________________ The binary format of the integers is bigendian so the bytes need to be swapped if you are running on an Intel processor. Here is a typical command for swapping bytes: dd if=topo_9.1.img of=topo_9.1.img.swab bs=21600 conv=swab.
- GMT's img2grd + grd2xyz shows FP elevation values to the nearest cm not meter. Are these from contributed datasets? How does that fit with the odd/even real/interpolated soundings?
Import using GMT
- see also the GRASS and GMT wiki page.
Process with GMT's img2grd to convert from spherical Mercator projection to geographic coordinates, then import into GRASS
img2grd topo_10.1.img -T1 -S1 -V -R0/360/-80.738/80.738 -m1 -D -Gtopo_all.grd
# (out of memory, needs 1.4gb)
# try just for NZ (W/E/S/N bounds)
REGION=160/180/-50/-30
img2grd topo_10.1.img -T1 -S1 -V -R"$REGION" -m1 -D -Gtopo_NZ.grd
grd2xyz topo_NZ.grd -S > topo_NZ.xyz
# get adjusted region bounds and resolution from img2grd output
# ** check that rows and columns match **
g.region n=-29.9945810754 s=-50.0056468984 w=160E e=180 \
ewres=0:01 nsres=0.0126094 -p
r.in.xyz in=topo_NZ.xyz out=topo_NZ_1min fs=tab
r.colors output=topo_NZ_1min color=etopo2
To save a step or some disk space, in the above you could set the region first then pipe grd2xyz directly into r.in.xyz instead of creating the .xyz file.
# create a r.in.xyz "n" map to test input point coverage r.in.xyz in=topo_NZ.xyz out=topo_NZ_1min_n fs=tab method=n # check rast map stats, min=max=1 and there should be no null cells r.univar topo_NZ_1min_n # cleanup g.remove topo_NZ_1min_n
or, import GMT .grd file directly (old GMT grd format introduces FP +0.005 elev shift error??). New GMT netCDF format .grd files can be imported with the r.in.gdal module.
# convert COARDS-compliant netCDF grdfile to old GMT native .grd grdreformat topo_NZ.grd topo_NZ_old.grd=bf # import r.in.bin -hf in=topo_NZ_old.grd out=topo_NZ_old
Import directly
To load it into GRASS lat/lon location (spherical):
- Location setup:
- http://thread.gmane.org/gmane.comp.gis.gmt.user/918
- http://article.gmane.org/gmane.comp.gis.proj-4.devel/192/
Is it even possible to load directly into GRASS?
Set up Mercator/Sphere location:
- g.setproj commands for manual projection settings
Projection type> D "other" proj> merc No datum ellipsoid> sphere radius> default (doesn't matter) Scale Factor> 1.0 Latitude of True Scale> 0 Central Meridian> 0
Which creates:
G63> g.proj -j
+proj=merc
+k_0=1.0000000000
+lat_ts=0.0000000000
+lon_0=0.0000000000
+a=6370997
+b=6370997
+no_defs
+to_meter=1.0
G63> g.proj -w
PROJCS["Mercator",
GEOGCS["unnamed",
DATUM["unknown",
SPHEROID["unnamed",6370997,"inf"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]],
PROJECTION["Mercator_2SP"],
PARAMETER["standard_parallel_1",0],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",0],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["meter",1]]
- Some web searches show people using the MRWORLD projection??
- "MRWORLD" also seen in ftp://topex.ucsd.edu/pub/global_topo_1min/topo_9.1b.img.ers
- gdal's ermapper support files has params: (/usr/share/gdal/ecw_cs.dat)
MRWORLD:PROJCS["unnamed",PROJECTION["Mercator_1SP"], PARAMETER["latitude_of_origin",0], PARAMETER["central_meridian",0], PARAMETER["scale_factor",1], PARAMETER["false_easting",20000000], PARAMETER["false_northing",0]]
Note Mercator_1SP vs. Mercator_2SP in the above. (does 2 std parallels merc with only one defined == 1 std par merc?)
- Load using r.in.bin
# the following does not work correctly, just a trial
# offset n,s,e,w by 1/2 a grid cell?
r.in.bin input=topo_9.1b.img output=topo_9.1b \
title="1' worldwide relief (1.852 km-sq)" \
-b -s bytes=2 rows=17280 cols=21600 \
n=80.738 s=-80.738 w=0 e=360
r.colors output=topo_9.1b color=etopo2
Official coloring
Download the "official" GMT color rules from:
wget ftp://topex.ucsd.edu/pub/global_topo_1min/gmt_examples/map/topo.cpt
Convert HSV GMT cpt color rules to RGB GRASS color rules with the r.cpt2grass add-on script.
r.cpt2grass in=topo.cpt out=palette_topo.gcolors
(HSV -> RGB conversion in that script is now partially functional)
Multibeam sonar processing
- See the LIDAR and Multi-beam Swath bathymetry data wiki help page
- r.in.xyz can import multibeam sonar swaths from a raw x,y,z stream and bin them on the fly using statistical filters.
MB-System
- MB-System is Free software for the processing and display of swath and sidescan sonar data. It can handle both multibeam bathymetry and sidescan sonar image data.
Import into GRASS
- Gridded data can be loaded into GRASS either as a GMT NetCDF grid via the r.in.gdal module (see GDAL supported formats), old-style GMT grid using r.in.bin -h, or as an Arc ASCII grid via the r.in.arc module.
- See the GRASS and GMT wiki help page for more information.
- Ungridded data points may be piped directly from mblist to GRASS's v.in.ascii module. See the v.colors addon script for colorizing point data in GRASS (v.colors may be unsuitable for massive datasets).
Examples
- Export Lat/Lon + depth data from XTF datafile into a GRASS Lat/Lon location
mblist -I 074.XTF -OXYz | v.in.ascii out=track074 x=1 y=2 fs=tab
- Export Lat/Lon from the XTF datafile, reproject into the current GRASS location's projection, and import into GRASS with v.in.ascii
mblist -I 074.XTF -OXY | m.proj -i | cut -f1 -d' ' | \ v.in.ascii out=track074 x=1 y=2 fs=tab
Idea: write a v.in.cdl script that will parse a NetCDF/CDL file and automatically set v.in.ascii's column= option with column names and types.
- Read a series of .fbt pre-processed bathymetry files, reproject into the current projection, and grid into raster maps.
First create .fbt .fnv and .inf summary files
find . | grep '.[xX][tT][fF]$' > tmplist
mbdatalist -F-1 -I tmplist > datalist-1
mbdatalist -F-1 -I datalist-1 -N
Next create quick GMT plot of the nav lines to check coverage
# scan multi-dir for bounds
for lltype in Longitude Latitude ; do
for tb in head tail ; do
grep $lltype `find | grep '\.inf$'` | \
awk '{print $3 "\n" $6}' | sort -n | $tb -n1
done
done
EXTENT=167.68/169.09/40.48/40.82
mbdatalist -F-1 -I datalist-1 -R$EXTENT > survey-datalist
mbm_plot -F-1 -I survey-datalist -N
./survey-datalist.cmd
Based on the above PostScript file set GRASS region by eye:
#LL: n=40:48N s=40:37N w=168:40E e=169:06E
Convert to local projection with m.proj:
echo "168d40E 40d37N 169d06E 40d48N" | m.proj -i
Plug those numbers into g.region:
g.region s=4616845.02 n=4636836.33 w=520872.47 e=524862.14
Set the grid cell size to 5m, and align to whole numbers, and check that the rows x columns is reasonable (smaller than 40000x40000):
g.region res=5 -a -p
Adjust resolution to smaller grid size if needed (1000x1000 is fine for a summary image).
g.region res=10 -a -p
Scan, read, and reproject all .fbt data into a single x,y,z text file:
- note that this file can get very large, perhaps too large for a 32bit OS/filesystem. In these cases you can pipe directly into r.in.xyz so no file is written to disk.
- note that -Rw/e/s/n is used to skip over out-of-region data. MB-System and GMT will accept DDD:MM:SS.SSSh format as well as decimal degrees, just like GRASS.
- The following assumes the output is projected, for import into WGS84 lat/lon skip the cs2cs command or at least change "%.3f" to "%.8f".
OUTPROJ="`g.proj -jf`"
(
for FILE in `find . | grep '\.fbt$'` ; do
echo "Reading <$FILE> ..." 1>&2
mblist -I$FILE -D3 -R168:40E/169:06E/40:37N/40:48N
done
) | cs2cs -f '%.3f' +init=epsg:4326 +to $OUTPROJ | \
tr ' ' '\t' > fbt_UTM58_zoom.dat
Based on scan of the .inf files and personal knowledge set some reasonable bounds for the z-data:
# scan multi-dir for depth bounds
for tb in head tail ; do
grep Depth `find | grep '\.inf$'` | \
awk '{print $3 "\n" $6}' | sort -n | $tb -n1
done
#-2.7339
#200.6072
ZRANGE="-250,-2"
Run r.in.xyz to make aggregate raster maps. (Running 3 of them in the background to take full advantage of a quad-core CPU)
r.in.xyz fs=tab x=1 x=2 z=3 zrange=$ZRANGE \
in=zoom_fbt.dat out=zoom_fbt.mean method=mean --q &
r.in.xyz fs=tab x=1 x=2 z=3 zrange=$ZRANGE \
in=zoom_fbt.dat out=zoom_fbt.median method=median --q &
r.in.xyz fs=tab x=1 x=2 z=3 zrange=$ZRANGE \
in=zoom_fbt.dat out=zoom_fbt.trim20 method=trimmean trim=20 --q &
r.in.xyz fs=tab x=1 x=2 z=3 zrange=$ZRANGE \
in=zoom_fbt.dat out=zoom_fbt.n method=n
In my sample data, median and trim-mean 20% look the best.
Analyze "n" map for number of pings per square-meter.
- no-data cells not randomly distributed so we set n=0 to NULL so we aren't biasing the data as much
r.null zoom_fbt.n setnull=0 r.univar zoom_fbt.n
"mean:" is number of pings per square-meter for cells which have data.
Now analyze bathymetry data:
r.univar zoom_fbt.median
and set some nice colors:
r.colors zoom_fbt.median color=bcyr
Finally display it:
d.mon x0 d.rast zoom_fbt.median d.legend zoom_fbt.median d.barscale at=5,5
Mirone
- From the Google Code project description:
- "Mirone is a Windows MATLAB-based framework tool that allows the display and manipulation of a large number of grid/images formats through its interface with the GDAL library. Its main purpose is to provide users with an easy-to-use graphical interface to manipulate GMT grids. In addition it offers a wide range of tools dedicated to topics in the earth sciences, including tools for multibeam mission planning, elastic deformation studies, tsunami propagation modeling, earth magnetic field computations and magnetic Parker inversions, Euler rotations and poles computations, plate tectonic reconstructions, and seismicity and focal mechanism plotting. The high quality mapping and cartographic capabilities for which GMT is renowned is guaranteed through Mirone’s ability to automatically generate GMT cshell scripts and dos batch files."
You can interface with it via GDAL/GMT/netCDF formats, or directly transfer Matlab arrays with the r.out.mat and r.in.mat modules.
Sidescan sonar processing
- MB-System, as above.
- i.warp script for georectifying and mosaicking scanned paper rolls into a GeoTIFF with GDAL's gdalwarp program
- v.swathwidth script for planning swath bathymetry surveys
Wave exposure
- Using GRASS to prepare and process data for the SWAN Wave Model
Circulation models
- Preparing input grids
Tutorials
Remote Sensing
- Importing MODIS Aqua SST and chlorophyll-a data, SeaWiFS chlorophyll-a, and Pathfinder AVHRR SST satellite images.