david__init__~

#!/usr/bin/env python

import os
import sys
import subprocess
import numpy
import osgeo
import time
from osgeo import ogr

# Linux
grass7bin_lin = 'grass70'

# add your path to grassdata (the GRASS GIS database) directory
gisbase = '/usr/lib/grass70'

gisdb = os.path.join(os.path.expanduser("~"), "grassdata")

# specify (existing) location and mapset, program has to run once to set these up
location = "newLocation"
mapset   = "test"

#gisbase = 'C:\Program Files (x86)\GRASS GIS 6.4.5svn' # query GRASS 7 itself for its GISBASE
gisbase = '/usr/lib/grass70'

# Set GISBASE environment variable
os.environ['GISBASE'] = gisbase
# the following not needed with trunk
os.environ['PATH'] += os.pathsep + os.path.join(gisbase, 'extrabin')

# define GRASS-Python environment
gpydir = os.path.join(gisbase, "etc", "python")
sys.path.append(gpydir)

# Set GISDBASE environment variable
os.environ['GISDBASE'] = gisdb

# import GRASS Python bindings (see also pygrass)
import grass.script as gscript
import grass.script.setup as gsetup
import random
from grass.script import raster as grassR
from osgeo import ogr

#time variable for file names

timestr = time.strftime("%Y%m%d_%H%M%S")

#Consider adding time and date to name + maybe random integer
original = 'original' + timestr
ocean_point = 'ocean_point' + timestr
Xland = 'X_land' + timestr

#output lists
outputnames_vect=[]
outputnames_rast=[]


#ask for flood info
maxlevel_input = float(input("Please enter the maximum desired water level in meters: "))
interval_input = float(input("Please enter the desired flood intervals in meters: "))

loop_no = int(maxlevel_input/interval_input)

#area statistic
areaList=[]

areaPtg=[]

def main():

		#ask for coordinates of the ocean point
		#x_ocean = float(input("Please provide with the X coordinate of the ocean point: "))
		#y_ocean = float(input("Please provide with the Y coordinate of the ocean point: "))

		y_ocean = 43.201712

		x_ocean = 4.450673


#		y_ocean = 55.786290

