#!/usr/bin/env python
#****************************************************************
#*
#* MODULE:     tin.lago, v 1.0.0
#*
#* AUTHOR(S):  Antonio Alliegro
#*
#* PURPOSE:    Calcolo profondita' media e massima di un lago
#*             mediante campionamento
#*
#* COPYRIGHT:  (C) 2012 Antonio Alliegro Civil Engineer
#*             Salerno, Italy
#*             antonioall(at)libero.it
#*
#*             First Version: 2012/03/12
#*             Last  Version: 2012/03/12
#*
#*             This program is free software under the
#*             GNU General Public License (>=v2).
#*             Read the file COPYING that comes with GRASS
#*             for details.
#*
#****************************************************************

#%module
#% description: Calcolo profondita' media e massima di un lago mediante campionamento
#% keywords: TIN, vector
#%end

#%option
#% key: map
#% type: string
#% gisprompt: old,vector,vector
#% key_desc: tin
#% description: Nome copertura TIN
#% required: yes
#%end

#%option
#% key: lago
#% type: string
#% gisprompt: old,vector,vector
#% key_desc: vector
#% description: Nome copertura contorno lago
#% required: yes
#%end

#%option
#% key: quota
#% type: double
#% answer: 400.0
#% key_desc: float
#% description: Quota superfice lago (m s.m.)
#% required: yes
#%end

#%option
#% key: delta
#% type: double
#% answer: 1.0
#% key_desc: float
#% description: Delta campionamento
#% required: yes
#%end

import os, sys
import grass.script as grass
from grass.lib.gis    import *
from grass.lib.vector import *

def main():
    input = options['map']
    basin = options['lago']
    quota = float(options['quota'])
    delta = float(options['delta'])
    hmin = 99999.9
    hmed = 0.0
    n = 0
    
    # Inializza la libreria GRASS
    G_gisinit('')
    
    # Controlla se le coperture sono nel percorso di ricerca
    mapset = G_find_vector2(input, "")
    if not mapset:
        grass.fatal("Il file vettoriale <%s> non e' stato trovato" % input)
    mapset_b = G_find_vector2(basin, "")
    if not mapset_b:
        grass.fatal("Il file vettoriale <%s> non e' stato trovato" % basin)

    # Definisce la struttura delle coperture
    map_info = pointer(Map_info())
    map_info_b = pointer(Map_info())

    # Setta le coperture al livello 2 (level 2: topology)
    Vect_set_open_level(2)

    # Apre i files vettoriali
    Vect_open_old(map_info, input, mapset)
    Vect_open_old(map_info_b, basin, mapset_b)

    # Verifica se il TIN e' un vettoriale 3D
    if Vect_is_3d(map_info):
        grass.message("Il file vettoriale <%s> e' 3D" % input)
    else:
        grass.fatal("Il file vettoriale <%s> non e' 3D" % input)
    
    # Verifica se il vettoriale del bacino ha tipologia areale
    nareas = Vect_get_num_areas(map_info_b)
    if nareas == 0:
        grass.fatal("Il file vettoriale <%s> non ha tipologia areale" % basin)
    
    # Estrazione valori della regione corrente
    window = pointer(Cell_head())
    G_get_window(window)
    nrows = window.contents.rows
    ncols = window.contents.cols
    xref = window.contents.west
    yref = window.contents.south
    xres = window.contents.ew_res
    yres = window.contents.ns_res
    
    # Calcolo numero righe e colonne nuova griglia
    nx = int(nrows * xres / delta) + 1
    ny = int(ncols * yres / delta) + 1
    
    # Blocco elaborazione
    grass.message("Passo 1/2: Pre-elaborazione TIN...")
    z = c_double()
    G_percent(0, nx, 2)
    Vect_tin_get_z(map_info, xref, yref, byref(z), None, None)
    grass.message("...Fatto")
    
    grass.message("Passo 2/2: Calcolo quote...")
    for i in range(nx):
        for j in range(ny):
            x = xref + j * delta
            y = yref + i * delta
            if Vect_point_in_area(map_info_b, 1, x, y) == 1:
                Vect_tin_get_z(map_info, x, y, byref(z), None, None)
                if z.value <= quota:
                    hmed += z.value
                    if z.value < hmin:
                        hmin = z.value
                    n += 1
        G_percent(i, nx, 2)
    grass.message("...Fatto")
    
    # Risultati
    if n > 0:
        hmed = hmed / n
        grass.message("----------------------------------")
        grass.message("Delta campionamento = %12.3f" % delta)
        grass.message("Numero celle totali = %12d" % (nx * ny))
        grass.message("Numero celle lago = %14d" % n)
        grass.message("Profondita' media =%15.3f" % (quota - hmed))
        grass.message("Profondita' massima =%13.3f" % (quota - hmin))
        grass.message("Volume acqua lago =%15.3f" % ((quota - hmed) * n * delta * delta))
    
    # Chiusura coperture
    Vect_close(map_info)
    Vect_close(map_info_b)
    
    return 0
    #End main

if __name__ == "__main__":
    options, flags = grass.parser()
    sys.exit(main())