#!/bin/sh
#############################################################################
#
# reference.sh
#
# SCRIPT:   	create a reference vector map for raster maps
# AUTHOR(S):	Paolo Zatelli - 2006/05/14
# PURPOSE:  	create a vector file containing the boundaries of a set
#		of raster maps and build an associated table with
#		map names and extents. Optionally save it to a
#		ESRI shape file.
#   	    	
# COPYRIGHT:	2006 Paolo Zatelli
# LICENSE:	GPL 2.0 or (at your option) any later version
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
#
#############################################################################
# IMPORTANT: the script must be run from the GRASS shell in the directory 
# containing the raster file.

# pattern of the file to be scanned
pattern='*asc'

# Type of vector primitives to be created: line or area.
# Uncomment below the feature type to be written into the output map.
#primitive=line
primitive=area

# The name of the output map. Beware of the DBF file name limitations.
output_map=quadro_dtm

# Save to an ESRI Shapefile with the nema of the output map, in the current
# directory. 'yes' or 'no'
save_shape=yes

#############################################################################
# nothing to modify below this line

START_DATE=`date`
num=0

echo $pattern

# remove the output map if exists, GRASS vector modules do not overwrite
g.remove vect=$output_map

for namein in $pattern
	do

let "num=$num+1"

# strips directory and extension of the raster map file
col1=`expr index "$namein" /`
col2=`expr length "$namein"`

let "col2=$col2-4"
let "lungh=$col2-$col1"

# strips the file name from the complete path+filename string
nameout=${namein:col1:lungh}

# remove the map in GRASS if exists, r.in.gdal does not overwrite maps
g.remove rast=$nameout
#  import
r.in.gdal input=$namein output=$nameout -e -o

# set the region to the map extent
g.region rast=$nameout

# remove the GRASS map
g.remove rast=$nameout

# create a vector area or line (depending on the $parameter value) primitive
g.remove vect=temp_ref
v.in.region output=temp_ref type=$primitive

# set category to current map number
g.remove vect=temp_1_ref
v.category in=temp_ref out=temp_1_ref option=del
g.remove vect=temp_ref
v.category in=temp_1_ref out=temp_ref option=add cat=$num type=centroid

if (( $num>1 ))
then

# create the DB table
v.db.addtable temp_ref columns="map varchar(20), north double, south double, west double, east double"
# fill the first table row
echo "UPDATE temp_ref SET map=$nameout WHERE cat=$num" | db.execute

g.remove vect=temp_tot
g.copy vect=$output_map,temp_tot

else
# first iteration
g.copy vect=temp_ref,$output_map

# create the DB table
v.db.addtable $output_map columns="map varchar(20), north double, south double, west double, east double"

# fill the first table row
echo "UPDATE $output_map SET map=$nameout WHERE cat=$num" | db.execute

fi

# boundaries to the DB table
nlimit=0
for limit in north south west east
	do

	let "nlimit=$nlimit+1"	
	
	extents=`g.region -g|head -$nlimit|tail -1`
	lenght=`expr length "$extents"`
	extent=${extents:2:lenght}

	if (( $num == 1 ))
	then
		echo "UPDATE $output_map SET $limit=$extent WHERE cat=$num" | db.execute
	else
		echo "UPDATE temp_ref SET $limit=$extent WHERE cat=$num" | db.execute
	fi
done

if (( $num>1 ))
then
# patch maps
v.patch input=temp_tot,temp_ref output=$output_map  -e
fi

done 

# cleanup
g.remove vect=temp_1_ref
g.remove vect=temp_ref
g.remove vect=temp_tot

# show the result in monitor x0
g.region vect=$output_map
d.mon start=x0
d.vect map=$output_map color=255:0:0 lcolor=0:0:255 fcolor=170:170:170 display=shape,attr type=point,line,boundary,centroid,area icon=basic/x size=5 layer=1 att=map lsize=8 xref=left yref=center llayer=1

# show the table
echo "select * from $output_map"|db.select

if [ $save_shape == "yes" ]
then
# output to a shapefile in the current directory
v.out.ogr input=$output_map type=$primitive dsn=$PWD olayer=$output_map layer=1 format=ESRI_Shapefile
fi

STOP_DATE=`date`
echo $num " files processed"
echo "Start: "$START_DATE
echo "Stop: "$STOP_DATE