#		x_ocean = 12.064200

		#initial loop number

		actual_loop = 0

		#point creation

		import os.path
		if os.path.isfile('test.geojson'):
			os.remove('test.geojson')

		create_point(x_ocean, y_ocean)

		print "You have chosen a point with the following coordinates: X %s, Y %s" % (x_ocean, y_ocean)

		# WGS84/initial grass location/mapset environment
		gsetup.init(gisbase, gisdb, location, mapset)

		#cleaning up the data

		cleanup_data()

		gscript.run_command('v.in.ogr',  flags='o', input='test.geojson', output=ocean_point)

		#importing the investigated DEM
		gscript.run_command('r.in.gdal', flags='', input ='eudem_N43E0041.tif', output=original)

		#Deleting  values keep
		b_area = 'b_area' + timestr
		gscript.run_command('r.mapcalc', expression= '%s = if(%s <=0, 0, 0)' % ('b_area', original))

		#converting the original tif into bounding vector polygon

		gscript.run_command('r.to.vect', input = 'b_area', output = 'b_polygon', type = 'area')

		#Extract the land from DEM
		gscript.run_command('r.mapcalc', expression= '%s = if(%s != 0, %s, null())' % (Xland, original, original))

		#Area statistic of the original DEM
		print gscript.run_command('r.stats', flags='aA', input=original)

		print "Import done"


		# Flood modeling
		while actual_loop <= loop_no:

			#Defining of the region of the DEM
			gscript.run_command('g.region', rast=original)

			expressionout = 'out_flood' + timestr + '_' + str(actual_loop)

			#Extracting ocean from the DEM

			gscript.run_command('r.mapcalc', expression= '%s = if(%s <= (%s * %s), 1, null())' % (expressionout, original, interval_input, actual_loop))

			print "Mapcalc done"

			# Vector point creation

			ocean_vector = 'ocean_vector_' + timestr + '_' + str(actual_loop)

			gscript.run_command('r.to.vect', input = expressionout, output = ocean_vector, type = 'area')

			print "Vector conversion done"

			#Waterrise as a vector output

			selected_ocean = 'ocean_select_' + timestr + '_' + str(actual_loop)

			gscript.run_command('v.select', ainput=ocean_vector, binput=ocean_point, output=selected_ocean, operator='intersects')

			intermediate = 'intermediate_'+ timestr + "_" + str(actual_loop)

			selected_ocean_rast = 'ocean_select_rast_' + timestr + '_' + str(actual_loop)

			gscript.run_command('v.to.rast', input=selected_ocean, output=intermediate, use='val', value=0)

			gscript.run_command('r.mapcalc', expression= '%s = if( %s == 0, (%s * %s) , null())' % (selected_ocean_rast, intermediate, interval_input, actual_loop))

			gscript.run_command('g.remove', flags='f', type = 'raster', pattern='inter*')

			outputnames_rast.append(selected_ocean_rast)

			outputnames_vect.append(selected_ocean)

			print "select done"
			import os.path
			if os.path.isfile('temp.txt'):
				os.remove('temp.txt')

			gscript.run_command('r.report', map=selected_ocean_rast, units='k', output='temp.txt')

			print selected_ocean_rast

			tempTxt = open('temp.txt', 'r')

			actual_increment = str(actual_loop*interval_input)


			for row in tempTxt:
				print row
				if row[1]== actual_increment[0]:
					AreaOut=row.split('|')
					areaList.append(AreaOut[3])

			tempTxt.close()

			old_loop = actual_loop - 1

			print areaList

			if float(areaList[actual_loop]) - float(areaList[old_loop]) > 0.02:

					gscript.run_command('v.overlay', ainput = selected_ocean, binput = 'old_flood', operator='not', output='subtructed_flood')

					gscript.run_command('v.to.db', map = 'subtructed_flood', option = 'area', columns = 'a_value')

					gscript.run_command('v.extract', input = 'subtructed_flood' , output = 'subtracted_flood_extract', where='a_value > 1000')

					gscript.run_command('v.to.rast', input = 'subtracted_flood_extract', output = 'subtracted_raster', use = 'cat')

					gscript.run_command('r.mapcalc', expression= '%s = if( isnull(%s), %s, %s)' % ('merged', 'subtracted_raster', 'old_flood_rast', selected_ocean_rast))

					gscript.run_command('r.neighbors', input = 'merged', output = 'diversity', method = 'diversity')

					pourpoint = 'pourpoint'+ "_" + str(actual_loop)

					gscript.run_command('r.mapcalc', expression= '%s = if(%s == 2, %s, null())' % (pourpoint, 'diversity','diversity'))

					gscript.run_command('r.to.vect', input = pourpoint, output = pourpoint, feature = 'point')

					gscript.run_command('g.remove', flags='f', type = 'raster', pattern='sub*')

					gscript.run_command('g.remove', flags='f', type = 'vector', pattern='sub*')

					gscript.run_command('g.remove', flags='f', type = 'raster', pattern='merged*')

					gscript.run_command('g.remove', flags='f', type = 'raster', pattern='diversity*')

			gscript.run_command('g.remove', flags='f', type = 'vector', pattern='old_flood')

			gscript.run_command('g.remove', flags='f', type = 'raster', pattern='old_flood_rast')

			gscript.run_command('g.copy', vect = (selected_ocean, 'old_flood'))

			gscript.run_command('g.copy', rast = (selected_ocean_rast, 'old_flood_rast'))


			actual_loop = actual_loop + 1

		print "Here we are"




		#creating watershed
		#watershed()

		#projection of the whole dataset to the required projected coordinate system
		#projection()

		#execute the calculatio in projected environment:
		#projected_calc()

#def area_analysis():






