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
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]:
| projection | 1 (UTM) |
| zone | 19 |
| datum | wgs84 |
| ellipsoid | wgs84 |
| north | 4718670 |
| south | 4716550 |
| west | 392840 |
| east | 397080 |
| nsres | 10 |
| ewres | 10 |
| rows | 212 |
| cols | 424 |
| cells | 89888 |
In [7]:
HTML(rasterinfo('dtm1')[0])
Out[7]:
| north | 4718670 |
| south | 4716550 |
| east | 397080 |
| west | 392840 |
| nsres | 10 |
| ewres | 10 |
| rows | 212 |
| cols | 424 |
| cells | 89888 |
| datatype | FCELL |
| 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
Next notebooks ...¶
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 [ ]: