Seafloor Habitat Mapping

Massimo Di Stefano

Software Engineer

Rensselaer Polytechnic Institute, Troy, NY

Woods Hole Oceanographic Institution, Woods Hole, MA

HABCAM

WorkfloW

In this notebook :

In [1]:
import os
import sys
import grass.script as grass
from ecoop.gis import setGenv, region, reproj, makemorfo, morf2png, tif2kml2, makeGE, gearth, rasterinfo, clusterize, makepng
from IPython.core.display import Image, HTML, display_html
from ecoop.esrutil import getID, getResults, ensure_dir, getZip, uploadfile
from ecoop.img2slider import img2slide1
import locale
# problem with my Italian ;)

locale.setlocale(locale.LC_ALL, 'en_US');
In [2]:
# get an unic ID based on the timestamp
# and generate a directory where to store data products
ID = getID('morfo')
ensure_dir(ID)
your session ID is :  morfo_Friday_28_June_2013_09_45_56_PM
  • Import Data ####GRASS working environment : #####World Geodetic System EPSG:4326
In [3]:
setGenv(location='lonlat')
{'MAPSET': 'PERMANENT', 'GISDBASE': '/Users/epi/grass7data', 'LOCATION_NAME': 'lonlat'}
In [4]:
grass.run_command('r.in.gdal', input='stell_10m.tiff', output='stell10m', flags='oe', overwrite=True)
grass.run_command('g.region', rast='stell10m', flags='ap', res=0.0001)
grass.run_command('r.mapcalc', expression="stell_na10m=if(abs(stell10m) < 198.3, stell10m)", overwrite=True)
grass.run_command('r.null', map='stell_na10m', setnull=0)
grass.run_command('v.in.ogr', dsn='kdata/points/site1.shp', output='site1', flags='oe', overwrite=True)
grass.run_command('g.region', vect='site1', flags='ap', res=0.0001)
grass.run_command('v.in.region', type='area', output='site1area', overwrite=True)
grass.run_command('r.out.png', input='stell_na10m', output=os.path.join(ID,'stell_na10m.png'), overwrite=True)
Image(filename=os.path.join(ID,'stell_na10m.png'))
Out[4]:
  • Reproject subset to UTM ####GRASS working environment : #####Universal Transverse Meractor - WGS84 fuse 19 North EPSG: 32619
In [5]:
#grass.run_command('g.proj', epsg='32619', location='utm19N_wgs84')
setGenv(location='utm19N_wgs84')
#grass.run_command('v.proj', input='site1area', location='lonlat', mapset='PERMANENT', overwrite=True)
#grass.run_command('v.proj', input='site1', location='lonlat', mapset='PERMANENT', overwrite=True)
grass.run_command('g.remove', rast='dtm1')
reproj(reg='site1', inp='stell_na10m', out='dtm1', res=10, buff=0, location='lonlat', mapset='PERMANENT')
grass.run_command('r.fillnulls', input='dtm1', output='dtm1', overwrite=True)
grass.run_command('r.out.png', input='dtm1', output=os.path.join(ID,'dtm1.png'), overwrite=True)
Image(filename=os.path.join(ID,'dtm1.png'))
{'MAPSET': 'PERMANENT', 'GISDBASE': '/Users/epi/grass7data', 'LOCATION_NAME': 'utm19N_wgs84'}
Out[5]:

Check extension and resolution

In [6]:
HTML(region('dtm1')[0])
Out[6]:
projection1 (UTM)
zone19
datumwgs84
ellipsoidwgs84
north4718670
south4716550
west392840
east397080
nsres10
ewres10
rows212
cols424
cells89888
In [7]:
HTML(rasterinfo('dtm1')[0])
Out[7]:
north4718670
south4716550
east397080
west392840
nsres10
ewres10
rows212
cols424
cells89888
datatypeFCELL
min-133.76
max-56.76
title (dtm1)
units"none"
vertical_datum"none"
timestamp"none"

Morphological Features Extraction

In [8]:
elev = 'dtm1'
sf = 5
cl = 'elevation_cl5'
clf='elevation_cl5nn5'
vpoint='site1'
In [9]:
img = makemorfo(input=elev,nnwin=9, pmwin=15, resolution=10, overwrite=True, remove=False) # makemorfo(input=layer, overwrite=True)
pnglist = morf2png(img)
from ecoop.uploadimage import uploadimg
from ecoop.img2slider import img2slide2
ID = !pwd
#for i in pnglist[3:]:
#    uploadfile(inputfile=i,outputfile='/var/www/esr/%s' % i)
html = img2slide1(pnglist)
f = open('img.html','w')
f.write(html)
f.close()
htmlslide = 'img.html'
#uploadfile(inputfile=htmlslide,outputfile='/var/www/esr/%s' % htmlslide)
written  dtm1_avg.png
written  dtm1_min.png
written  dtm1_max.png
written  dtm1_er.png
written  dtm1_slope.png
written  dtm1_profc.png
written  dtm1_crosc.png
written  dtm1_minic.png
written  dtm1_maxic.png
written  dtm1_longc.png
In [10]:
display_html('<iframe src="http://geofemengineering.it/data/img.html" width=750 height=650> </iframe>', raw=True)

