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');

# get an unic ID based on the timestamp
# and generate a directory where to store data products
ID = getID('morfo')
ensure_dir(ID)

setGenv(location='lonlat')

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'))

#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'))

HTML(region('dtm1')[0])

HTML(rasterinfo('dtm1')[0])

elev = 'dtm1'
sf = 5
cl = 'elevation_cl5'
clf='elevation_cl5nn5'
vpoint='site1'

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)

display_html('<iframe src="http://geofemengineering.it/data/img.html" width=750 height=650> </iframe>', raw=True)

clusterize(img[3:], cl, k=5, ss=10, ns=500)
Image(filename=cl+'.png')

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')

!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')

print clf3d
Image(filename=clf3d+'.png')

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')

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')

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')

#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)

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

htmlslide