#!/usr/bin/env python
#****************************************************************
#*
#* MODULE:     v.quota, v 1.0.0
#*
#* AUTHOR(S):  Antonio Alliegro
#*
#* PURPOSE:    Calcola la quota di un punto di un vettoriale 3D
#*
#* COPYRIGHT:  (C) 2011 Antonio Alliegro Civil Engineer
#*             Salerno, Italy
#*             antonioall(at)libero.it
#*
#*             First Version: 2011/08/16
#*             Last  Version: 2011/08/16
#*
#*             This program is free software under the
#*             GNU General Public License (>=v2).
#*             Read the file COPYING that comes with GRASS
#*             for details.
#*
#****************************************************************

#%module
#% description: Calcola la quota di un punto di un vettoriale 3D
#% keywords: vector
#%end

#%option
#% key: map
#% type: string
#% gisprompt: old,vector,vector
#% key_desc: nome
#% description: Nome copertura vettoriale 3D
#% required: yes
#%end

#%option
#% key: x
#% type: double
#% key_desc: Coordinata X
#% description: Coordinata X
#% required: yes
#% answer: 2515000
#%end

#%option
#% key: y
#% type: double
#% key_desc: Coordinata Y
#% description: Coordinata Y
#% required: yes
#% answer: 4456000
#%end

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

def main():
    input = options['map']
    x = options['x']
    y = options['y']
    
    # Inializza la libreria GRASS
    G_gisinit('')
    
    # Controlla se il vettoriale e' nel percorso di ricerca
    mapset = G_find_vector2(input, "")
    if not mapset:
        grass.fatal("Il file vettoriale <%s> non e' stato trovato" % input)

    # Definisce la struttura della copertura
    map_info = pointer(Map_info())

    # Setta la copertura al livello 2 (level 2: topology)
    Vect_set_open_level(2)

    # Apre il file vettoriale
    Vect_open_old(map_info, input, mapset)

    # Verifica se il vettoriale e' 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)

    # Conversione del tipo float di Python nel tipo double di ctypes
    dx = c_double(float(x))
    dy = c_double(float(y))
    z = c_double()
    
    # CALCOLA LA QUOTA DEL PUNTO DI COORDINATE (x, y)
    # si noti che z viene passata per riferimento, difatti, la funzione
    # richiede che questo paramentro sia un puntatore a un double ctypes
    ret = Vect_tin_get_z(map_info, dx, dy, byref(z), None, None);
    if ret > 0:
        # Si osservi che per leggere il valore da una variable ctypes di tipo
        # double si utilizza la proprieta' .value
        grass.message("Coordinata X del punto:%15f" % dx.value)
        grass.message("Coordinata Y del punto:%15f" % dy.value)
        grass.message("Quota del punto (Z):%18f" % z.value)
    else:
        grass.fatal("Il punto non ricade nel vettoriale <%s>" % input)
        
    # Chiude la copertura
    Vect_close(map_info)
    
    return 0
    #End main

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