[GRASS-SVN] r60364 - grass-addons/grass6/raster/r.out.mbtiles
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon May 19 22:06:00 PDT 2014
Author: hamish
Date: 2014-05-19 22:06:00 -0700 (Mon, 19 May 2014)
New Revision: 60364
Added:
grass-addons/grass6/raster/r.out.mbtiles/r.out.mbtiles
Removed:
grass-addons/grass6/raster/r.out.mbtiles/r.out.mbtiles_prep
Log:
rename module, part 2
Copied: grass-addons/grass6/raster/r.out.mbtiles/r.out.mbtiles (from rev 60363, grass-addons/grass6/raster/r.out.mbtiles/r.out.mbtiles_prep)
===================================================================
--- grass-addons/grass6/raster/r.out.mbtiles/r.out.mbtiles (rev 0)
+++ grass-addons/grass6/raster/r.out.mbtiles/r.out.mbtiles 2014-05-20 05:06:00 UTC (rev 60364)
@@ -0,0 +1,550 @@
+#!/bin/sh
+############################################################################
+#
+# MODULE: r.out.mbtiles
+# AUTHOR(S): M. Hamish Bowman, Dunedin, New Zealand
+# PURPOSE: Export a raster map into a TMS bundle ready for conversion
+# to MBTiles SQLite format
+# COPYRIGHT: (C) 2014 Hamish Bowman, and the GRASS Development Team
+#
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 2 of the License, 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.
+#
+# *** NOT FOR NAVIGATIONAL USE - USE AT YOUR OWN RISK ***
+#
+############################################################################
+
+# https://github.com/geopaparazzi/geopaparazzi/wiki/.mapurl-Converting-Tile-Set-to-.mbtiles
+# https://github.com/geopaparazzi/geopaparazzi/wiki/.mapurl-parameters
+# https://github.com/mapbox/mbtiles-spec/blob/master/1.1/spec.md
+# https://github.com/geopaparazzi/geopaparazzi/wiki/mbtiles-Implementation
+
+#%Module
+#% description: Export GRASS raster as a TMS tree ready for converting to mbtiles format.
+#% keywords: raster, export, tms, mbtiles
+#%End
+#%Option
+#% key: input
+#% type: string
+#% required: yes
+#% multiple: no
+#% key_desc: name
+#% description: Name of input raster map
+#% gisprompt: old,cell,raster
+#%End
+#%Option
+#% key: output
+#% type: string
+#% required: no
+#% multiple: no
+#% key_desc: name
+#% label: Name for output project directory
+#% description: default: input map name
+#% gisprompt: new_file,file,output
+#%End
+#%Option
+#% key: format
+#% type: string
+#% required: no
+#% multiple: no
+#% key_desc: string
+#% description: Image format to store the tiles
+#% options: png,jpeg
+#% answer: png
+#%End
+#%Option
+#% key: zoom
+#% type: string
+#% required: no
+#% multiple: no
+#% key_desc: string
+#% label: TMS zoom levels to create
+#% description: Format: '0-19', or '10' (default: automatic)
+#%End
+#%Option
+#% key: dpi
+#% type: integer
+#% label: DPI of the target display
+#% description: Value for a smartphone might be 225, a desktop perhaps 96.
+#% answer: 225
+#%End
+#%Flag
+#% key: m
+#% description: Create a MBTiles database (requires SQLite3)
+#%End
+#%Flag
+#% key: t
+#% description: Create a tarball of the tile tree
+#%End
+
+
+
+if [ -z "$GISBASE" ] ; then
+ echo "You must be in GRASS GIS to run this program." >&2
+ exit 1
+fi
+
+if [ "$1" != "@ARGS_PARSED@" ] ; then
+ exec g.parser "$0" "$@"
+fi
+
+#### check if we have awk, gdalwarp, gdal2tiles, seq
+pgms="awk seq gdalwarp gdalinfo gdal2tiles.py"
+for pgm in $pgms ; do
+ if [ ! -x "`which $pgm`" ] ; then
+ g.message -e "$pgm is required, please install it first."
+ exit 1
+ fi
+done
+
+
+MAP_NAME="$GIS_OPT_INPUT"
+
+if [ -n "$GIS_OPT_OUTPUT" ] ; then
+ OUTFILE="$GIS_OPT_OUTPUT"
+else
+ OUTFILE=`echo "$MAP_NAME" | cut -f1 -d'@'`
+fi
+
+
+g.findfile element=cell file="$MAP_NAME" > /dev/null
+if [ $? -eq 0 ] ; then
+ MAP_TYPE=raster
+else
+ g.findfile element=group file="$MAP_NAME" > /dev/null
+ if [ $? -eq 0 ] ; then
+ MAP_TYPE=group
+ else
+ g.message -e "Could not find <$MAP_NAME>"
+ exit 1
+ fi
+fi
+
+cleanup()
+{
+ #remove temporary files
+ if [ -e "$TMPFILE" ] ; then
+ rm -rf "$TMPFILE"*
+ fi
+}
+
+# what to do in case of user break:
+exitprocedure()
+{
+ g.message -e message='User break!'
+ cleanup
+ exit 1
+}
+
+# shell check for user break (signal list: trap -l)
+trap "exitprocedure" 2 3 15
+
+
+#### setup temporary file
+TMPFILE="`g.tempfile pid=$$`"
+if [ $? -ne 0 ] || [ -z "$TMPFILE" ] ; then
+ g.message -e "Unable to create temporary file"
+ exit 1
+fi
+
+
+TMP_GTIFF="${TMPFILE}_export.tif"
+TMP_DIR=`dirname "$TMP_GTIFF"`
+TMP_GMERC="$TMP_DIR"/`basename "$TMP_GTIFF" .tif`_gmerc.tif
+
+
+if [ -n "$GRASS_VERBOSE" ] && [ "$GRASS_VERBOSE" -eq 0 ] ; then
+ QUIET="-q"
+else
+ QUIET=""
+fi
+
+# ?? do the source of TMS tiles sometimes need to be reprojected into the
+# funny google merc projection first?
+
+if [ "$MAP_TYPE" = "raster" ] ; then
+ r.out.tiff -t input="$MAP_NAME" output="$TMP_GTIFF"
+ #gdalwarp -s_srs "`g.proj -jf`" -t_srs EPSG:900913 -dstalpha "$TMP_GTIFF" "$TMP_GMERC"
+ gdal_translate $QUIET -a_srs "`g.proj -jf`" "$TMP_GTIFF" "$TMP_GMERC"
+else
+ if [ `echo "$MAP_NAME" | grep -c '@'` -gt 0 ] ; then
+ g.message -e "Imagery groups must exist in the current mapset, sorry."
+ cleanup
+ exit 1
+ fi
+
+ r.out.gdal input="$MAP_NAME" output="$TMP_GTIFF" format=GTiff
+ #gdalwarp -t_srs EPSG:900913 -dstalpha "$TMP_GTIFF" "$TMP_GMERC"
+ mv "$TMP_GTIFF" "$TMP_GMERC"
+fi
+## gdalwarp -srcnodata -dstnodata -dstalpha ?
+
+if [ $? -ne 0 ] ; then
+ g.message -e "Problem exporting image"
+ cleanup
+ exit 1
+fi
+
+# free the disk space ASAP
+rm -f "$TMP_GTIFF"
+
+
+if [ 0 -eq 1 ] ; then
+
+## I don't think we need to do the gdalwarp, but it's worth keeping this
+## next bit of code for parsing the lat/lon bbox out of gdalinfo for a
+## rainy day.
+
+#### we need to find the bounding box in lat/lon of the post-warped geotiff
+# -- parse gdalinfo --
+gdalinfo -noct -nomd "$TMP_GMERC" > "$TMPFILE.warpinfo"
+
+IN_PROJ="+proj=longlat +datum=WGS84"
+OUT_PROJ="+proj=longlat +datum=WGS84"
+#?OUT_PROJ="+proj=longlat +ellps=sphere +a=6378137 +b=6378137 +nadgrids=@null"
+
+
+UL_DMS=`grep '^Upper Left ' "$TMPFILE.warpinfo" | cut -f3- -d'(' | \
+ sed -e 's/^ //' -e 's/, / /' -e 's/)$//'`
+LL_DMS=`grep '^Lower Left ' "$TMPFILE.warpinfo" | cut -f3- -d'(' | \
+ sed -e 's/^ //' -e 's/, / /' -e 's/)$//'`
+UR_DMS=`grep '^Upper Right ' "$TMPFILE.warpinfo" | cut -f3- -d'(' | \
+ sed -e 's/^ //' -e 's/, / /' -e 's/)$//'`
+LR_DMS=`grep '^Lower Right ' "$TMPFILE.warpinfo" | cut -f3- -d'(' | \
+ sed -e 's/^ //' -e 's/, / /' -e 's/)$//'`
+# FIXME: will this work for geographic coords?
+CENTER_DMS=`grep '^Center ' "$TMPFILE.warpinfo" | cut -f3- -d'(' | \
+ sed -e 's/^ //' -e 's/, / /' -e 's/)$//'`
+
+N1=`echo "$UL_DMS" | m.proj -dg proj_in="$IN_PROJ" proj_out="$OUT_PROJ" --quiet | cut -f2 -d'|'`
+N2=`echo "$UR_DMS" | m.proj -dg proj_in="$IN_PROJ" proj_out="$OUT_PROJ" --quiet | cut -f2 -d'|'`
+S1=`echo "$LL_DMS" | m.proj -dg proj_in="$IN_PROJ" proj_out="$OUT_PROJ" --quiet | cut -f2 -d'|'`
+S2=`echo "$LR_DMS" | m.proj -dg proj_in="$IN_PROJ" proj_out="$OUT_PROJ" --quiet | cut -f2 -d'|'`
+E1=`echo "$UR_DMS" | m.proj -dg proj_in="$IN_PROJ" proj_out="$OUT_PROJ" --quiet | cut -f1 -d'|'`
+E2=`echo "$LR_DMS" | m.proj -dg proj_in="$IN_PROJ" proj_out="$OUT_PROJ" --quiet | cut -f1 -d'|'`
+W1=`echo "$UL_DMS" | m.proj -dg proj_in="$IN_PROJ" proj_out="$OUT_PROJ" --quiet | cut -f1 -d'|'`
+W2=`echo "$LL_DMS" | m.proj -dg proj_in="$IN_PROJ" proj_out="$OUT_PROJ" --quiet | cut -f1 -d'|'`
+
+CENTER_DDD=`echo "$CENTER_DMS" | \
+ m.proj -dg proj_in="$IN_PROJ" proj_out="$OUT_PROJ" fs=' ' --quiet | \
+ cut -f1,2 -d ' '`
+
+WARP_N=`echo "$N1 $N2" | awk '{ if ($1 > $2) {print $1} else {print $2} }'`
+WARP_S=`echo "$S1 $S2" | awk '{ if ($1 < $2) {print $1} else {print $2} }'`
+WARP_E=`echo "$E1 $E2" | awk '{ if ($1 > $2) {print $1} else {print $2} }'`
+WARP_W=`echo "$W1 $W2" | awk '{ if ($1 < $2) {print $1} else {print $2} }'`
+# -- parse gdalinfo --
+
+# then this would go in the mapurl file:
+#center=$CENTER_DDD
+#bounds=$WARP_W $WARP_S $WARP_E $WARP_N
+fi
+
+
+# Screen DPI is monitor/device dependent, and therefore a crude one
+# size fits none.
+# a: major radius of Earth (WGS84)
+# pixel factor: meters/pixel for a 256 pixel wide tile
+#
+# zoom_level, map_scale results at 96 dpi, 45deg:
+# [0 1 : 418365887]
+# [1 1 : 209182943]
+# [2 1 : 104591472]
+# [3 1 : 52295736]
+# [4 1 : 26147868]
+# [5 1 : 13073934]
+# [6 1 : 6536967]
+# [7 1 : 3268483]
+# [8 1 : 1634242]
+# [9 1 : 817121]
+# [10 1 : 408560]
+# [11 1 : 204280]
+# [12 1 : 102140]
+# [13 1 : 51070]
+# [14 1 : 25535]
+# [15 1 : 12768]
+# [16 1 : 6384]
+# [17 1 : 3192]
+# [18 1 : 1596]
+# [19 1 : 798]
+# [20 1 : 399]
+#
+calc_webtile_scale()
+{
+ lat=$1
+ scale=$2
+ dpi=$3
+ echo "$lat $scale $dpi" | awk \
+ '{
+ a = 6378137.0;
+ PI = 3.14159265358979323846;
+ pixelfact = a * 2*PI / 256;
+ resolution = pixelfact * cos($1 * PI/180) / 2^$2;
+
+ printf("%.0f\n", $3 * resolution / 0.0254)
+ }'
+}
+
+
+if [ -n "$GIS_OPT_ZOOM" ] ; then
+ ZOOM="$GIS_OPT_ZOOM"
+
+ if [ `echo "$ZOOM" | grep -c '\-'` -gt 0 ] ; then
+ # pick a mid-level zoom for the mapurl file
+ default_zoom=`echo "$ZOOM" | awk -F '-' '{print $1 + int(($2 - $1) / 2)}'`
+ maxzoom=`echo "$ZOOM" | cut -f2 -d'-'`
+ minzoom=`echo "$ZOOM" | cut -f1 -d'-'`
+ else
+ default_zoom="$ZOOM"
+ maxzoom="$ZOOM"
+ minzoom="$ZOOM"
+ fi
+
+else
+ # auto-calc map scale from map extent and dpi
+ eval `g.region -gu`
+ eval `g.region -legu`
+
+ # adjust as needed depending if you are working on a desktop monitor or mobile device:
+ # a value of 96 might be appropriate for desktop, 225 for tablet or smartphone. YMMV.
+ SCREEN_DPI="$GIS_OPT_DPI"
+
+ #SCALE_NS=$ns_extent/(($rows/$DPI)*$INCH2METER)
+ #SCALE_EW=$ew_extent/(($cols/$DPI)*$INCH2METER)
+ SCALE_NS=`echo "$ns_extent $rows $SCREEN_DPI" | awk '{print $1/(($2/$3)*0.0254)}'`
+ SCALE_EW=`echo "$ew_extent $cols $SCREEN_DPI" | awk '{print $1/(($2/$3)*0.0254)}'`
+ # for use with TMS the two above should be the same.
+ MAP_SCALE=`echo "$SCALE_NS $SCALE_EW" | awk '{print int(0.5 + ($1 + $2) / 2)}'`
+
+ # beware zoom levels lower than 9 (~1:500k) badly distort the Mercator
+ minzoom=0
+ maxzoom=3
+
+ for zoom in `seq 0 20` ; do
+ zoomscale=`calc_webtile_scale $center_lat $zoom $SCREEN_DPI`
+ #echo "[$zoom 1 : $zoomscale]" # @ ${SCREEN_DPI}dpi]"
+ if [ "$zoomscale" -gt "$MAP_SCALE" ] ; then
+ maxzoom=`expr "$zoom" + 2`
+ minzoom=`expr "$zoom" - 2`
+ fi
+ done
+
+ if [ "$maxzoom" -gt 20 ] ; then
+ maxzoom=20
+ fi
+ if [ "$minzoom" -lt 0 ] ; then
+ minzoom=0
+ fi
+
+ ZOOM="$minzoom-$maxzoom"
+
+ default_zoom=`expr $maxzoom - 2`
+fi
+
+
+if [ "$MAP_TYPE" = "raster" ] ; then
+ TITLE=`r.info -m "$MAP_NAME" | cut -f2- -d=`
+else
+ TITLE=`r.info -m $(i.group -g $MAP_NAME | head -n 1) | cut -f2- -d=`
+fi
+if [ -z "$TITLE" ] ; then
+ TITLE="$MAP_NAME"
+fi
+
+if [ -n "$GRASS_VERBOSE" ] && [ "$GRASS_VERBOSE" -gt 1 ] ; then
+ VERBOSE="-v"
+else
+ VERBOSE=""
+fi
+
+OUT_DIR="$TMP_DIR"/`basename "$TMP_GMERC" .tif`
+
+g.message "Generating tiles for zoom level(s) $ZOOM"
+
+gdal2tiles.py $VERBOSE -r bilinear -z "$ZOOM" -w none -t "$TITLE" \
+ "$TMP_GMERC" "$OUT_DIR"
+
+if [ $? -ne 0 ] ; then
+ g.message -e "gdal2tiles failed."
+ cleanup
+ exit 1
+fi
+
+rm -f "$TMP_GMERC"
+
+g.message "Dropping unneeded auxiliary metadata files ..."
+
+for file in `find $(dirname "$OUT_DIR") | grep .aux.xml` ; do
+ rm -f "$file"
+done
+
+#### write out .mapurl file
+
+# get max extent bbox in lat/lon
+eval `g.region -bgu`
+
+URLFILE="$OUT_DIR"/`basename "$OUTFILE"`.mapurl
+
+
+if [ "$GIS_OPT_FORMAT" = "jpeg" ] ; then
+ ext=jpg
+else
+ ext=png
+fi
+
+cat << EOF > "$URLFILE"
+url=$OUTFILE/ZZZ/XXX/YYY.$ext
+minzoom=$minzoom
+maxzoom=$maxzoom
+center=$ll_clon $ll_clat
+bounds=$ll_w $ll_s $ll_e $ll_n
+type=tms
+mbtiles=$OUTFILE/$OUTFILE.mbtiles
+name=$TITLE
+description=Exported zoom level(s) $ZOOM from GRASS GIS r.out.mbtiles
+defaultzoom=$default_zoom
+format=$ext
+request_type=run,fill,load
+EOF
+
+cd "$OUT_DIR"
+
+g.message "Scanning for blank tiles ..."
+# warping can leave a lot of left over border around the edges due to skew
+# MBTiles links these all to a single tile, but they're unsightly.
+# This takes a while so consider running in parallel, but it is probably
+# near-saturation already due to the large number of tiny processes creaded
+# and destroyed.
+i=0
+for file in `find . | grep '.png$'` ; do
+ # check if image is empty
+ COUNT=`gdalinfo -mm -noct -nomd "$file" | grep 'Min/Max' | \
+ cut -f2 -d= | tr ',' '\n' | uniq | wc -l`
+ if [ "$COUNT" -eq 1 ] ; then
+ rm -f "$file"
+ i=`expr $i + 1`
+ fi
+done
+if [ "$i" -gt 0 ] ; then
+ g.message "Removed $i blank tiles."
+fi
+
+
+if [ "$GIS_OPT_FORMAT" = "jpeg" ] ; then
+ g.message "Converting to JPEG ..."
+
+ pgms="pngtopnm pnmtojpeg"
+ for pgm in $pgms ; do
+ if [ ! -x "`which $pgm`" ] ; then
+ g.message -e "$pgm is required, please install the NetPBM tools."
+ cleanup
+ exit 1
+ fi
+ done
+
+ for file in `find . | grep '.png$'` ; do
+ outfile=`echo "$file" | sed -e 's/\.png$/.jpg/'`
+ pngtopnm -mix -background '#FFFFFF' "$file" | pnmtojpeg > "$outfile"
+ if [ $? -eq 0 ] ; then
+ rm "$file"
+ fi
+ done
+
+ # sed -i is not portable to BSD/Mac
+ mv tilemapresource.xml tilemapresource.xml.tmp
+ sed -e 's+mime-type="image/png"+mime-type="image/jpeg"+' \
+ -e 's+extension="png"+extension="jpg"+' \
+ tilemapresource.xml.tmp > tilemapresource.xml
+ rm tilemapresource.xml.tmp
+fi
+
+cd - > /dev/null
+
+mv "$OUT_DIR" "$OUTFILE"
+mv "$OUTFILE"/*.mapurl "$OUTFILE"/..
+
+
+#### Set up SQLite commands to build MBTiles database
+# must use sqlite3+
+if [ "$GIS_FLAG_M" -eq 1 ] ; then
+
+ g.message "Creating SQLite instructions ..."
+
+ if [ ! -x "`which sqlite3`" ] ; then
+ g.message -e "sqlite3 is required for creating the MBTiles database, please install it first."
+ cleanup
+ exit 1
+ fi
+
+ TITLE_CLEAN=`echo "$TITLE" | sed -e "s|'|''|g"`
+
+ cat << EOF > "$TMPFILE.sql"
+CREATE TABLE metadata (name text, value text);
+INSERT INTO metadata VALUES ('name', '`basename "$OUTFILE"`');
+INSERT INTO metadata VALUES ('type', 'baselayer');
+INSERT INTO metadata VALUES ('version', '1.1');
+INSERT INTO metadata VALUES ('description', '$TITLE_CLEAN');
+INSERT INTO metadata VALUES ('format', '$ext');
+INSERT INTO metadata VALUES ('bounds', '$ll_w,$ll_s,$ll_e,$ll_n');
+INSERT INTO metadata VALUES ('center', '$ll_clon,$ll_clat');
+INSERT INTO metadata VALUES ('minzoom', '$minzoom');
+INSERT INTO metadata VALUES ('maxzoom', '$maxzoom');
+INSERT INTO metadata VALUES ('tile_row_type', 'tms');
+CREATE TABLE tiles (zoom_level integer, tile_column integer, tile_row integer, tile_data blob);
+BEGIN TRANSACTION;
+EOF
+
+
+ # FIXME: keep DB small by checking for duplicate tile images (ocean, wasteland)
+ # Make the 'tiles' table a view, then replace 'blob' with (tile_id text),
+ # and in the 'images' table keep unique tiles in (tile_data blob, tile_id text).
+ # perhaps try 'md5sum | sort | uniq -d' to weed out duplicates?
+ # see https://github.com/geopaparazzi/geopaparazzi/wiki/mbtiles-Implementation
+ (
+ for file in `find "$OUTFILE" -name "*.$ext" | sort -n -t / -k 2` ; do
+ zoomcolrow=`echo "$file" | sed -e "s|.*$OUTFILE/||" -e "s|\.$ext||" -e 's|/|, |g'`
+
+ echo "INSERT INTO tiles VALUES ($zoomcolrow, X'`hexdump -v -e '1/1 "%.2x"' $file`');"
+ done
+ ) >> "$TMPFILE.sql"
+
+ echo "COMMIT TRANSACTION;" >> "$TMPFILE.sql"
+
+ if [ `wc -l < "${TMPFILE}_tiles.sql"` -gt 65535 ] ; then
+ echo "FIXME: insert some 'COMMIT TRANSACTION; BEGIN TRANSACTION;' inline the sql every 65535 inserts or so"
+ fi
+
+ g.message "Populating SQLite database ..."
+ cat "$TMPFILE.sql" | sqlite3 "$OUTFILE.mbtiles"
+
+fi
+
+
+cleanup
+
+if [ "$GIS_FLAG_T" -eq 1 ] ; then
+ if [ ! -x "`which tar`" ] ; then
+ g.message -e "tar is required, please install it first."
+ exit 1
+ fi
+
+ cd `dirname "$OUTFILE"`
+ tar cf `basename "$OUTFILE"`.tar `basename "$OUTFILE"` `basename "$OUTFILE"`.mapurl
+fi
+
+
+g.message message=""
+g.message "Done."
+
+#g.message "Next copy the '$OUTFILE' directory and mapurl file to \
+# /sdcard/maps/ on your Android device and select the mapurl file as the \
+# tile source in Geopaparazzi. The MBTiles SQLite database will then \
+# be populated in the background."
+
Deleted: grass-addons/grass6/raster/r.out.mbtiles/r.out.mbtiles_prep
===================================================================
--- grass-addons/grass6/raster/r.out.mbtiles/r.out.mbtiles_prep 2014-05-20 05:04:11 UTC (rev 60363)
+++ grass-addons/grass6/raster/r.out.mbtiles/r.out.mbtiles_prep 2014-05-20 05:06:00 UTC (rev 60364)
@@ -1,550 +0,0 @@
-#!/bin/sh
-############################################################################
-#
-# MODULE: r.out.mbtiles
-# AUTHOR(S): M. Hamish Bowman, Dunedin, New Zealand
-# PURPOSE: Export a raster map into a TMS bundle ready for conversion
-# to MBTiles SQLite format
-# COPYRIGHT: (C) 2014 Hamish Bowman, and the GRASS Development Team
-#
-# This program is free software; you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation; either version 2 of the License, 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.
-#
-# *** NOT FOR NAVIGATIONAL USE - USE AT YOUR OWN RISK ***
-#
-############################################################################
-
-# https://github.com/geopaparazzi/geopaparazzi/wiki/.mapurl-Converting-Tile-Set-to-.mbtiles
-# https://github.com/geopaparazzi/geopaparazzi/wiki/.mapurl-parameters
-# https://github.com/mapbox/mbtiles-spec/blob/master/1.1/spec.md
-# https://github.com/geopaparazzi/geopaparazzi/wiki/mbtiles-Implementation
-
-#%Module
-#% description: Export GRASS raster as a TMS tree ready for converting to mbtiles format.
-#% keywords: raster, export, tms, mbtiles
-#%End
-#%Option
-#% key: input
-#% type: string
-#% required: yes
-#% multiple: no
-#% key_desc: name
-#% description: Name of input raster map
-#% gisprompt: old,cell,raster
-#%End
-#%Option
-#% key: output
-#% type: string
-#% required: no
-#% multiple: no
-#% key_desc: name
-#% label: Name for output project directory
-#% description: default: input map name
-#% gisprompt: new_file,file,output
-#%End
-#%Option
-#% key: format
-#% type: string
-#% required: no
-#% multiple: no
-#% key_desc: string
-#% description: Image format to store the tiles
-#% options: png,jpeg
-#% answer: png
-#%End
-#%Option
-#% key: zoom
-#% type: string
-#% required: no
-#% multiple: no
-#% key_desc: string
-#% label: TMS zoom levels to create
-#% description: Format: '0-19', or '10' (default: automatic)
-#%End
-#%Option
-#% key: dpi
-#% type: integer
-#% label: DPI of the target display
-#% description: Value for a smartphone might be 225, a desktop perhaps 96.
-#% answer: 225
-#%End
-#%Flag
-#% key: m
-#% description: Create a MBTiles database (requires SQLite3)
-#%End
-#%Flag
-#% key: t
-#% description: Create a tarball of the tile tree
-#%End
-
-
-
-if [ -z "$GISBASE" ] ; then
- echo "You must be in GRASS GIS to run this program." >&2
- exit 1
-fi
-
-if [ "$1" != "@ARGS_PARSED@" ] ; then
- exec g.parser "$0" "$@"
-fi
-
-#### check if we have awk, gdalwarp, gdal2tiles, seq
-pgms="awk seq gdalwarp gdalinfo gdal2tiles.py"
-for pgm in $pgms ; do
- if [ ! -x "`which $pgm`" ] ; then
- g.message -e "$pgm is required, please install it first."
- exit 1
- fi
-done
-
-
-MAP_NAME="$GIS_OPT_INPUT"
-
-if [ -n "$GIS_OPT_OUTPUT" ] ; then
- OUTFILE="$GIS_OPT_OUTPUT"
-else
- OUTFILE=`echo "$MAP_NAME" | cut -f1 -d'@'`
-fi
-
-
-g.findfile element=cell file="$MAP_NAME" > /dev/null
-if [ $? -eq 0 ] ; then
- MAP_TYPE=raster
-else
- g.findfile element=group file="$MAP_NAME" > /dev/null
- if [ $? -eq 0 ] ; then
- MAP_TYPE=group
- else
- g.message -e "Could not find <$MAP_NAME>"
- exit 1
- fi
-fi
-
-cleanup()
-{
- #remove temporary files
- if [ -e "$TMPFILE" ] ; then
- rm -rf "$TMPFILE"*
- fi
-}
-
-# what to do in case of user break:
-exitprocedure()
-{
- g.message -e message='User break!'
- cleanup
- exit 1
-}
-
-# shell check for user break (signal list: trap -l)
-trap "exitprocedure" 2 3 15
-
-
-#### setup temporary file
-TMPFILE="`g.tempfile pid=$$`"
-if [ $? -ne 0 ] || [ -z "$TMPFILE" ] ; then
- g.message -e "Unable to create temporary file"
- exit 1
-fi
-
-
-TMP_GTIFF="${TMPFILE}_export.tif"
-TMP_DIR=`dirname "$TMP_GTIFF"`
-TMP_GMERC="$TMP_DIR"/`basename "$TMP_GTIFF" .tif`_gmerc.tif
-
-
-if [ -n "$GRASS_VERBOSE" ] && [ "$GRASS_VERBOSE" -eq 0 ] ; then
- QUIET="-q"
-else
- QUIET=""
-fi
-
-# ?? do the source of TMS tiles sometimes need to be reprojected into the
-# funny google merc projection first?
-
-if [ "$MAP_TYPE" = "raster" ] ; then
- r.out.tiff -t input="$MAP_NAME" output="$TMP_GTIFF"
- #gdalwarp -s_srs "`g.proj -jf`" -t_srs EPSG:900913 -dstalpha "$TMP_GTIFF" "$TMP_GMERC"
- gdal_translate $QUIET -a_srs "`g.proj -jf`" "$TMP_GTIFF" "$TMP_GMERC"
-else
- if [ `echo "$MAP_NAME" | grep -c '@'` -gt 0 ] ; then
- g.message -e "Imagery groups must exist in the current mapset, sorry."
- cleanup
- exit 1
- fi
-
- r.out.gdal input="$MAP_NAME" output="$TMP_GTIFF" format=GTiff
- #gdalwarp -t_srs EPSG:900913 -dstalpha "$TMP_GTIFF" "$TMP_GMERC"
- mv "$TMP_GTIFF" "$TMP_GMERC"
-fi
-## gdalwarp -srcnodata -dstnodata -dstalpha ?
-
-if [ $? -ne 0 ] ; then
- g.message -e "Problem exporting image"
- cleanup
- exit 1
-fi
-
-# free the disk space ASAP
-rm -f "$TMP_GTIFF"
-
-
-if [ 0 -eq 1 ] ; then
-
-## I don't think we need to do the gdalwarp, but it's worth keeping this
-## next bit of code for parsing the lat/lon bbox out of gdalinfo for a
-## rainy day.
-
-#### we need to find the bounding box in lat/lon of the post-warped geotiff
-# -- parse gdalinfo --
-gdalinfo -noct -nomd "$TMP_GMERC" > "$TMPFILE.warpinfo"
-
-IN_PROJ="+proj=longlat +datum=WGS84"
-OUT_PROJ="+proj=longlat +datum=WGS84"
-#?OUT_PROJ="+proj=longlat +ellps=sphere +a=6378137 +b=6378137 +nadgrids=@null"
-
-
-UL_DMS=`grep '^Upper Left ' "$TMPFILE.warpinfo" | cut -f3- -d'(' | \
- sed -e 's/^ //' -e 's/, / /' -e 's/)$//'`
-LL_DMS=`grep '^Lower Left ' "$TMPFILE.warpinfo" | cut -f3- -d'(' | \
- sed -e 's/^ //' -e 's/, / /' -e 's/)$//'`
-UR_DMS=`grep '^Upper Right ' "$TMPFILE.warpinfo" | cut -f3- -d'(' | \
- sed -e 's/^ //' -e 's/, / /' -e 's/)$//'`
-LR_DMS=`grep '^Lower Right ' "$TMPFILE.warpinfo" | cut -f3- -d'(' | \
- sed -e 's/^ //' -e 's/, / /' -e 's/)$//'`
-# FIXME: will this work for geographic coords?
-CENTER_DMS=`grep '^Center ' "$TMPFILE.warpinfo" | cut -f3- -d'(' | \
- sed -e 's/^ //' -e 's/, / /' -e 's/)$//'`
-
-N1=`echo "$UL_DMS" | m.proj -dg proj_in="$IN_PROJ" proj_out="$OUT_PROJ" --quiet | cut -f2 -d'|'`
-N2=`echo "$UR_DMS" | m.proj -dg proj_in="$IN_PROJ" proj_out="$OUT_PROJ" --quiet | cut -f2 -d'|'`
-S1=`echo "$LL_DMS" | m.proj -dg proj_in="$IN_PROJ" proj_out="$OUT_PROJ" --quiet | cut -f2 -d'|'`
-S2=`echo "$LR_DMS" | m.proj -dg proj_in="$IN_PROJ" proj_out="$OUT_PROJ" --quiet | cut -f2 -d'|'`
-E1=`echo "$UR_DMS" | m.proj -dg proj_in="$IN_PROJ" proj_out="$OUT_PROJ" --quiet | cut -f1 -d'|'`
-E2=`echo "$LR_DMS" | m.proj -dg proj_in="$IN_PROJ" proj_out="$OUT_PROJ" --quiet | cut -f1 -d'|'`
-W1=`echo "$UL_DMS" | m.proj -dg proj_in="$IN_PROJ" proj_out="$OUT_PROJ" --quiet | cut -f1 -d'|'`
-W2=`echo "$LL_DMS" | m.proj -dg proj_in="$IN_PROJ" proj_out="$OUT_PROJ" --quiet | cut -f1 -d'|'`
-
-CENTER_DDD=`echo "$CENTER_DMS" | \
- m.proj -dg proj_in="$IN_PROJ" proj_out="$OUT_PROJ" fs=' ' --quiet | \
- cut -f1,2 -d ' '`
-
-WARP_N=`echo "$N1 $N2" | awk '{ if ($1 > $2) {print $1} else {print $2} }'`
-WARP_S=`echo "$S1 $S2" | awk '{ if ($1 < $2) {print $1} else {print $2} }'`
-WARP_E=`echo "$E1 $E2" | awk '{ if ($1 > $2) {print $1} else {print $2} }'`
-WARP_W=`echo "$W1 $W2" | awk '{ if ($1 < $2) {print $1} else {print $2} }'`
-# -- parse gdalinfo --
-
-# then this would go in the mapurl file:
-#center=$CENTER_DDD
-#bounds=$WARP_W $WARP_S $WARP_E $WARP_N
-fi
-
-
-# Screen DPI is monitor/device dependent, and therefore a crude one
-# size fits none.
-# a: major radius of Earth (WGS84)
-# pixel factor: meters/pixel for a 256 pixel wide tile
-#
-# zoom_level, map_scale results at 96 dpi, 45deg:
-# [0 1 : 418365887]
-# [1 1 : 209182943]
-# [2 1 : 104591472]
-# [3 1 : 52295736]
-# [4 1 : 26147868]
-# [5 1 : 13073934]
-# [6 1 : 6536967]
-# [7 1 : 3268483]
-# [8 1 : 1634242]
-# [9 1 : 817121]
-# [10 1 : 408560]
-# [11 1 : 204280]
-# [12 1 : 102140]
-# [13 1 : 51070]
-# [14 1 : 25535]
-# [15 1 : 12768]
-# [16 1 : 6384]
-# [17 1 : 3192]
-# [18 1 : 1596]
-# [19 1 : 798]
-# [20 1 : 399]
-#
-calc_webtile_scale()
-{
- lat=$1
- scale=$2
- dpi=$3
- echo "$lat $scale $dpi" | awk \
- '{
- a = 6378137.0;
- PI = 3.14159265358979323846;
- pixelfact = a * 2*PI / 256;
- resolution = pixelfact * cos($1 * PI/180) / 2^$2;
-
- printf("%.0f\n", $3 * resolution / 0.0254)
- }'
-}
-
-
-if [ -n "$GIS_OPT_ZOOM" ] ; then
- ZOOM="$GIS_OPT_ZOOM"
-
- if [ `echo "$ZOOM" | grep -c '\-'` -gt 0 ] ; then
- # pick a mid-level zoom for the mapurl file
- default_zoom=`echo "$ZOOM" | awk -F '-' '{print $1 + int(($2 - $1) / 2)}'`
- maxzoom=`echo "$ZOOM" | cut -f2 -d'-'`
- minzoom=`echo "$ZOOM" | cut -f1 -d'-'`
- else
- default_zoom="$ZOOM"
- maxzoom="$ZOOM"
- minzoom="$ZOOM"
- fi
-
-else
- # auto-calc map scale from map extent and dpi
- eval `g.region -gu`
- eval `g.region -legu`
-
- # adjust as needed depending if you are working on a desktop monitor or mobile device:
- # a value of 96 might be appropriate for desktop, 225 for tablet or smartphone. YMMV.
- SCREEN_DPI="$GIS_OPT_DPI"
-
- #SCALE_NS=$ns_extent/(($rows/$DPI)*$INCH2METER)
- #SCALE_EW=$ew_extent/(($cols/$DPI)*$INCH2METER)
- SCALE_NS=`echo "$ns_extent $rows $SCREEN_DPI" | awk '{print $1/(($2/$3)*0.0254)}'`
- SCALE_EW=`echo "$ew_extent $cols $SCREEN_DPI" | awk '{print $1/(($2/$3)*0.0254)}'`
- # for use with TMS the two above should be the same.
- MAP_SCALE=`echo "$SCALE_NS $SCALE_EW" | awk '{print int(0.5 + ($1 + $2) / 2)}'`
-
- # beware zoom levels lower than 9 (~1:500k) badly distort the Mercator
- minzoom=0
- maxzoom=3
-
- for zoom in `seq 0 20` ; do
- zoomscale=`calc_webtile_scale $center_lat $zoom $SCREEN_DPI`
- #echo "[$zoom 1 : $zoomscale]" # @ ${SCREEN_DPI}dpi]"
- if [ "$zoomscale" -gt "$MAP_SCALE" ] ; then
- maxzoom=`expr "$zoom" + 2`
- minzoom=`expr "$zoom" - 2`
- fi
- done
-
- if [ "$maxzoom" -gt 20 ] ; then
- maxzoom=20
- fi
- if [ "$minzoom" -lt 0 ] ; then
- minzoom=0
- fi
-
- ZOOM="$minzoom-$maxzoom"
-
- default_zoom=`expr $maxzoom - 2`
-fi
-
-
-if [ "$MAP_TYPE" = "raster" ] ; then
- TITLE=`r.info -m "$MAP_NAME" | cut -f2- -d=`
-else
- TITLE=`r.info -m $(i.group -g $MAP_NAME | head -n 1) | cut -f2- -d=`
-fi
-if [ -z "$TITLE" ] ; then
- TITLE="$MAP_NAME"
-fi
-
-if [ -n "$GRASS_VERBOSE" ] && [ "$GRASS_VERBOSE" -gt 1 ] ; then
- VERBOSE="-v"
-else
- VERBOSE=""
-fi
-
-OUT_DIR="$TMP_DIR"/`basename "$TMP_GMERC" .tif`
-
-g.message "Generating tiles for zoom level(s) $ZOOM"
-
-gdal2tiles.py $VERBOSE -r bilinear -z "$ZOOM" -w none -t "$TITLE" \
- "$TMP_GMERC" "$OUT_DIR"
-
-if [ $? -ne 0 ] ; then
- g.message -e "gdal2tiles failed."
- cleanup
- exit 1
-fi
-
-rm -f "$TMP_GMERC"
-
-g.message "Dropping unneeded auxiliary metadata files ..."
-
-for file in `find $(dirname "$OUT_DIR") | grep .aux.xml` ; do
- rm -f "$file"
-done
-
-#### write out .mapurl file
-
-# get max extent bbox in lat/lon
-eval `g.region -bgu`
-
-URLFILE="$OUT_DIR"/`basename "$OUTFILE"`.mapurl
-
-
-if [ "$GIS_OPT_FORMAT" = "jpeg" ] ; then
- ext=jpg
-else
- ext=png
-fi
-
-cat << EOF > "$URLFILE"
-url=$OUTFILE/ZZZ/XXX/YYY.$ext
-minzoom=$minzoom
-maxzoom=$maxzoom
-center=$ll_clon $ll_clat
-bounds=$ll_w $ll_s $ll_e $ll_n
-type=tms
-mbtiles=$OUTFILE/$OUTFILE.mbtiles
-name=$TITLE
-description=Exported zoom level(s) $ZOOM from GRASS GIS r.out.mbtiles
-defaultzoom=$default_zoom
-format=$ext
-request_type=run,fill,load
-EOF
-
-cd "$OUT_DIR"
-
-g.message "Scanning for blank tiles ..."
-# warping can leave a lot of left over border around the edges due to skew
-# MBTiles links these all to a single tile, but they're unsightly.
-# This takes a while so consider running in parallel, but it is probably
-# near-saturation already due to the large number of tiny processes creaded
-# and destroyed.
-i=0
-for file in `find . | grep '.png$'` ; do
- # check if image is empty
- COUNT=`gdalinfo -mm -noct -nomd "$file" | grep 'Min/Max' | \
- cut -f2 -d= | tr ',' '\n' | uniq | wc -l`
- if [ "$COUNT" -eq 1 ] ; then
- rm -f "$file"
- i=`expr $i + 1`
- fi
-done
-if [ "$i" -gt 0 ] ; then
- g.message "Removed $i blank tiles."
-fi
-
-
-if [ "$GIS_OPT_FORMAT" = "jpeg" ] ; then
- g.message "Converting to JPEG ..."
-
- pgms="pngtopnm pnmtojpeg"
- for pgm in $pgms ; do
- if [ ! -x "`which $pgm`" ] ; then
- g.message -e "$pgm is required, please install the NetPBM tools."
- cleanup
- exit 1
- fi
- done
-
- for file in `find . | grep '.png$'` ; do
- outfile=`echo "$file" | sed -e 's/\.png$/.jpg/'`
- pngtopnm -mix -background '#FFFFFF' "$file" | pnmtojpeg > "$outfile"
- if [ $? -eq 0 ] ; then
- rm "$file"
- fi
- done
-
- # sed -i is not portable to BSD/Mac
- mv tilemapresource.xml tilemapresource.xml.tmp
- sed -e 's+mime-type="image/png"+mime-type="image/jpeg"+' \
- -e 's+extension="png"+extension="jpg"+' \
- tilemapresource.xml.tmp > tilemapresource.xml
- rm tilemapresource.xml.tmp
-fi
-
-cd - > /dev/null
-
-mv "$OUT_DIR" "$OUTFILE"
-mv "$OUTFILE"/*.mapurl "$OUTFILE"/..
-
-
-#### Set up SQLite commands to build MBTiles database
-# must use sqlite3+
-if [ "$GIS_FLAG_M" -eq 1 ] ; then
-
- g.message "Creating SQLite instructions ..."
-
- if [ ! -x "`which sqlite3`" ] ; then
- g.message -e "sqlite3 is required for creating the MBTiles database, please install it first."
- cleanup
- exit 1
- fi
-
- TITLE_CLEAN=`echo "$TITLE" | sed -e "s|'|''|g"`
-
- cat << EOF > "$TMPFILE.sql"
-CREATE TABLE metadata (name text, value text);
-INSERT INTO metadata VALUES ('name', '`basename "$OUTFILE"`');
-INSERT INTO metadata VALUES ('type', 'baselayer');
-INSERT INTO metadata VALUES ('version', '1.1');
-INSERT INTO metadata VALUES ('description', '$TITLE_CLEAN');
-INSERT INTO metadata VALUES ('format', '$ext');
-INSERT INTO metadata VALUES ('bounds', '$ll_w,$ll_s,$ll_e,$ll_n');
-INSERT INTO metadata VALUES ('center', '$ll_clon,$ll_clat');
-INSERT INTO metadata VALUES ('minzoom', '$minzoom');
-INSERT INTO metadata VALUES ('maxzoom', '$maxzoom');
-INSERT INTO metadata VALUES ('tile_row_type', 'tms');
-CREATE TABLE tiles (zoom_level integer, tile_column integer, tile_row integer, tile_data blob);
-BEGIN TRANSACTION;
-EOF
-
-
- # FIXME: keep DB small by checking for duplicate tile images (ocean, wasteland)
- # Make the 'tiles' table a view, then replace 'blob' with (tile_id text),
- # and in the 'images' table keep unique tiles in (tile_data blob, tile_id text).
- # perhaps try 'md5sum | sort | uniq -d' to weed out duplicates?
- # see https://github.com/geopaparazzi/geopaparazzi/wiki/mbtiles-Implementation
- (
- for file in `find "$OUTFILE" -name "*.$ext" | sort -n -t / -k 2` ; do
- zoomcolrow=`echo "$file" | sed -e "s|.*$OUTFILE/||" -e "s|\.$ext||" -e 's|/|, |g'`
-
- echo "INSERT INTO tiles VALUES ($zoomcolrow, X'`hexdump -v -e '1/1 "%.2x"' $file`');"
- done
- ) >> "$TMPFILE.sql"
-
- echo "COMMIT TRANSACTION;" >> "$TMPFILE.sql"
-
- if [ `wc -l < "${TMPFILE}_tiles.sql"` -gt 65535 ] ; then
- echo "FIXME: insert some 'COMMIT TRANSACTION; BEGIN TRANSACTION;' inline the sql every 65535 inserts or so"
- fi
-
- g.message "Populating SQLite database ..."
- cat "$TMPFILE.sql" | sqlite3 "$OUTFILE.mbtiles"
-
-fi
-
-
-cleanup
-
-if [ "$GIS_FLAG_T" -eq 1 ] ; then
- if [ ! -x "`which tar`" ] ; then
- g.message -e "tar is required, please install it first."
- exit 1
- fi
-
- cd `dirname "$OUTFILE"`
- tar cf `basename "$OUTFILE"`.tar `basename "$OUTFILE"` `basename "$OUTFILE"`.mapurl
-fi
-
-
-g.message message=""
-g.message "Done."
-
-#g.message "Next copy the '$OUTFILE' directory and mapurl file to \
-# /sdcard/maps/ on your Android device and select the mapurl file as the \
-# tile source in Geopaparazzi. The MBTiles SQLite database will then \
-# be populated in the background."
-
More information about the grass-commit
mailing list