[GRASS-SVN] r30894 - grass-addons/raster/r.out.gmap

svn_grass at osgeo.org svn_grass at osgeo.org
Mon Apr 7 10:04:53 EDT 2008


Author: neteler
Date: 2008-04-07 10:04:53 -0400 (Mon, 07 Apr 2008)
New Revision: 30894

Added:
   grass-addons/raster/r.out.gmap/main_tmp
Removed:
   grass-addons/raster/r.out.gmap/main_tmp
   grass-addons/raster/r.out.gmap/r.out.gmap/
Log:
moved into right directory

Deleted: grass-addons/raster/r.out.gmap/main_tmp
===================================================================
--- grass-addons/raster/r.out.gmap/main_tmp	2008-04-07 14:03:40 UTC (rev 30893)
+++ grass-addons/raster/r.out.gmap/main_tmp	2008-04-07 14:04:53 UTC (rev 30894)
@@ -1,546 +0,0 @@
-#!/bin/sh
-
-############################################################################
-#
-# MODULE:	r.out.gmap
-# AUTHOR(S):	Tomas Cebecauer.
-# PURPOSE:	Raster data tiling for Google Maps and Microsoft Virtual Earth
-# COPYRIGHT:	(C) 2008 by the GRASS Development Team
-#
-#		This program is free software under the GNU General Public
-#		License (>=v2). Read the file COPYING that comes with GRASS
-#		for details.
-# REQUIRES:     - bc, PROJ4 (cs2cs), sed, grep, awk, wc, cat
-#
-#############################################################################
-
-#%Module
-#% description: Exports raster map in the Google Map or Microsoft Virtual Earth format
-#% keywords: raster, export
-#%End
-#%flag
-#% key: e
-#% description: Microsoft Virtual Earth naming style (default is Google Map naming style)
-#%END
-#%flag
-#% key: f
-#% description: Force processing if the current location projection test fails (not recommended)
-#%END
-#%flag
-#% key: h
-#% description: Create an example html file with google maps API for exported data (google maps only)
-#%END
-#%flag
-#% key: m
-#% description: Cut the output to specified edges
-#%END
-#%option
-#% key: input
-#% type: string
-#% gisprompt: old,cell,raster
-#% description: Input raster map
-#% required : yes
-#%END
-#%option
-#% key: location
-#% type: string
-#% description: Location of input map
-#% required : no
-#%END
-#%option
-#% key: mapset
-#% type: string
-#% description: Mapset of input map
-#% required : no
-#%END
-#%option
-#% key: interpolation
-#% type: string
-#% description: Interpolation method to use for re-projection
-#% answer: nearest
-#% options: nearest,bilinear,cubic
-#% required : no
-#% multiple: no
-#%end
-#%option
-#% key: zoom
-#% type: double
-#% description: Google Map zoom level to export
-#% options : 0-21
-#% required : yes
-#%END
-#%option
-#% key: n
-#% type: double
-#% description: Custom northern edge in lat/long on wgs84 ellipsoid
-#% options : -85.05-85.05
-#% required : no
-#%END
-#%option
-#% key: s
-#% type: double
-#% description: Custom southern edge in lat/long on wgs84 ellipsoid
-#% options : -85.05-85.05
-#% required : no
-#%END
-#%option
-#% key: e
-#% type: double
-#% description: Custom eastern edge in lat/long on wgs84 ellipsoid
-#% options : -180-180
-#% required : no
-#%END
-#%option
-#% key: w
-#% type: double
-#% description: Custom western edge in lat/long on wgs84 ellipsoid
-#% options : -180-180
-#% required : no
-#%END
-#%option
-#% key: outdir
-#% type: string
-#% description: directory to write output tiles
-#% answer: googlemap
-#% required : no
-#%END
-
-#constants:
-M2PX=$(echo "5340353.715440871795 / 6378137" | bc -l) #correction from meters to pixels(zoom 17)
-PX2M=$(echo "6378137 / 5340353.715440871795" | bc -l) #correction from meters to pixels(zoom 17)
-LEFT=$(echo "-16777216 * $PX2M" | bc -l)  #minimum left coordinate of projection
-
-cleanup()
-{
-#remove data
-   g.remove rast=out$$ --quiet
-}
-
-exitprocedure()
-{
- g.message -e "User break!"
- cleanup
- exit 1
-}
-trap "exitprocedure" 2 3 15
-
-if  [ -z $GISBASE ] ; then
-    echo "You must be in GRASS GIS to run this program." 1>&2
- exit 1
-fi
-
-
-
-# check if bc exists
-if [ ! -x "`which bc`" ] ; then
-   g.message -e "'bc' program not found, install it first"
-   exit 1
-fi
-# check if cs2cs exists
-if [ ! -x "`which cs2cs`" ] ; then
-   g.message -e "'cs2cs' program not found, install it first\
-      http://proj.maptools.org"
-   exit 1
-fi
-# check if grep exists
-if [ ! -x "`which grep`" ] ; then
-   g.message -e "'grep' program not found, install it first"
-   exit 1
-fi
-# check if awk exists
-if [ ! -x "`which awk`" ] ; then
-   g.message -e "'awk' program not found, install it first"
-   exit 1
-fi
-# check if sed exists
-if [ ! -x "`which sed`" ] ; then
-   g.message -e "'sed' program not found, install it first"
-   exit 1
-fi
-# check if tr exists
-if [ ! -x "`which tr`" ] ; then
-   g.message -e "'tr' program not found, install it first"
-   exit 1
-fi
-# check if wc exists
-if [ ! -x "`which wc`" ] ; then
-   g.message -e "'wc' program not found, install it first"
-   exit 1
-fi
-# check if cat exists
-if [ ! -x "`which cat`" ] ; then
-   g.message -e "'cat' program not found, install it first"
-   exit 1
-fi
-
-# save command line
-if [ "$1" != "@ARGS_PARSED@" ] ; then
-    exec g.parser "$0" "$@"
-fi
-
-map="$GIS_OPT_INPUT"
-location="$GIS_OPT_LOCATION"
-mapset="$GIS_OPT_MAPSET"
-level="$GIS_OPT_ZOOM"
-interpol="$GIS_OPT_INTERPOLATION"
-newdir="$GIS_OPT_OUTDIR"
-
-# test if input map is set in form map at mapset 
-hasat=`echo $map | grep '@' |wc -c`
-if [ $hasat -ge "3" ]; then #at least a at b
-	mapset=`echo $map | awk -F'@' '{print $2}'`
-	map=`echo $map | awk -F'@' '{print $1}'`
-fi
-
-# test zoom level for MVE export
-if [  $GIS_FLAG_E -eq "1" -a $level -eq "0" ]; then
-    g.message -e "zoom level $level not allowed for Microsoft Virtual Earth export"
-	exit 1
-fi
-
-# take env variables from current GRASS session
-eval `g.gisenv`
-: ${GISBASE?} ${GISDBASE?} ${LOCATION_NAME?} ${MAPSET?}
-
-
-	GMlocation=$LOCATION_NAME
-	GMmapset=$MAPSET
-	GMgisdbase=$GISDBASE
-	GMmapsetpath=$GISDBASE"/"$LOCATION_NAME"/"$MAPSET
-
-# if location and/or mapset not set use current location/mapset
- 	if [ -z $mapset ]; then mapset=$GMmapset; fi
- 	if [ -z $location ]; then location=$GMlocation; fi
-
-	mapsetpath=$GISDBASE"/"$location"/"$mapset
-
-# check wheter input map mapset differs from current mapset
-	diffinputmapset=1
-	if [ $GMmapset = $mapset -a  $GMlocation = $location ]; then diffinputmapset=0;  fi
-	if [ $GMmapset != $mapset -a  $GMlocation = $location ]; then 
-		g.mapsets addmapset=$mapset
-		diffinputmapset=0;
-	fi
-
-
-
-#test whether the script is run from mapset with google maps projection 
-	gmproj='+proj=merc +lat_ts=0.0000000000 +lon_0=0.0000000000 +k_0=1.0000000000 +x_0=0 +y_0=0 +no_defs +a=6378137 +b=6378137'
-	latlonproj='+proj=longlat +a=6378137 +rf=298.257223563 +no_defs'
-
-correctproj=1
-	# test presence of all proj parameters
-reqired_par="proj lat_ts lon_0 x_0 y_0 a b no_defs"
-for reqpar in $reqired_par
-do
-	par=`g.proj -j | grep $reqpar`
-	if [ -z "$par" ]; then 
-		correctproj=0
-	fi
-done
-
-if [ $correctproj -eq 1 ]; then
-	# get individual proj parameters
-	eval `g.proj -j | tr '+' ' ' | grep -v no_defs`
-	: ${proj?} ${lat_ts?} ${lon_0?} ${x_0?} ${y_0?} ${a?} ${b?}  >&1
-	nodefs=`g.proj -j | grep no_defs`
-	
-		# test presence/value of text parameters
-	if [ -z $nodefs -o $proj != 'merc' ]; then correctproj=0;  fi
-	
-		# test values of numeric parameters
-	testline="$lat_ts!=0.0 || $lon_0!=0.0 || $x_0!=0.0 || $y_0!=0.0 || $a!=6378137 || $b!=6378137"
-	if [ $(echo $testline | bc -l) -eq 1 ]; then correctproj=0; fi
-fi	
-
-	# inform user
-if [  $correctproj -ne 1 ]; then
-	if [  $GIS_FLAG_F -eq "1" ]; then
-		g.message -w "projection for location \"$GMlocation\" differs from required google map mercator projection"
-	else
-		g.message -e "projection for location \"$GMlocation\" differs from required google map mercator projection"
-		currproj=`g.proj -jf`
-		g.message message="current projection: $currproj"
-		g.message message="required projection: $gmproj"
-		exit 1
-	fi
-fi
-
-#test whether the input grass mapset is used/exists
-if [ $diffinputmapset -eq 1 ]; then
-	msg=`g.mapset mapset=$mapset location=$location gisdbase=$GISDBASE 2>&1`
-	msgERR=`echo $msg | grep 'ERROR' | wc -c`
-	if [ $msgERR -gt 0 ] 
-	then
-		g.message -e "input Mapset does not exist, or is used by GRASS"
-		exit 1
-	else # location/mapset successfuly changed, return back to the google map location/mapset
-		g.mapset mapset=$GMmapset location=$GMlocation gisdbase=$GMgisdbase >& /dev/null
-	fi
-fi
-
-#creates new directory for output
- 	if [ -z $newdir ]; then
-		newdir="."
-	elif [ ! -d $newdir ]; then
-		mkdir $newdir
-	elif [ $newdir = "googlemap" ]; then
-		g.message "exporting tiles to ./googlemap"
-	fi 
-
-# check existence of input map
-	infile=$mapsetpath"/cell/"$map
-	if [ ! -f "$infile" ]; then
-		g.message -e "cannot find input map " $map
-		exit 1
-	fi
-
-# force creation of color file if not exists
-	colorfile=$mapsetpath"/colr/"$map
-	if [ ! -f "$colorfile" ]; then
-		if [ $diffinputmapset -eq 1 ]; then
-			g.mapset mapset=$mapset location=$location gisdbase=$GISDBASE >& /dev/null
-			r.colors $map rast=$map
-			g.mapset mapset=$GMmapset location=$GMlocation gisdbase=$GMgisdbase >& /dev/null
-		else
-			r.colors $map rast=$map
-		fi
-		g.message -w 'forced creation of colr file...'
-	fi
-
-
-# calculate resolution
-	res=$(echo "(2^(17-$level)) * $PX2M" | bc -l)
-	reshalf=$(echo "$res/2" | bc -l)
-	segsize=$(echo "$res*256" | bc -l)
-	g.region res=$res
-
-
-# get N,E,W,S and proj in input LOCATION
-	if [ $diffinputmapset -eq 1 ]; then
-		g.mapset mapset=$mapset location=$location gisdbase=$GISDBASE >& /dev/null 
-		inproj=`g.proj -jf`
-		eval `r.info -g $map`
-		: ${north?} ${south?} ${west?} ${east?}
-		# go back to the GM loacation
-		g.mapset mapset=$GMmapset location=$GMlocation gisdbase=$GMgisdbase >& /dev/null
-	else
-		inproj=`g.proj -jf`
-		eval `r.info -g $map`
-		: ${north?} ${south?} ${west?} ${east?}
-		minx=$west
-		maxx=$east
-		miny=$south
-		maxy=$north
-	fi
-
-# calculate bounding box in longitude/latitude and check for the limits
-	cs2cs_cmd='cs2cs -f %.10f '$inproj' +to '$latlonproj
-
-	#get min longitude
-	if [ -n "$GIS_OPT_W" ]; then 
-		minlon="$GIS_OPT_W" 
-	else
-		minlon1=`echo $west $north 0 | $cs2cs_cmd | awk '{print $1}'`
-		minlon2=`echo $west $south 0 | $cs2cs_cmd | awk '{print $1}'`
-		minlon=$minlon2
-		if [ $(echo "$minlon1 < $minlon2" | bc -l) -eq 1 ] ; then minlon=$minlon1; fi
-		if [ $(echo "$minlon <= -180.0" | bc -l) -eq 1 ]; then minlon=-179.999999999999; fi
-	fi
-	#get max longitude
-	if [ -n "$GIS_OPT_E" ]; then 
-		maxlon="$GIS_OPT_E"
-	else
-		maxlon1=`echo $east $north 0 | $cs2cs_cmd | awk '{print $1}'`
-		maxlon2=`echo $east $south 0 | $cs2cs_cmd | awk '{print $1}'`
-		maxlon=$maxlon2
-		if [ $(echo "$maxlon1 > $maxlon2" | bc -l) -eq 1 ] ; then maxlon=$maxlon1; fi
-		if [ $(echo "$maxlon >= 180.0" | bc -l) -eq 1 ]; then maxlon=179.999999999999; fi
-	fi
-	#get min latitude
-	if [ -n "$GIS_OPT_S" ]; then 
-		minlat="$GIS_OPT_S"; 
-	else 
-		minlat1=`echo $west $south 0 | $cs2cs_cmd | awk '{print $2}'`
-		minlat2=`echo $east $south 0 | $cs2cs_cmd | awk '{print $2}'`
-		minlat=$minlat2
-		if [ $(echo "$minlat1 < $minlat2" | bc -l) -eq 1 ] ; then minlat=$minlat1; fi
-		if [ $(echo "$minlat < -85.05112877980659" | bc -l) -eq 1 ]; then minlat=-85.05112877980659; fi
-	fi
-	#get max latitude
-	if [ -n "$GIS_OPT_N" ]; then 
-		maxlat="$GIS_OPT_N"; 
-	else
-		maxlat1=`echo $west $north 0 | $cs2cs_cmd | awk '{print $2}'`
-		maxlat2=`echo $east $north 0 | $cs2cs_cmd | awk '{print $2}'`
-		maxlat=$maxlat2
-		if [ $(echo "$maxlat1 > $maxlat2" | bc -l) -eq 1 ] ; then maxlat=$maxlat1; fi
-		if [ $(echo "$maxlat > 85.05112877980659" | bc -l) -eq 1 ]; then maxlat=85.05112877980659; fi
-	fi
-
-
-	cs2cs_cmd='cs2cs -f %.10f '$latlonproj' +to '$gmproj
-
-	if [ $diffinputmapset -eq 1 -o -n "$GIS_OPT_W" ]; then
-		minx=`echo $minlon $minlat 0 | $cs2cs_cmd | awk '{print $1}'`
-	fi
-	if [ $diffinputmapset -eq 1 -o -n "$GIS_OPT_E" ]; then
-		maxx=`echo $maxlon $minlat 0 | $cs2cs_cmd | awk '{print $1}'`
-	fi
-	if [ $diffinputmapset -eq 1 -o -n "$GIS_OPT_S" ]; then
-		miny=`echo $minlon $minlat 0 | $cs2cs_cmd | awk '{print $2}'`
-	fi
-	if [ $diffinputmapset -eq 1 -o -n "$GIS_OPT_N" ]; then
-		maxy=`echo $minlon $maxlat 0 | $cs2cs_cmd | awk '{print $2}'`
-	fi
-
-# remove mask if exists
-	if [ -f $GMmapsetpath/cell/MASK ] ; then
-		g.message "removing mask"
-		g.rename rast=MASK,MASK$$ >& /dev/null
-	fi
-
-# tiles usually cover larger area than bounding box, to clip it to box uncomment this
-	if [  $GIS_FLAG_M -eq "1" ]; then 
-		g.region n=$maxy s=$miny e=$maxx w=$minx res=$res >& /dev/null
-		r.mapcalc MASK=1 >& /dev/null
-	fi
-
-#calcuate tiles
-	#reshalf is used to calculate centers of pixels from edges of pixels
-	colleft=$(echo "numb=(($minx - $LEFT + $reshalf) / $segsize);scale = 0;numb/1" | bc -l)
-	colright=$(echo "numb=(($maxx - $LEFT - $reshalf) / $segsize);scale = 0;numb/1" | bc -l)
-	rowtop=$(echo "numb=((0.0 - $maxy - $LEFT + $reshalf) / $segsize);scale = 0;numb/1" | bc -l)
-	rowbottom=$(echo "numb=((0.0 - $miny - $LEFT - $reshalf) / $segsize);scale = 0;numb/1" | bc -l)
-
-	g.message "Tiles: zoom $level cols $colleft - $colright rows $rowtop - $rowbottom"
-
-#finaly loop through the tiles and export 
-ncolorfile=$GMmapsetpath"/colr/out"$$
-
-i=$colleft
-while [ $i -le $colright ]
-do
-   leftt=$(echo "numb=($i * $segsize) + $LEFT;scale = 10;numb/1" | bc -l) 
-   right=$(echo "numb=($leftt + $segsize) ;scale = 10;numb/1" | bc -l)
-
-   j=$rowtop
-   while [ $j -le $rowbottom ]
-   do
-     top=$(echo "numb=0.0 - ($j * $segsize) - $LEFT;scale = 10;numb/1" | bc -l)
-     bottom=$(echo "numb=($top - $segsize) ;scale = 10;numb/1" | bc -l)
-
-	g.message "tile $i $j"
-
-	g.region w=$leftt e=$right s=$bottom n=$top res=$res
-
-#decide what will be exported
-	if [ $diffinputmapset -eq 1 ]; then
-		g.message -v "projecting..."
-		r.proj input=$map location=$location mapset=$mapset method=$interpol output=out$$ --quiet >& /dev/null
-	else	
-		r.mapcalc out$$=$map 
-	fi
-
-
-	#if nothing was projected (may happen if outside tha data) create empty file
-	nfile=$GMmapsetpath"/cell/out"$$
-	if [ ! -f $nfile ];	then
-		r.mapcalc out$$='null()' 2>&1 /dev/null
-	fi
-
-	cp $colorfile $ncolorfile    # just copy color file 
-
-
-	#export data
-	if [ $GIS_FLAG_E -eq "1" ] ; then
-			#file name in microsoft virtual earth style
-			mve_tile=$(echo "scale=0; res=0; lev=$level; anx=$i; any=$j; while(lev > 0){rx=0; ry=0; tileshalf=(2^(lev-1)); if (anx >= tileshalf) {anx=anx - tileshalf; rx=1};if (any >= tileshalf) {any=any - tileshalf; ry=1};res=res+((10^(lev-1))*((2*ry)+rx));lev=--lev }; res" | bc -l)
-			mve_tile=`printf "%0"$level"d" $mve_tile`
-			mve_file=$newdir"/"$mve_tile".png"
-			outfile=$mve_file
-	else
-			# file name in google map style
-			gm_file=$newdir"/z"$level"_"$i"x"$j".png"
-			outfile=$gm_file
-	fi	
-
-	g.message -v "exporting $outfile ..."
-
-	r.out.png --quiet input=out$$ output=$outfile
-
-#remove data
-	if [ -f $GMmapsetpath"/cell/out"$$ ]
-		then
-		g.remove out$$ >& /dev/null
-	fi
-
-     j=`expr $j + 1`
-   done  
-   i=`expr $i + 1`
-done
-
-
-#### html file
-if [  $GIS_FLAG_H -eq "1" -a $GIS_FLAG_E -ne "1"  ]; then 
-
-	midlon=$(echo "($maxlon + $minlon)/2" | bc -l)
-	midlat=$(echo "($maxlat + $minlat)/2" | bc -l)
-
-	cat << EOF > "${map}.html"
-<!DOCTYPE html "-//W3C//DTD XHTML 1.0 Strict//EN" 
-"http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd">
-<html xmlns="http://www.w3.org/1999/xhtml">
-<head>
-	<meta http-equiv="content-type" content="text/html; charset=utf-8"/>
-	<title>GRASS GIS r.out.gmap result</title>
-	<script src="http://maps.google.com/maps?file=api&amp;v=2&amp;key=abcdefg"
-			type="text/javascript"></script>
-	<script type="text/javascript">
-
-	function initialize() {
-	if (GBrowserIsCompatible()) {
-		var map = new GMap2(document.getElementById("map_canvas"));
-		/* add layer from grass*/
-		mapbounds = new GLatLngBounds(new GLatLng($minlat, $minlon), new GLatLng($maxlat, $maxlon));
-		var copyCollection = new GCopyrightCollection("");
-		copyCollection.addCopyright(new GCopyright(17, mapbounds, 2, "&#169 2008 GRASS GIS<br>"));
-
-		var anewLayer = new GTileLayer(copyCollection, 1,171);
-		anewLayer.getTileUrl = function(a,b){ return "${newdir}/z"+b+"_"+(a.x)+"x"+(a.y)+".png";};	
-
-		var anewMapType = new GMapType([anewLayer], G_NORMAL_MAP.getProjection(), "Grassmap");
-		anewMapType.getMinimumResolution = function() {return 0;};
-		anewMapType.getMaximumResolution = function() {return 17;};
-		map.addMapType(anewMapType);
-
-		map.addControl(new GLargeMapControl());
-		map.addControl(new GMapTypeControl());
-		map.setCenter(new GLatLng(${midlat}, ${midlon}), ${level}, anewMapType);
-	}
-	}
-
-	</script>
-</head>
-<body onload="initialize()" onunload="GUnload()">
-	<div id="map_canvas" style="width: 500px; height: 400px"></div>
-</body>
-</html>
-EOF
-
-	g.message "'${map}.html' created."
-fi
-
-#remove custom mask
-if [  $GIS_FLAG_M -eq "1" ]; then 
-	g.remove rast=MASK >& /dev/null
-fi
-
-#restore mask if exists
-if [ -f $GMmapsetpath/cell/MASK$$ ] ; then
-	g.message "restoring mask"
-	g.remove MASK  >& /dev/null
-	g.rename rast=MASK$$,MASK  >& /dev/null
-fi
-
-