Unsupervised Spatial Clustering

R - Cluster package - Clara()

In [11]:
clusterize(img[3:], cl, k=5, ss=10, ns=500)
Image(filename=cl+'.png')
Loading required package: morfoclara
Loading required package: sp
Loading required package: XML

Attaching package: 'XML'

The following object is masked from 'package:tools':

    toHTML

GRASS GIS interface loaded with GRASS version: GRASS 7.0.svn (2013)
and location: utm19N_wgs84
Loading required package: rgdal
rgdal: version: 0.8-9, (SVN revision Unversioned directory)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 1.10.0, released 2013/04/13
Path to GDAL shared files: /usr/local/Cellar/gdal/HEAD/share/gdal
Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
Path to PROJ.4 shared files: (autodetected)
Out[11]:

Apply a nearest neighbors filter to remove smaller classes and export the results

In [12]:
grass.run_command('r.neighbors', flags='c', input=cl, output=clf, method='mode', size=sf, overwrite=True)
grass.run_command('g.region', rast=clf, flags='ap')
grass.run_command('r.out.gdal', input=clf, output=clf+'.tif', format='GTiff', type='Byte', flags='f', overwrite=True)
makepng(rasterlist=[str(clf)],
        vectorlist=[str(vpoint)], 
        output=clf+'.png')
Image(filename=clf+'.png')
png done
Out[12]:

3D Perspective View

In [13]:
!rm -rf dtm1_3d.png
!rm -rf {clf3d}+'.png'
grass.run_command('m.nviz.image', elevation_map='dtm1', vpoint=vpoint, vpoint_size=30, output='dtm1_3d', perspective=15, height=2000, color_map='dtm1', resolution_fine=1, resolution_coarse=1, format='tif')
os.system('convert dtm1_3d.tif dtm1_3d.png')
clf3d = clf+'_3d'
grass.run_command('m.nviz.image', elevation_map=elev, vpoint=vpoint, vpoint_size=30, output=clf3d, perspective=15, height=2000, color_map=clf, resolution_fine=1, resolution_coarse=1, format='tif')
instruction = 'convert %s.tif %s.png' % (clf3d, clf3d)
os.system(instruction)
Image(filename='dtm1_3d.png')
Out[13]:
In [14]:
print clf3d
Image(filename=clf3d+'.png')
elevation_cl5nn5_3d
Out[14]:

Vector - Raster Overlay

In [15]:
grass.run_command('v.db.addcolumn', map=vpoint, columns="depth double precision,morfo int")
grass.run_command('v.what.rast', map=vpoint, rast=clf, col='morfo')
grass.run_command('v.what.rast', map=vpoint, rast=elev, col='depth')
Out[15]:
0

Tiles export & visualization in Cesium

In [16]:
elev = 'dtm1'
sf = 5
cl = 'elev1cl5'
clf='elev1cl5nn5'
vpoint='site1'
grass.run_command('g.region', rast=clf, flags='ap')
grass.run_command('r.out.gdal', input=clf, output=clf+'.tif', format='GTiff', type='Byte', overwrite=True, flags='f')
Out[16]:
0
In [17]:
import envoy
instruction = []
inp, out = clf+'.tif', clf
instruction.append('gdalbuildvrt %s_tmp.vrt %s' % (out, inp))
instruction.append('gdalwarp -of VRT -t_srs EPSG:4326 %s_tmp.vrt %s_geo.vrt' % (out, out)) 
instruction.append('gdal_translate -of vrt -expand rgba %s_geo.vrt %s.vrt' % (out, out))
instruction.append('gdal2tiles.py --zoom=0-12 -p geodetic -k %s.vrt' % out)
for i in instruction:
    print i
    r = envoy.run(i, timeout=12)
    print r.std_out.strip().split('\n')
gdalbuildvrt elev1cl5nn5_tmp.vrt elev1cl5nn5.tif
['0...10...20...30...40...50...60...70...80...90...100 - done.']
gdalwarp -of VRT -t_srs EPSG:4326 elev1cl5nn5_tmp.vrt elev1cl5nn5_geo.vrt
['']
gdal_translate -of vrt -expand rgba elev1cl5nn5_geo.vrt elev1cl5nn5.vrt
['Input file size is 447, 169']
gdal2tiles.py --zoom=0-12 -p geodetic -k elev1cl5nn5.vrt
['Generating Base Tiles:', '0...10...20...30...40...50...60...70...80...90...100 - done.', 'Generating Overview Tiles:', '0...10...20...30...40...50...60...70...80...90...100 - done.']
In [18]:
#cesium = makeCesium(template=default, tiles=clf)
# the data are embedded into the cesium widget
# you should zoom and pan to the Cape Code area (North East US, south of Boston)
display_html('<iframe src="http://geofemengineering.it/widget/globes/data.html" width=800 height=600> </iframe>', raw=True)

Share it!

In [ ]:
lines = !jist -p Seafloor_Habitat_Mapping.ipynb
print lines[0].replace("https://gist.github.com", "http://nbviewer.ipython.org")
In [ ]:
htmlslide
In [ ]: