[GRASS-SVN] r30893 - grass-addons/raster/r.out.gmap/r.out.gmap
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Apr 7 10:03:40 EDT 2008
Author: neteler
Date: 2008-04-07 10:03:40 -0400 (Mon, 07 Apr 2008)
New Revision: 30893
Added:
grass-addons/raster/r.out.gmap/r.out.gmap/main_tmp
Log:
moved into right directory
Copied: grass-addons/raster/r.out.gmap/r.out.gmap/main_tmp (from rev 30892, grass-addons/raster/r.out.gmap/main_tmp)
===================================================================
--- grass-addons/raster/r.out.gmap/r.out.gmap/main_tmp (rev 0)
+++ grass-addons/raster/r.out.gmap/r.out.gmap/main_tmp 2008-04-07 14:03:40 UTC (rev 30893)
@@ -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&v=2&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, "© 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