Copied: grass-addons/raster/r.out.gmap/main_tmp (from rev 30893, grass-addons/raster/r.out.gmap/r.out.gmap/main_tmp)
===================================================================
--- grass-addons/raster/r.out.gmap/main_tmp	                        (rev 0)
+++ grass-addons/raster/r.out.gmap/main_tmp	2008-04-07 14:04:53 UTC (rev 30894)
@@ -0,0 +1,546 @@
+#!/bin/sh
+
+############################################################################
+#
+# MODULE:	r.out.gmap
+# AUTHOR(S):	Tomas Cebecauer.
+# PURPOSE:	Raster data tiling for Google Maps and Microsoft Virtual Earth
+# COPYRIGHT:	(C) 2008 by the GRASS Development Team
+#
+#		This program is free software under the GNU General Public
+#		License (>=v2). Read the file COPYING that comes with GRASS
+#		for details.
+# REQUIRES:     - bc, PROJ4 (cs2cs), sed, grep, awk, wc, cat
+#
+#############################################################################
+
+#%Module
+#% description: Exports raster map in the Google Map or Microsoft Virtual Earth format
+#% keywords: raster, export
+#%End
+#%flag
+#% key: e
+#% description: Microsoft Virtual Earth naming style (default is Google Map naming style)
+#%END
+#%flag
+#% key: f
+#% description: Force processing if the current location projection test fails (not recommended)
+#%END
+#%flag
+#% key: h
+#% description: Create an example html file with google maps API for exported data (google maps only)
+#%END
+#%flag
+#% key: m
+#% description: Cut the output to specified edges
+#%END
+#%option
+#% key: input
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Input raster map
+#% required : yes
+#%END
+#%option
+#% key: location
+#% type: string
+#% description: Location of input map
+#% required : no
+#%END
+#%option
+#% key: mapset
+#% type: string
+#% description: Mapset of input map
+#% required : no
+#%END
+#%option
+#% key: interpolation
+#% type: string
+#% description: Interpolation method to use for re-projection
+#% answer: nearest
+#% options: nearest,bilinear,cubic
+#% required : no
+#% multiple: no
+#%end
+#%option
+#% key: zoom
+#% type: double
+#% description: Google Map zoom level to export
+#% options : 0-21
+#% required : yes
+#%END
+#%option
+#% key: n
+#% type: double
+#% description: Custom northern edge in lat/long on wgs84 ellipsoid
+#% options : -85.05-85.05
+#% required : no
+#%END
+#%option
+#% key: s
+#% type: double
+#% description: Custom southern edge in lat/long on wgs84 ellipsoid
+#% options : -85.05-85.05
+#% required : no
+#%END
+#%option
+#% key: e
+#% type: double
+#% description: Custom eastern edge in lat/long on wgs84 ellipsoid
+#% options : -180-180
+#% required : no
+#%END
+#%option
+#% key: w
+#% type: double
+#% description: Custom western edge in lat/long on wgs84 ellipsoid
+#% options : -180-180
+#% required : no
+#%END
+#%option
+#% key: outdir
+#% type: string
+#% description: directory to write output tiles
+#% answer: googlemap
+#% required : no
+#%END
+
+#constants:
+M2PX=$(echo "5340353.715440871795 / 6378137" | bc -l) #correction from meters to pixels(zoom 17)
+PX2M=$(echo "6378137 / 5340353.715440871795" | bc -l) #correction from meters to pixels(zoom 17)
+LEFT=$(echo "-16777216 * $PX2M" | bc -l)  #minimum left coordinate of projection
+
+cleanup()
+{
+#remove data
+   g.remove rast=out$$ --quiet
+}
+
+exitprocedure()
+{
+ g.message -e "User break!"
+ cleanup
+ exit 1
+}
+trap "exitprocedure" 2 3 15
+
+if  [ -z $GISBASE ] ; then
+    echo "You must be in GRASS GIS to run this program." 1>&2
+ exit 1
+fi
+
+
+
+# check if bc exists
+if [ ! -x "`which bc`" ] ; then
+   g.message -e "'bc' program not found, install it first"
+   exit 1
+fi
+# check if cs2cs exists
+if [ ! -x "`which cs2cs`" ] ; then
+   g.message -e "'cs2cs' program not found, install it first\
+      http://proj.maptools.org"
+   exit 1
+fi
+# check if grep exists
+if [ ! -x "`which grep`" ] ; then
+   g.message -e "'grep' program not found, install it first"
+   exit 1
+fi
+# check if awk exists
+if [ ! -x "`which awk`" ] ; then
+   g.message -e "'awk' program not found, install it first"
+   exit 1
+fi
+# check if sed exists
+if [ ! -x "`which sed`" ] ; then
+   g.message -e "'sed' program not found, install it first"
+   exit 1
+fi
+# check if tr exists
+if [ ! -x "`which tr`" ] ; then
+   g.message -e "'tr' program not found, install it first"
+   exit 1
+fi
+# check if wc exists
+if [ ! -x "`which wc`" ] ; then
+   g.message -e "'wc' program not found, install it first"
+   exit 1
+fi
+# check if cat exists
+if [ ! -x "`which cat`" ] ; then
+   g.message -e "'cat' program not found, install it first"
+   exit 1
+fi
+
+# save command line
+if [ "$1" != "@ARGS_PARSED@" ] ; then
+    exec g.parser "$0" "$@"
+fi
+
+map="$GIS_OPT_INPUT"
+location="$GIS_OPT_LOCATION"
+mapset="$GIS_OPT_MAPSET"
+level="$GIS_OPT_ZOOM"
+interpol="$GIS_OPT_INTERPOLATION"
+newdir="$GIS_OPT_OUTDIR"
+
+# test if input map is set in form map at mapset 
+hasat=`echo $map | grep '@' |wc -c`
+if [ $hasat -ge "3" ]; then #at least a at b
+	mapset=`echo $map | awk -F'@' '{print $2}'`
+	map=`echo $map | awk -F'@' '{print $1}'`
+fi
+
+# test zoom level for MVE export
+if [  $GIS_FLAG_E -eq "1" -a $level -eq "0" ]; then
+    g.message -e "zoom level $level not allowed for Microsoft Virtual Earth export"
+	exit 1
+fi
+
+# take env variables from current GRASS session
+eval `g.gisenv`
+: ${GISBASE?} ${GISDBASE?} ${LOCATION_NAME?} ${MAPSET?}
+
+
+	GMlocation=$LOCATION_NAME
+	GMmapset=$MAPSET
+	GMgisdbase=$GISDBASE
+	GMmapsetpath=$GISDBASE"/"$LOCATION_NAME"/"$MAPSET
+
+# if location and/or mapset not set use current location/mapset
+ 	if [ -z $mapset ]; then mapset=$GMmapset; fi
+ 	if [ -z $location ]; then location=$GMlocation; fi
+
+	mapsetpath=$GISDBASE"/"$location"/"$mapset
+
+# check wheter input map mapset differs from current mapset
+	diffinputmapset=1
+	if [ $GMmapset = $mapset -a  $GMlocation = $location ]; then diffinputmapset=0;  fi
+	if [ $GMmapset != $mapset -a  $GMlocation = $location ]; then 
+		g.mapsets addmapset=$mapset
+		diffinputmapset=0;
+	fi
+
+
+
+#test whether the script is run from mapset with google maps projection 
+	gmproj='+proj=merc +lat_ts=0.0000000000 +lon_0=0.0000000000 +k_0=1.0000000000 +x_0=0 +y_0=0 +no_defs +a=6378137 +b=6378137'
+	latlonproj='+proj=longlat +a=6378137 +rf=298.257223563 +no_defs'
+
+correctproj=1
+	# test presence of all proj parameters
+reqired_par="proj lat_ts lon_0 x_0 y_0 a b no_defs"
+for reqpar in $reqired_par
+do
+	par=`g.proj -j | grep $reqpar`
+	if [ -z "$par" ]; then 
+		correctproj=0
+	fi
+done
+
+if [ $correctproj -eq 1 ]; then
+	# get individual proj parameters
+	eval `g.proj -j | tr '+' ' ' | grep -v no_defs`
+	: ${proj?} ${lat_ts?} ${lon_0?} ${x_0?} ${y_0?} ${a?} ${b?}  >&1
+	nodefs=`g.proj -j | grep no_defs`
+	
+		# test presence/value of text parameters
+	if [ -z $nodefs -o $proj != 'merc' ]; then correctproj=0;  fi
+	
+		# test values of numeric parameters
+	testline="$lat_ts!=0.0 || $lon_0!=0.0 || $x_0!=0.0 || $y_0!=0.0 || $a!=6378137 || $b!=6378137"
+	if [ $(echo $testline | bc -l) -eq 1 ]; then correctproj=0; fi
+fi	
+
+	# inform user
+if [  $correctproj -ne 1 ]; then
+	if [  $GIS_FLAG_F -eq "1" ]; then
+		g.message -w "projection for location \"$GMlocation\" differs from required google map mercator projection"
+	else
+		g.message -e "projection for location \"$GMlocation\" differs from required google map mercator projection"
+		currproj=`g.proj -jf`
+		g.message message="current projection: $currproj"
+		g.message message="required projection: $gmproj"
+		exit 1
+	fi
+fi
+
+#test whether the input grass mapset is used/exists
+if [ $diffinputmapset -eq 1 ]; then
+	msg=`g.mapset mapset=$mapset location=$location gisdbase=$GISDBASE 2>&1`
+	msgERR=`echo $msg | grep 'ERROR' | wc -c`
+	if [ $msgERR -gt 0 ] 
+	then
+		g.message -e "input Mapset does not exist, or is used by GRASS"
+		exit 1
+	else # location/mapset successfuly changed, return back to the google map location/mapset
+		g.mapset mapset=$GMmapset location=$GMlocation gisdbase=$GMgisdbase >& /dev/null
+	fi
+fi
+
+#creates new directory for output
+ 	if [ -z $newdir ]; then
+		newdir="."
+	elif [ ! -d $newdir ]; then
+		mkdir $newdir
+	elif [ $newdir = "googlemap" ]; then
+		g.message "exporting tiles to ./googlemap"
+	fi 
+
+# check existence of input map
+	infile=$mapsetpath"/cell/"$map
+	if [ ! -f "$infile" ]; then
+		g.message -e "cannot find input map " $map
+		exit 1
+	fi
+
+# force creation of color file if not exists
+	colorfile=$mapsetpath"/colr/"$map
+	if [ ! -f "$colorfile" ]; then
+		if [ $diffinputmapset -eq 1 ]; then
+			g.mapset mapset=$mapset location=$location gisdbase=$GISDBASE >& /dev/null
+			r.colors $map rast=$map
+			g.mapset mapset=$GMmapset location=$GMlocation gisdbase=$GMgisdbase >& /dev/null
+		else
+			r.colors $map rast=$map
+		fi
+		g.message -w 'forced creation of colr file...'
+	fi
+
+
+# calculate resolution
+	res=$(echo "(2^(17-$level)) * $PX2M" | bc -l)
+	reshalf=$(echo "$res/2" | bc -l)
+	segsize=$(echo "$res*256" | bc -l)
+	g.region res=$res
+
+
+# get N,E,W,S and proj in input LOCATION
+	if [ $diffinputmapset -eq 1 ]; then
+		g.mapset mapset=$mapset location=$location gisdbase=$GISDBASE >& /dev/null 
+		inproj=`g.proj -jf`
+		eval `r.info -g $map`
+		: ${north?} ${south?} ${west?} ${east?}
+		# go back to the GM loacation
+		g.mapset mapset=$GMmapset location=$GMlocation gisdbase=$GMgisdbase >& /dev/null
+	else
+		inproj=`g.proj -jf`
+		eval `r.info -g $map`
+		: ${north?} ${south?} ${west?} ${east?}
+		minx=$west
+		maxx=$east
+		miny=$south
+		maxy=$north
+	fi
+
+# calculate bounding box in longitude/latitude and check for the limits
+	cs2cs_cmd='cs2cs -f %.10f '$inproj' +to '$latlonproj
+
+	#get min longitude
+	if [ -n "$GIS_OPT_W" ]; then 
+		minlon="$GIS_OPT_W" 
+	else
+		minlon1=`echo $west $north 0 | $cs2cs_cmd | awk '{print $1}'`
+		minlon2=`echo $west $south 0 | $cs2cs_cmd | awk '{print $1}'`
+		minlon=$minlon2
+		if [ $(echo "$minlon1 < $minlon2" | bc -l) -eq 1 ] ; then minlon=$minlon1; fi
+		if [ $(echo "$minlon <= -180.0" | bc -l) -eq 1 ]; then minlon=-179.999999999999; fi
+	fi
+	#get max longitude
+	if [ -n "$GIS_OPT_E" ]; then 
+		maxlon="$GIS_OPT_E"
+	else
+		maxlon1=`echo $east $north 0 | $cs2cs_cmd | awk '{print $1}'`
+		maxlon2=`echo $east $south 0 | $cs2cs_cmd | awk '{print $1}'`
+		maxlon=$maxlon2
+		if [ $(echo "$maxlon1 > $maxlon2" | bc -l) -eq 1 ] ; then maxlon=$maxlon1; fi
+		if [ $(echo "$maxlon >= 180.0" | bc -l) -eq 1 ]; then maxlon=179.999999999999; fi
+	fi
+	#get min latitude
+	if [ -n "$GIS_OPT_S" ]; then 
+		minlat="$GIS_OPT_S"; 
+	else 
+		minlat1=`echo $west $south 0 | $cs2cs_cmd | awk '{print $2}'`
+		minlat2=`echo $east $south 0 | $cs2cs_cmd | awk '{print $2}'`
+		minlat=$minlat2
+		if [ $(echo "$minlat1 < $minlat2" | bc -l) -eq 1 ] ; then minlat=$minlat1; fi
+		if [ $(echo "$minlat < -85.05112877980659" | bc -l) -eq 1 ]; then minlat=-85.05112877980659; fi
+	fi
+	#get max latitude
+	if [ -n "$GIS_OPT_N" ]; then 
+		maxlat="$GIS_OPT_N"; 
+	else
+		maxlat1=`echo $west $north 0 | $cs2cs_cmd | awk '{print $2}'`
+		maxlat2=`echo $east $north 0 | $cs2cs_cmd | awk '{print $2}'`
+		maxlat=$maxlat2
+		if [ $(echo "$maxlat1 > $maxlat2" | bc -l) -eq 1 ] ; then maxlat=$maxlat1; fi
+		if [ $(echo "$maxlat > 85.05112877980659" | bc -l) -eq 1 ]; then maxlat=85.05112877980659; fi
+	fi
+
+
+	cs2cs_cmd='cs2cs -f %.10f '$latlonproj' +to '$gmproj
+
+	if [ $diffinputmapset -eq 1 -o -n "$GIS_OPT_W" ]; then
+		minx=`echo $minlon $minlat 0 | $cs2cs_cmd | awk '{print $1}'`
+	fi
+	if [ $diffinputmapset -eq 1 -o -n "$GIS_OPT_E" ]; then
+		maxx=`echo $maxlon $minlat 0 | $cs2cs_cmd | awk '{print $1}'`
+	fi
+	if [ $diffinputmapset -eq 1 -o -n "$GIS_OPT_S" ]; then
+		miny=`echo $minlon $minlat 0 | $cs2cs_cmd | awk '{print $2}'`
+	fi
+	if [ $diffinputmapset -eq 1 -o -n "$GIS_OPT_N" ]; then
+		maxy=`echo $minlon $maxlat 0 | $cs2cs_cmd | awk '{print $2}'`
+	fi
+
+# remove mask if exists
+	if [ -f $GMmapsetpath/cell/MASK ] ; then
+		g.message "removing mask"
+		g.rename rast=MASK,MASK$$ >& /dev/null
+	fi
+
+# tiles usually cover larger area than bounding box, to clip it to box uncomment this
+	if [  $GIS_FLAG_M -eq "1" ]; then 
+		g.region n=$maxy s=$miny e=$maxx w=$minx res=$res >& /dev/null
+		r.mapcalc MASK=1 >& /dev/null
+	fi
+
+#calcuate tiles
+	#reshalf is used to calculate centers of pixels from edges of pixels
+	colleft=$(echo "numb=(($minx - $LEFT + $reshalf) / $segsize);scale = 0;numb/1" | bc -l)
+	colright=$(echo "numb=(($maxx - $LEFT - $reshalf) / $segsize);scale = 0;numb/1" | bc -l)
+	rowtop=$(echo "numb=((0.0 - $maxy - $LEFT + $reshalf) / $segsize);scale = 0;numb/1" | bc -l)
+	rowbottom=$(echo "numb=((0.0 - $miny - $LEFT - $reshalf) / $segsize);scale = 0;numb/1" | bc -l)
+
+	g.message "Tiles: zoom $level cols $colleft - $colright rows $rowtop - $rowbottom"
+
+#finaly loop through the tiles and export 
+ncolorfile=$GMmapsetpath"/colr/out"$$
+
+i=$colleft
+while [ $i -le $colright ]
+do
+   leftt=$(echo "numb=($i * $segsize) + $LEFT;scale = 10;numb/1" | bc -l) 
+   right=$(echo "numb=($leftt + $segsize) ;scale = 10;numb/1" | bc -l)
+
+   j=$rowtop
+   while [ $j -le $rowbottom ]
+   do
+     top=$(echo "numb=0.0 - ($j * $segsize) - $LEFT;scale = 10;numb/1" | bc -l)
+     bottom=$(echo "numb=($top - $segsize) ;scale = 10;numb/1" | bc -l)
+
+	g.message "tile $i $j"
+
+	g.region w=$leftt e=$right s=$bottom n=$top res=$res
+
+#decide what will be exported
+	if [ $diffinputmapset -eq 1 ]; then
+		g.message -v "projecting..."
+		r.proj input=$map location=$location mapset=$mapset method=$interpol output=out$$ --quiet >& /dev/null
+	else	
+		r.mapcalc out$$=$map 
+	fi
+
+
+	#if nothing was projected (may happen if outside tha data) create empty file
+	nfile=$GMmapsetpath"/cell/out"$$
+	if [ ! -f $nfile ];	then
+		r.mapcalc out$$='null()' 2>&1 /dev/null
+	fi
+
+	cp $colorfile $ncolorfile    # just copy color file 
+
+
+	#export data
+	if [ $GIS_FLAG_E -eq "1" ] ; then
+			#file name in microsoft virtual earth style
+			mve_tile=$(echo "scale=0; res=0; lev=$level; anx=$i; any=$j; while(lev > 0){rx=0; ry=0; tileshalf=(2^(lev-1)); if (anx >= tileshalf) {anx=anx - tileshalf; rx=1};if (any >= tileshalf) {any=any - tileshalf; ry=1};res=res+((10^(lev-1))*((2*ry)+rx));lev=--lev }; res" | bc -l)
+			mve_tile=`printf "%0"$level"d" $mve_tile`
+			mve_file=$newdir"/"$mve_tile".png"
+			outfile=$mve_file
+	else
+			# file name in google map style
+			gm_file=$newdir"/z"$level"_"$i"x"$j".png"
+			outfile=$gm_file
+	fi	
+
+	g.message -v "exporting $outfile ..."
+
+	r.out.png --quiet input=out$$ output=$outfile
+
+#remove data
+	if [ -f $GMmapsetpath"/cell/out"$$ ]
+		then
+		g.remove out$$ >& /dev/null
+	fi
+
+     j=`expr $j + 1`
+   done  
+   i=`expr $i + 1`
+done
+
+
+#### html file
+if [  $GIS_FLAG_H -eq "1" -a $GIS_FLAG_E -ne "1"  ]; then 
+
+	midlon=$(echo "($maxlon + $minlon)/2" | bc -l)
+	midlat=$(echo "($maxlat + $minlat)/2" | bc -l)
+
+	cat << EOF > "${map}.html"
+<!DOCTYPE html "-//W3C//DTD XHTML 1.0 Strict//EN" 
+"http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd">
+<html xmlns="http://www.w3.org/1999/xhtml">
+<head>
+	<meta http-equiv="content-type" content="text/html; charset=utf-8"/>
+	<title>GRASS GIS r.out.gmap result</title>
+	<script src="http://maps.google.com/maps?file=api&amp;v=2&amp;key=abcdefg"
+			type="text/javascript"></script>
+	<script type="text/javascript">
+
+	function initialize() {
+	if (GBrowserIsCompatible()) {
+		var map = new GMap2(document.getElementById("map_canvas"));
+		/* add layer from grass*/
+		mapbounds = new GLatLngBounds(new GLatLng($minlat, $minlon), new GLatLng($maxlat, $maxlon));
+		var copyCollection = new GCopyrightCollection("");
+		copyCollection.addCopyright(new GCopyright(17, mapbounds, 2, "&#169 2008 GRASS GIS<br>"));
+
+		var anewLayer = new GTileLayer(copyCollection, 1,171);
+		anewLayer.getTileUrl = function(a,b){ return "${newdir}/z"+b+"_"+(a.x)+"x"+(a.y)+".png";};	
+
+		var anewMapType = new GMapType([anewLayer], G_NORMAL_MAP.getProjection(), "Grassmap");
+		anewMapType.getMinimumResolution = function() {return 0;};
+		anewMapType.getMaximumResolution = function() {return 17;};
+		map.addMapType(anewMapType);
+
+		map.addControl(new GLargeMapControl());
+		map.addControl(new GMapTypeControl());
+		map.setCenter(new GLatLng(${midlat}, ${midlon}), ${level}, anewMapType);
+	}
+	}
+
+	</script>
+</head>
+<body onload="initialize()" onunload="GUnload()">
+	<div id="map_canvas" style="width: 500px; height: 400px"></div>
+</body>
+</html>
+EOF
+
+	g.message "'${map}.html' created."
+fi
+
+#remove custom mask
+if [  $GIS_FLAG_M -eq "1" ]; then 
+	g.remove rast=MASK >& /dev/null
+fi
+
+#restore mask if exists
+if [ -f $GMmapsetpath/cell/MASK$$ ] ; then
+	g.message "restoring mask"
+	g.remove MASK  >& /dev/null
+	g.rename rast=MASK$$,MASK  >& /dev/null
+fi
+
+



More information about the grass-commit mailing list