def projection():

	# selection of the projected coordinate system
	epsg_code=str(raw_input('Choose your Projected Coordinate System:'))

	# chosen name of the location where the PERMANENT mapset will be stored
	location_name=str(raw_input('Choose a name for your the location of your mapset(one word):'))

	# creation of the new location
	gscript.run_command('g.proj', flags='c', epsg=epsg_code, location=location_name)

	# changing to the new location and to its PERMANENT mapset
	gscript.run_command('g.mapset', mapset='PERMANENT', location=location_name) #% (epsg_code, gisdb))

	#projecting of the bounding polygon
	gscript.run_command('v.proj', input='b_polygon', location=location, mapset=mapset, dbase=gisdb)

	gscript.run_command('v.proj', input=ocean_point, location=location, mapset=mapset, dbase=gisdb)

	gscript.run_command('g.region', vect='b_polygon')

	gscript.run_command('r.proj', flags='n', input=original, location=location, mapset=mapset, dbase=gisdb, output=original, memory='800')

	print 'projecting is done'

def projected_calc():

	actual_loop_proj = 0
	# Flood modeling
	while actual_loop_proj <= loop_no:

		#Defining of the region of the DEM
		gscript.run_command('g.region', rast=original)

		expressionout = 'out_flood' + timestr + '_' + str(actual_loop_proj)

		expressionout_land = 'out_land' + timestr + '_' + str(actual_loop_proj)

		#Extracting ocean from the DEM

		gscript.run_command('r.mapcalc', expression= '%s = if(%s <= %s, 0, null())' % (expressionout, original, actual_loop_proj))

		#Extracting land from the DEM
		gscript.run_command('r.mapcalc', expression= '%s = if(%s > %s, 1, null())' % (expressionout_land, original, actual_loop_proj))

		print "Mapcalc done"

		# Vector point creation

		ocean_vector = 'ocean_vector_' + timestr + '_' + str(actual_loop_proj)

		gscript.run_command('r.to.vect', input = expressionout, output = ocean_vector, type = 'area')

		print "Vector conversion done"

		#Waterrise as a vector output

		selected_ocean = 'ocean_select_' + timestr + '_' + str(actual_loop_proj)

		gscript.run_command('v.select', ainput=ocean_vector, binput=ocean_point, output=selected_ocean, operator='intersects')

		selected_ocean_rast = 'ocean_select_rast_' + timestr + '_' + str(actual_loop_proj)

		gscript.run_command('v.to.rast', input=selected_ocean, output=selected_ocean_rast, use='val', value='0')

		#create temporary report about the area

		os.remove('temp.txt')

		gscript.run_command('r.report', map=selected_ocean_rast, units='k', output='temp.txt')

		tempTxt = open('temp.txt', 'r')

		for row in tempTxt:
			if row[1]=='0':
				AreaOut=row.split()
				Areatext=AreaOut[30].split('|')
				areaList.append(Areatext[0])

		tempTxt.close()

		print "select done"

		actual_loop_proj = actual_loop_proj + 1

	print "Here we are"

	#area_analysis()

def watershed():

	filled = 'DEM_filled' + timestr
	Direction = 'direction' + timestr
	Accumulation = 'accumulation' + timestr
	label_stat = 'label_stat' + timestr
	watershed = 'watershed' + timestr
	half_basin = 'half_basin' + timestr
	drainage = 'drainage' + timestr
	zonal_stat = 'zonal_stat' + timestr
	prepour_point = 'prepour_point' + timestr
	pour_point = 'pour_point' + timestr
	pour_point_vect=  'pour_point_vect' + timestr
	basin_perimeter = 'basin_perimeter' + timestr
	stream = 'stream' + timestr
	focal = 'focal' + timestr
	Accumulation_int = 'accumulation_int' + timestr
	stat_int = 'stat_int' + timestr

	gscript.run_command('g.region', rast=original)

	gscript.run_command('r.fill.dir', input = original, output = filled, direction = Direction)

	#gscript.run_command('r.terraflow', elev = original, filled = filled, direction = Direction, swatershed = s_watershed, accumulation=tAccumulation, tci = tci)

	gscript.run_command('r.watershed', flags = '4s', elevation = filled, accumulation = Accumulation)

	accMax = gscript.raster_info(Accumulation)['max']

	thresHold = accMax * 0.01

	print thresHold

	print accMax

	gscript.run_command('r.watershed', flags = '4s', elevation = filled, drainage = drainage, half_basin = half_basin, threshold = thresHold, stream = stream, basin = watershed)

	gscript.run_command('r.neighbors', input = watershed, output = focal, method = 'diversity')

	gscript.run_command('r.mapcalc', expression='%s = if(%s >= 2 , %s, null())' % (basin_perimeter, focal, watershed))

	gscript.run_command('r.mapcalc', expression='%s = int(%s)' % (Accumulation_int, Accumulation))

	gscript.run_command('r.statistics', base = basin_perimeter, cover = Accumulation_int, method = 'max', output = label_stat)

	gscript.run_command('r.mapcalc', expression='%s=@%s' % (zonal_stat, label_stat))

	gscript.run_command('r.mapcalc', expression='%s = int(%s)' % (stat_int, zonal_stat))

	gscript.run_command('r.mapcalc', expression='%s = if(%s == %s, 10, null())' % (prepour_point, Accumulation_int, stat_int))

	gscript.run_command('r.mapcalc', expression='%s = if(%s > 0, %s, null())' % (pour_point, original, prepour_point))

	gscript.run_command('r.to.vect', input = pour_point, output = pour_point_vect, feature = 'point')

#second try
	#gscript.run_command('r.neighbors', input = watershed, output = focal, method = 'diversity')

	#gscript.run_command('r.mapcalc', expression='%s = if(%s >= 2 && %s >= 0, 10, null())' % (prepour_point, focal, stream))

	#gscript.run_command('r.mapcalc', expression='%s = if(%s > 0, %s, null())' % (pour_point, original, prepour_point))

	#gscript.run_command('r.to.vect', input = pour_point, output = pour_point_vect, feature = 'point')

#First try
	#gscript.run_command('r.stats.zonal', base = watershed, cover = Accumulation, method = 'max', output = statoutput)

	#gscript.run_command('r.mapcalc', expression='%s = int(%s)' % (flow_int, Accumulation))

	#gscript.run_command('r.mapcalc', expression='%s = int(%s)' % (stat_int, statoutput))

	#gscript.run_command('r.mapcalc', expression='%s = if(%s == %s, 10, null())' % (pour_point, flow_int, stat_int))



def create_point(x_ocean, y_ocean):

	try:
		point = ogr.Geometry(ogr.wkbPoint)
		point.AddPoint(x_ocean, y_ocean)

		print '%d, %d' % (point.GetX(), point.GetY())

		#geojson = point.ExportToJson()
		#print geojson

		# Create the output Driver
		outDriver = ogr.GetDriverByName('GeoJSON')

		# Create the output GeoJSON
		outDataSource = outDriver.CreateDataSource('test.geojson')
		outLayer = outDataSource.CreateLayer('test.geojson', geom_type=ogr.wkbPoint )

		# Get the output Layer's Feature Definition
		featureDefn = outLayer.GetLayerDefn()

		# create a new feature
		outFeature = ogr.Feature(featureDefn)

		# Set new geometry
		outFeature.SetGeometry(point)

		# Add new feature to output Layer
		outLayer.CreateFeature(outFeature)
	except:
		"Cannot create point"

def cleanup_data():
	try:
		gscript.run_command('g.remove', flags='f', type = 'raster', pattern='*')

		gscript.run_command('g.remove', flags='f', type = 'vector', pattern='*')

		print "Cleanup done"

	except:
		print "Cleanup could not be performed"

if __name__ == "__main__":
	main()