r.tileset pythonized

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@37171 15284696-431f-4ddb-bdfa-cd5b030d7da7

scripts/r.tileset/r.tileset

@@ -1,564 +0,0 @@
-#!/bin/bash
-
-############################################################################
-#
-# MODULE:       r.tileset for GRASS 6
-# AUTHOR(S):    Cedric Shock
-# PURPOSE:      To produce tilings of regions in other projections.
-# COPYRIGHT:    (C) 2006 by Cedric Shoc
-#
-#               This program is free software under the GNU General Public
-#               License (>=v2). Read the file COPYING that comes with GRASS
-#               for details.
-#
-#############################################################################
-
-# Requires:
-#  bc
-#  cs2cs
-#  grep
-#  sed
-
-# Bugs:
-#  Does not know about meridians in projections. However, unlike the usual
-#   hack used to find meridians, this code is perfectly happy with arbitrary
-#   rotations and flips
-#  The following are planned to be fixed in a future version, if it turns out
-#   to be necessary for someone:
-#   Does not generate optimal tilings. This means that between an appropriate
-#    projection for a region and latitude longitude projection, in the
-#    degenerate case, it may create tiles demanding up to twice the necessary
-#    information. Requesting data from cylindrical projections near their poles
-#    results in divergence. You really don't want to use source data from
-#    someone who put it in a cylindrical projection near a pole, do you?
-#   Not generating "optimal" tilings has another side effect; the sanity
-#    of the destination region will not carry over to generating tiles of
-#    realistic aspect ratio. This might be a problem for some WMS servers
-#    presenting data in a highly inappropriate projection. Do you really
-#    want their data?
-
-# Usage example
-#  r.tileset -w sourceproj="+init=epsg:4326" maxcols=4096 maxrows=4096
-
-# Changes:
-# Version 2:
-#   Added destination projection options
-#   Added extra cells for overlap
-
-#%Module
-#% description: Produces tilings of the source projection for use in the destination region and projection.
-#% keywords: raster, tiling
-#%End
-#%flag
-#% key: g
-#% description: Produces shell script output
-#%end
-#%flag
-#% key: w
-#% description: Produces web map server query string output
-#%end
-#%option
-#% key: region
-#% type: string
-#% description: Name of region to use instead of current region for bounds and resolution
-#%end
-#%option
-#% key: sourceproj
-#% type: string
-#% description: Source projection
-#% required : yes
-#%end
-#%option
-#% key: sourcescale
-#% type: string
-#% description: Conversion factor from units to meters in source projection
-#% answer : 1
-#%end
-#%option
-#% key: destproj
-#% type: string
-#% description: Destination projection, defaults to this location's projection
-#% required : no
-#%end
-#%option
-#% key: destscale
-#% type: string
-#% description: Conversion factor from units to meters in source projection
-#% required : no
-#%end
-#%option
-#% key: maxcols
-#% type: integer
-#% description: Maximum number of columns for a tile in the source projection
-#% answer: 1024
-#%end
-#%option
-#% key: maxrows
-#% type: integer
-#% description: Maximum number of rows for a tile in the source projection
-#% answer: 1024
-#%end
-#%option
-#% key: overlap
-#% type: integer
-#% description: Number of cells tiles should overlap in each direction
-#% answer: 0
-#%end
-#%option
-#% key: fs
-#% type: string
-#% description: Output field separator
-#% answer:|
-#%end
-#%option
-#% key: v
-#% type: integer
-#% description: Verbosity level
-#% answer: 0
-#%end
-
-if  [ -z "$GISBASE" ] ; then
-    echo "You must be in GRASS GIS to run this program." 1>&2
-    exit 1
-fi
-
-if [ "$1" != "@ARGS_PARSED@" ] ; then
-  exec g.parser "$0" "$@"
-fi
-
-# Program locations:
-BC="bc"
-BCARGS="-l"
-SED="sed"
-GREP="grep"
-CS2CS="cs2cs"
-
-
-# check if we have bc
-if [ ! -x "`which $BC`" ] ; then
-    g.message -e "$BC is required, please install it first"
-    exit 1
-fi
-
-# check if we have sed
-if [ ! -x "`which $SED`" ] ; then
-    g.message -e "$SED is required, please install it first"
-    exit 1
-fi
-
-# check if we have grep
-if [ ! -x "`which $GREP`" ] ; then
-    g.message -e "$GREP is required, please install it first"
-    exit 1
-fi
-
-# check if we have cs2cs
-if [ ! -x "`which $CS2CS`" ] ; then
-    g.message -e "$CS2CS is required, please install it first"
-    exit 1
-fi
-
-
-# Default IFS
-defaultIFS=$IFS
-
-# Data structures used in this program:
-# A bounding box:
-# 0 -> left, 1-> bottom, 2->right, 3-> top
-# A border:
-# An array of points indexed by 0 for "x" and 4 for "y" + by number 0, 1, 2, and 3
-# A reprojector [0] is name of source projection, [1] is name of destination
-# A projection - [0] is proj.4 text, [1] is scale
-
-#####################
-# name:     message
-# purpose:  displays messages to the user
-# usage: message level text
-
-message () {
-	if [ $1 -lt $GIS_OPT_V ] ; then
-		shift
-		echo "$@"  
-	fi
-}
-
-#########################
-# name: lookup
-# purpose: lookup values in passed by reference arrays
-
-Lookup_Mac() {
-# Assignment Command Statement Builder
-# "$1" is Source name, "$3" is Destination name
-    echo "eval $3"'=${'"$1[$2]}"
-}
-
-# Lookup_Mac
-
-declare -f LookupP               # Function "Pointer"
-LookupP=Lookup_Mac               # Statement Builder
-
-lookup () {
-	$($LookupP $1 $2 $3)
-}
-
-#####################3
-# name: calculate
-# purpose: perform calculations
-# usage: varname "expr"
-
-calculate() {
-	message 3 "$2"
-	c_tmp=`echo "$2" | $BC $BCARGS`
-	eval $1=$c_tmp
-}
-
-#########################
-# name: makeProjector
-# purpose: bundle up reprojection information
-#          and provide enough abstraction to use pipes
-# usage: makeReprojector source source-scale destination destination-scale name
-
-makeProjector() {
-	eval $5[0]=$1
-	eval $5[1]=$2
-	eval $5[2]=$3
-	eval $5[3]=$4
-}
-
-#######################################################################
-# name:     project
-# purpose:  projects x1 y1 using projector
-# usage: project -12 36 x2 y2 source-to-dest
-project() {
-	# echo "$5"
-	# Pick apart the projector
-	lookup $5 0 p_source_proj_p
-	lookup $5 1 p_source_units
-	lookup $5 2 p_dest_proj_p
-	lookup $5 3 p_dest_units
-
-	eval "p_source_proj=\$$p_source_proj_p"
-	eval "p_dest_proj=\$$p_dest_proj_p"
-
-	calculate p_x1 "$1 * $p_source_units"
-	calculate p_y1 "$2 * $p_source_units"
-
-	message 3 "echo \"$p_x1 $p_y1\" | $CS2CS -f \"%.8f\" $p_source_proj +to $p_dest_proj"
-
-	answer=`echo "$p_x1 $p_y1" | $CS2CS -f "%.8f" $p_source_proj +to $p_dest_proj`
-
-	if [ $? -ne 0 ] ; then
-	    g.message -e message="Problem running $CS2CS. <$answer>"
-	    exit 1
-	fi
-
-	message 3 "$answer"
-	# do not quote $answer in the following line
-	eval `echo $answer | $SED -e "s/^/p_x2=/"  -e "s/ /;p_y2=/"  -e "s/ /;p_z2=/"`
-
-	calculate p_x2 "$p_x2 / $p_dest_units"
-	calculate p_y2 "$p_y2 / $p_dest_units"
-	
-	eval $3=$p_x2
-	eval $4=$p_y2
-}
-
-min() {
-	lookup "#$1" @ min_n
-	lookup $1 0 min_a
-	for ((min_i=0;min_i<min_n;min_i+=1)) ; do
-		lookup $1 $min_i min_b
-		calculate min_a "($min_b < $min_a) * $min_b + ($min_a <= $min_b) * $min_a" 
-	done
-	eval "$2=$min_a"
-}
-
-max() {
-	lookup "#$1" @ max_n
-	lookup $1 0 max_a
-	for ((max_i=0;max_i<max_n;max_i+=1)) ; do
-		lookup $1 $max_i max_b
-		calculate max_a "($max_b > $max_a) * $max_b + ($max_a >= $max_b) * $max_a" 
-	done
-	eval "$2=$max_a"
-}
-
-########################
-# name:    bboxToPoints
-# purpose: make points that are the corners of a bounding box
-
-
-bboxToPoints() {
-	lookup $1 0 $2[0]
-	lookup $1 1 $2[1]
-	lookup $1 0 $2[2]
-	lookup $1 3 $2[3]
-	lookup $1 2 $2[4]
-	lookup $1 3 $2[5]
-	lookup $1 2 $2[6]
-	lookup $1 1 $2[7]
-}
-
-pointsToBbox() {
-	lookup $1 0 ptb_xs[0]
-	lookup $1 2 ptb_xs[1]
-	lookup $1 4 ptb_xs[2]
-	lookup $1 6 ptb_xs[3]
-	lookup $1 1 ptb_ys[0]
-	lookup $1 3 ptb_ys[1]
-	lookup $1 5 ptb_ys[2]
-	lookup $1 7 ptb_ys[3]
-	min ptb_xs $2[0] 
-	min ptb_ys $2[1]
-	max ptb_xs $2[2]
-	max ptb_ys $2[3]
-}
-
-#######################################################################
-# name:     projectPoints
-# purpose:  projects a list of points
-#
-projectPoints() {
-	lookup "#$1" @ pp_max
-	for ((pp_i=0;pp_i<pp_max;pp_i+=2)) ; do
-		calculate pp_i2 "$pp_i + 1"
-		lookup $1 $pp_i pp_x1
-		lookup $1 $pp_i2 pp_y1
-		# echo "project $pp_x1 $pp_y1 $2[$pp_i] $2[$pp_i2] $3" 1>&2
-		project $pp_x1 $pp_y1 $2[$pp_i] $2[$pp_i2] $3
-	done
-}
-
-######################
-# name:     Side Length
-# purpose:  Find the length of sides of a set of points from one to the next
-# usage: sideLengths points xmetric ymetric answer
-
-sideLengths() {
-	lookup "#$1" @ sl_max
-	for ((sl_i=0;sl_i<sl_max;sl_i+=2)) ; do
-		calculate sl_i2 "$sl_i + 1"
-		calculate sl_j "scale=0; ($sl_i + 2) % $sl_max"
-		calculate sl_j2 "$sl_j + 1"
-		calculate sl_k "scale=0; $sl_i / 2"
-		# echo "$sl_i, $sl_i2, $sl_j, $sl_j2"
-		lookup $1 $sl_i sl_x1
-		lookup $1 $sl_i2 sl_y1
-		lookup $1 $sl_j sl_x2
-		lookup $1 $sl_j2 sl_y2
-		sl_x="(($sl_x2 - $sl_x1) * $2)"
-		sl_y="(($sl_y2 - $sl_y1) * $3)"
-		# echo "$sl_x, $sl_y"
-		calculate sl_d "sqrt ( $sl_x * $sl_x + $sl_y * $sl_y )"
-		# echo "$sl_i ($sl_k): $sl_d"
-		eval "$4[$sl_k]=$sl_d"
-	done
-}
-
-
-#######################################################################
-# name:     bboxesIntersect
-# purpose:  determine if two bounding boxes intersect
-# usage: bboxesIntersect box1 box2 returnvar
-
-bboxesIntersect() {
-	lookup $1 0 bi_A1[0]
-	lookup $1 1 bi_A1[1]
-	lookup $1 2 bi_A2[0]
-	lookup $1 3 bi_A2[1]
-	lookup $2 0 bi_B1[0]
-	lookup $2 1 bi_B1[1]
-	lookup $2 2 bi_B2[0]
-	lookup $2 3 bi_B2[1]
-
-	# Fails at meridians
-	# I have no way to have general knowledge of meridians
-	# This means rewrite in C I imagine
-	for bi_i in 0 1 ; do
-		calculate IN[$bi_i] \
-		"(${bi_A1[$bi_i]} <= ${bi_B1[$bi_i]} && ${bi_A2[$bi_i]} >= ${bi_B2[$bi_i]}) || \
-		 (${bi_A1[$bi_i]} >= ${bi_B1[$bi_i]} && ${bi_A2[$bi_i]} <= ${bi_B2[$bi_i]}) || \
-		 (${bi_A1[$bi_i]} >= ${bi_B1[$bi_i]} && ${bi_A1[$bi_i]} <= ${bi_B2[$bi_i]}) || \
-		 (${bi_A2[$bi_i]} >= ${bi_B1[$bi_i]} && ${bi_A2[$bi_i]} <= ${bi_B2[$bi_i]})"
-	done
-
-	calculate $3 "${IN[0]} && ${IN[1]}"
-}
-
-### Fetch destination projection and region
-
-SOURCE_PROJ=$GIS_OPT_SOURCEPROJ
-SOURCE_SCALE=$GIS_OPT_SOURCESCALE
-
-# Take into account those extra pixels we'll be a addin'
-#MAX_COLS=$GIS_OPT_MAXCOLS
-#MAX_ROWS=$GIS_OPT_MAXROWS
-OVERLAP=$GIS_OPT_OVERLAP
-calculate MAX_COLS "$GIS_OPT_MAXCOLS - $OVERLAP"
-calculate MAX_ROWS "$GIS_OPT_MAXROWS - $OVERLAP"
-
-g.message "Getting destination projection"
-
-if [ -z "$GIS_OPT_DESTPROJ" ] ; then
-	DEST_PROJ=`g.proj -j`
-else
-	DEST_PROJ=`echo $GIS_OPT_DESTPROJ`
-fi
-
-	g.message "Getting projection scale"
-
-if [ -z "$GIS_OPT_DESTSCALE" ] ; then
-	eval `g.proj -p | $GREP '^meters' | $SED -e "s/:/=/" -e "s/ //g"`
-	DEST_SCALE=$meters;
-else
-	DEST_SCALE=$GIS_OPT_DESTSCALE
-fi
-
-# Strip newlines from both the source and destination projections
-DEST_PROJ=`echo $DEST_PROJ`
-SOURCE_PROJ=`echo $SOURCE_PROJ`
-
-# Set up the projections
-makeProjector SOURCE_PROJ $SOURCE_SCALE DEST_PROJ $DEST_SCALE SOURCE_TO_DEST
-makeProjector DEST_PROJ $DEST_SCALE SOURCE_PROJ $SOURCE_SCALE DEST_TO_SOURCE
-
-g.message "Getting destination region"
-if [ -z $GIS_OPT_REGION ] ; then
-	eval `g.region -g`
-else
-	eval `g.region region=$GIS_OPT_REGION -g`
-fi
-
-
-DEST_BBOX[0]=$w
-DEST_BBOX[1]=$s
-DEST_BBOX[2]=$e
-DEST_BBOX[3]=$n
-DEST_NSRES=$nsres
-DEST_EWRES=$ewres
-calculate DEST_ROWS "($n - $s) / $nsres"
-calculate DEST_COLS "($e - $w) / $ewres"
-
-message 1 "Rows: $DEST_ROWS, Cols: $DEST_COLS"
-
-### Project the destination region into the source:
-g.message "Projecting destination region into source"
-
-bboxToPoints DEST_BBOX DEST_BBOX_POINTS
-
-message 1 "${DEST_BBOX_POINTS[@]}"
-
-# echo "${DEST_TO_SOURCE[@]}"
-# echo "${DEST_TO_SOURCE[0]}"
-# echo "${DEST_TO_SOURCE[1]}"
-# echo "${DEST_TO_SOURCE[2]}"
-# echo "${DEST_TO_SOURCE[3]}"
-
-# exit 1
-
-projectPoints DEST_BBOX_POINTS DEST_BBOX_SOURCE_POINTS DEST_TO_SOURCE 
-
-message 1 "${DEST_BBOX_SOURCE_POINTS[@]}"
-
-pointsToBbox DEST_BBOX_SOURCE_POINTS SOURCE_BBOX
-
-message 1 "${SOURCE_BBOX[@]}"
-
-g.message "Projecting source bounding box into destination"
-
-bboxToPoints SOURCE_BBOX SOURCE_BBOX_POINTS
-
-message 1 "${SOURCE_BBOX_POINTS[@]}"
-
-projectPoints SOURCE_BBOX_POINTS SOURCE_BBOX_DEST_POINTS SOURCE_TO_DEST 
-
-message 1 "${SOURCE_BBOX_DEST_POINTS[@]}"
-
-calculate X_METRIC "1 / $DEST_EWRES"
-calculate Y_METRIC "1 / $DEST_NSRES"
-
-g.message "Computing length of sides of source bounding box"
-
-sideLengths  SOURCE_BBOX_DEST_POINTS $X_METRIC $Y_METRIC SOURCE_BBOX_DEST_LENGTHS
-
-message 1 "${SOURCE_BBOX_DEST_LENGTHS[@]}"
-
-# Find the skewedness of the two directions.
-# Define it to be greater than one
-# In the direction (x or y) in which the world is least skewed (ie north south in lat long)
-# Divide the world into strips. These strips are as big as possible contrained by max_
-# In the other direction do the same thing.
-# Theres some recomputation of the size of the world that's got to come in here somewhere.
-
-# For now, however, we are going to go ahead and request more data than is necessary.
-# For small regions far from the critical areas of projections this makes very little difference
-# in the amount of data gotten.
-# We can make this efficient for big regions or regions near critical points later.
-
-HORIZONTALS=( ${SOURCE_BBOX_DEST_LENGTHS[1]} ${SOURCE_BBOX_DEST_LENGTHS[3]} )
-VERTICALS=( ${SOURCE_BBOX_DEST_LENGTHS[0]} ${SOURCE_BBOX_DEST_LENGTHS[2]} )
-
-max HORIZONTALS BIGGER[0]
-max VERTICALS BIGGER[1]
-
-MAX[0]=$MAX_COLS
-MAX[1]=$MAX_ROWS
-
-# Compute the number and size of tiles to use in each direction
-# I'm making fairly even sized tiles
-# They differer from each other in height and width only by one cell
-# I'm going to make the numbers all simpler and add this extra cell to
-# every tile.
-
-g.message "Computing tiling"
-
-for i in 0 1 ; do
-	# Make these into integers.
-	# Round up
-	calculate BIGGER[$i] "scale=0; (${BIGGER[$i]} + 1) / 1"
-	calculate TILES[$i] "scale=0; (${BIGGER[$i]} / ${MAX[$i]}) + 1"
-	# BC truncates and doesn't round, which is perfect here:
-	calculate TILE_BASE_SIZE[$i] "scale=0; (${BIGGER[$i]} / ${TILES[$i]})"
-	calculate TILES_EXTRA_1[$i] "scale=0; (${BIGGER[$i]} % ${TILES[$i]})"
-	# This is adding the extra pixel (remainder) to all of the tiles:
-	calculate TILE_SIZE[$i] "scale=0; ${TILE_BASE_SIZE[$i]} + (${TILES_EXTRA_1[$i]} > 0)" 
-	calculate TILESET_SIZE[$i] "scale=0; ${TILE_SIZE[$i]} * ${TILES[$i]}"
-	# Add overlap to tiles (doesn't effect tileset_size
-	calculate TILE_SIZE_OVERLAP[$i] "${TILE_SIZE[$i]} + $OVERLAP"
-done;
-
-g.message "There will be ${TILES[0]} by ${TILES[1]} tiles each ${TILE_SIZE[0]} by ${TILE_SIZE[1]} cells."
-
-ximax=${TILES[0]}
-yimax=${TILES[1]}
-
-MIN_X=${SOURCE_BBOX[0]}
-MIN_Y=${SOURCE_BBOX[1]}
-MAX_X=${SOURCE_BBOX[2]}
-MAX_Y=${SOURCE_BBOX[3]}
-SPAN_X="($MAX_X - $MIN_X)"
-SPAN_Y="($MAX_Y - $MIN_Y)"
-
-
-for ((xi=0;xi<ximax;xi+=1)); do
-	calculate TILE_BBOX[0] "$MIN_X + ( $xi * ${TILE_SIZE[0]} / ${TILESET_SIZE[0]} ) * $SPAN_X"
-	calculate TILE_BBOX[2] "$MIN_X + ( ($xi + 1) * ${TILE_SIZE_OVERLAP[0]} / ${TILESET_SIZE[0]} ) * $SPAN_X"
-	for ((yi=0;yi<yimax;yi+=1)); do
-		calculate TILE_BBOX[1] "$MIN_Y + ( $yi * ${TILE_SIZE[1]} / ${TILESET_SIZE[1]} ) * $SPAN_Y"
-		calculate TILE_BBOX[3] "$MIN_Y + ( ($yi + 1) * ${TILE_SIZE_OVERLAP[1]} / ${TILESET_SIZE[1]} ) * $SPAN_Y"
-		bboxToPoints TILE_BBOX TILE_BBOX_POINTS
-		projectPoints TILE_BBOX_POINTS TILE_DEST_BBOX_POINTS SOURCE_TO_DEST
-		pointsToBbox TILE_DEST_BBOX_POINTS TILE_DEST_BBOX
-		bboxesIntersect TILE_DEST_BBOX DEST_BBOX intersect
-		if [ $intersect ] ; then
-			if [ $GIS_FLAG_W -eq 1 ] ; then
-				echo "bbox=${TILE_BBOX[0]},${TILE_BBOX[1]},${TILE_BBOX[2]},${TILE_BBOX[3]}&width=${TILE_SIZE_OVERLAP[0]}&height=${TILE_SIZE_OVERLAP[1]}"
-			elif [ $GIS_FLAG_G -eq 1 ] ; then
-				echo "w=${TILE_BBOX[0]};s=${TILE_BBOX[1]};e=${TILE_BBOX[2]};n=${TILE_BBOX[3]};cols=${TILE_SIZE_OVERLAP[0]};rows=${TILE_SIZE_OVERLAP[1]};"
-			else
-				echo "${TILE_BBOX[0]}$GIS_OPT_FS${TILE_BBOX[1]}$GIS_OPT_FS${TILE_BBOX[2]}$GIS_OPT_FS${TILE_BBOX[3]}$GIS_OPT_FS${TILE_SIZE_OVERLAP[0]}$GIS_OPT_FS${TILE_SIZE_OVERLAP[1]}"
-			fi
-		fi 
-	done
-done
-

scripts/r.tileset/r.tileset.html

@@ -1,16 +1,32 @@
-<H2>DESCRIPTION</H2>
+<h2>DESCRIPTION</h2>
 
-<EM>r.tileset</EM> generates sets of tiles in another projection that cover a region in this projection with adequate resolution. By default the current region and its resolution are used, the bounds and resolution of another region can be used via the region option.
+<em>r.tileset</em> generates sets of tiles in another projection that
+cover a region in this projection with adequate resolution. By default
+the current region and its resolution are used, the bounds and
+resolution of another region can be used via the region option.
 
-<H2>NOTES</H2>
+<h2>NOTES</h2>
 
-<EM>r.tileset</EM> does not make "optimal" tilings (as few tiles of the largest size less than the maximums). This means that from latitude longitude projection to an appropriate projection for a region, in the degenerate case, it may create tiles demanding up to twice the necessary information. Furthermore, generating a tiling near a divergant point of a source projection, such as the poles of a cylindrical source projections, results in divergence of the tile set.
+<em>r.tileset</em> does not make "optimal" tilings (as few tiles of
+the largest size less than the maximums). This means that from
+latitude longitude projection to an appropriate projection for a
+region, in the degenerate case, it may create tiles demanding up to
+twice the necessary information. Furthermore, generating a tiling near
+a divergant point of a source projection, such as the poles of a
+cylindrical source projections, results in divergence of the tile set.
 
-<p>Not generating "optimal" tilings may have another consequence; the aspect ratio of cells in the destination region will not necessarily carry over to the source region and generated tiles may have cells of strange aspect ratios. This might be a problem for some map request services presenting data in an inappropriate projection or with strict constraints on cell aspect ratio.
+<p>
+Not generating "optimal" tilings may have another consequence; the
+aspect ratio of cells in the destination region will not necessarily
+carry over to the source region and generated tiles may have cells of
+strange aspect ratios. This might be a problem for some map request
+services presenting data in an inappropriate projection or with strict
+constraints on cell aspect ratio.
 
-<H2>OUTPUT FORMAT</H2>
+<h2>OUTPUT FORMAT</h2>
 
-Each tile is listed on a seperate line in the output. The lines are formatted as follows:
+Each tile is listed on a seperate line in the output. The lines are
+formatted as follows:
 
 <dl>
 <dt>
@@ -18,86 +34,70 @@ Each tile is listed on a seperate line in the output. The lines are formatted as
 5|125|45|175|80|100
 </tt></span>
 
-<dd>
-This is the default output format. It is the tile's minimum x coordinate, minimum y coordinate, maximum x coordinate, maximum y coordinate, width in cells, and height in cells seperated by the "|" character. The fields can be seperated by a different character by changing the fs option.
-
-<p>
-
+<dd>This is the default output format. It is the tile's minimum x
+coordinate, minimum y coordinate, maximum x coordinate, maximum y
+coordinate, width in cells, and height in cells seperated by the "|"
+character. The fields can be seperated by a different character by
+changing the fs option.
 <dt>
 <span class="code"><tt>
 w=5;s=125;e=45;n=175;cols=80;rows=100;
 </tt></span>
 
-<dd>
-This is output in a format convinent for setting variables in a shell script.
-
-<p>
-
+<dd>This is output in a format convinent for setting variables in a shell
+script.
 <dt>
 <span class="code"><tt>
 bbox=5,125,45,175&amp;width=80&amp;height=100
 </tt></span>
 
-<dd>
-This is output in a format convinent for requesting data from some http services.
-
-<p>
+<dd>This is output in a format convinent for requesting data from some
+http services.
 </dl>
 
-
-<H2>EXAMPLES</H2>
+<h2>EXAMPLES</j2>
 
 <dl>
 <dt>
 <span class="code"><tt>
 r.tileset sourceproj=+init=epsg:4326 maxrows=1024 maxcols=2048
 </tt></span>
-<dd> Generates tiles in latitude longitude that cover the current region, each tile will be less than 1024 cells high and 2048 cells across. The bounds and sizes of tiles in the output are seperated by |
 
-<p>
+<dd>Generates tiles in latitude longitude that cover the current
+region, each tile will be less than 1024 cells high and 2048 cells
+across. The bounds and sizes of tiles in the output are seperated by |
 
 <dt>
 <span class="code"><tt>
 r.tileset sourceproj=+init=epsg:4326 overlap=2 -w region=ne-rio
 </tt></span>
-<dd>Generates tiles in latitude longitude projection that cover the named region "ne-rio". The tiles will have 2 cells of overlap. The output format will be strings like the bbox requests for WMS servers.
 
-<p>
+<dd>Generates tiles in latitude longitude projection that cover the
+named region "ne-rio". The tiles will have 2 cells of overlap. The
+output format will be strings like the bbox requests for WMS servers.
 
 <dt>
 <span class="code"><tt>
 r.tileset sourceproj=`g.proj -j location=IrishGrid` maxrows=400 maxcols=300 overlap=3 -g
 </tt></span>
-<dd>Generates tiles in the projection of the location "IrishGrid". Each tile will be less than 300x400 cells in size, with 3 cells of overlap in the top and right sides of each tile. The output is in a format where each line is in shell script style. The substitution <code>`g.proj -j location=IrishGrid`</code> will only work in a unix style shell.
 
+<dd>Generates tiles in the projection of the location
+"IrishGrid". Each tile will be less than 300x400 cells in size, with 3
+cells of overlap in the top and right sides of each tile. The output
+is in a format where each line is in shell script style. The
+substitution <code>`g.proj -j location=IrishGrid`</code> will only
+work in a unix style shell.
 </dl>
 
-<H2>REQUIRED PROGRAMS</H2>
-
-<EM>r.tileset</EM> requires the following programs to work:
-
-<dl>
-<dt>cs2cs
-<dd>The coordinate system transformation program from Proj4.
-
-<dt>bc
-<dd>A calculator program
-
-<dt>sed, grep
-<dd>Unix string processing and search programs
-
-</dl>
-
-
-<H2>BUGS</H2>
+<h2>BUGS</h2>
 
 <ul>
-<li><EM>r.tileset</EM> does not know about meridians that "wrap-around" in projections.
-</ul>
-
+<li><en>r.tileset</em> does not know about meridians that
+"wrap-around" in projections.  </ul>
 
-<H2>AUTHORS</H2>
+<h2>AUTHORS</h2>
 
-Cedric Shock
+Cedric Shock<br>
+Updated for GRASS 7 by Martin Landa, CTU in Prague, Czech Republic
 
 <p><i>Last changed: $Date$</i>

scripts/r.tileset/r.tileset.py

@@ -0,0 +1,389 @@
+#!/usr/bin/env python
+
+############################################################################
+#
+# MODULE:       r.tileset
+#
+# AUTHOR(S):    Cedric Shock
+#               Updated for GRASS7 by Martin Landa, 2009
+#
+# PURPOSE:      To produce tilings of regions in other projections.
+#
+# COPYRIGHT:    (C) 2006-2009 by Cedric Shoc, Martin Landa, and 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.
+#
+#############################################################################
+
+# Bugs:
+#  Does not know about meridians in projections. However, unlike the usual
+#   hack used to find meridians, this code is perfectly happy with arbitrary
+#   rotations and flips
+#  The following are planned to be fixed in a future version, if it turns out
+#   to be necessary for someone:
+#   Does not generate optimal tilings. This means that between an appropriate
+#    projection for a region and latitude longitude projection, in the
+#    degenerate case, it may create tiles demanding up to twice the necessary
+#    information. Requesting data from cylindrical projections near their poles
+#    results in divergence. You really don't want to use source data from
+#    someone who put it in a cylindrical projection near a pole, do you?
+#   Not generating "optimal" tilings has another side effect; the sanity
+#    of the destination region will not carry over to generating tiles of
+#    realistic aspect ratio. This might be a problem for some WMS servers
+#    presenting data in a highly inappropriate projection. Do you really
+#    want their data?
+
+#%Module
+#% description: Produces tilings of the source projection for use in the destination region and projection.
+#% keywords: raster, tiling
+#%End
+#%flag
+#% key: g
+#% description: Produces shell script output
+#%end
+#%flag
+#% key: w
+#% description: Produces web map server query string output
+#%end
+#%option
+#% key: region
+#% type: string
+#% description: Name of region to use instead of current region for bounds and resolution
+#%end
+#%option
+#% key: sourceproj
+#% type: string
+#% description: Source projection
+#% required : yes
+#%end
+#%option
+#% key: sourcescale
+#% type: string
+#% description: Conversion factor from units to meters in source projection
+#% answer : 1
+#%end
+#%option
+#% key: destproj
+#% type: string
+#% description: Destination projection, defaults to this location's projection
+#% required : no
+#%end
+#%option
+#% key: destscale
+#% type: string
+#% description: Conversion factor from units to meters in source projection
+#% required : no
+#%end
+#%option
+#% key: maxcols
+#% type: integer
+#% description: Maximum number of columns for a tile in the source projection
+#% answer: 1024
+#%end
+#%option
+#% key: maxrows
+#% type: integer
+#% description: Maximum number of rows for a tile in the source projection
+#% answer: 1024
+#%end
+#%option
+#% key: overlap
+#% type: integer
+#% description: Number of cells tiles should overlap in each direction
+#% answer: 0
+#%end
+#%option
+#% key: fs
+#% type: string
+#% description: Output field separator
+#% answer: |
+#%end
+#%option
+#% key: v
+#% type: integer
+#% description: Verbosity level
+#% answer: 0
+#%end
+
+# Data structures used in this program:
+# A bounding box:
+# 0 -> left, 1-> bottom, 2->right, 3-> top
+# A border:
+# An array of points indexed by 0 for "x" and 4 for "y" + by number 0, 1, 2, and 3
+# A reprojector [0] is name of source projection, [1] is name of destination
+# A projection - [0] is proj.4 text, [1] is scale
+
+import sys
+import subprocess
+import tempfile
+import math
+
+import grass
+
+def bboxToPoints(bbox):
+    """Make points that are the corners of a bounding box"""
+    points = []
+    points.append((bbox['w'], bbox['s']))
+    points.append((bbox['w'], bbox['n']))
+    points.append((bbox['e'], bbox['n']))
+    points.append((bbox['e'], bbox['s']))
+    
+    return points
+
+def pointsToBbox(points):
+    bbox = {}
+    min_x = min_y = max_x = max_y = None
+    for point in points:
+        if not min_x:
+            min_x = max_x = point[0]
+        if not min_y:
+            min_y = max_y = point[1]
+        
+        if min_x > point[0]:
+            min_x = point[0]
+        if max_x < point[0]:
+            max_x = point[0]
+        if min_y > point[1]:
+            min_y = point[1]
+        if max_y < point[1]:
+            max_y = point[1]
+    
+    bbox['n'] = max_y
+    bbox['s'] = min_y
+    bbox['w'] = min_x
+    bbox['e'] = max_x
+    
+    return bbox
+
+def project(file, source, dest):
+    """Projects point (x, y) using projector"""
+    points = []
+    ret = grass.read_command('m.proj',
+                             quiet = True,
+                             flags = 'd',
+                             proj_in = source['proj'],
+                             proj_out = dest['proj'],
+                             fs = ';,;',
+                             input = file)
+    
+    if not ret:
+        grass.fatal(cs2cs + ' failed')
+    
+    for line in ret.splitlines():
+        p_x2, p_y2, p_z2 = map(float, line.split(';'))
+        points.append((p_x2 / dest['scale'], p_y2 / dest['scale']))
+    
+    return points
+    
+def projectPoints(points, source, dest):
+    """Projects a list of points"""
+    dest_points = []
+    
+    input = tempfile.NamedTemporaryFile(mode="wt")
+    for point in points:
+        input.file.write('%f;%f\n' % \
+                             (point[0] * source['scale'],
+                              point[1] * source['scale']))
+    input.file.flush()
+    
+    dest_points = project(input.name, source, dest)
+    
+    return dest_points
+
+def sideLengths(points, xmetric, ymetric):
+    """Find the length of sides of a set of points from one to the next"""
+    ret = []
+    for i in range(len(points)):
+        x1, y1 = points[i]
+        j = i + 1
+        if j >= len(points):
+            j = 0
+        sl_x = (points[j][0] - points[i][0]) * xmetric
+        sl_y = (points[j][1] - points[i][1]) * ymetric
+        sl_d = math.sqrt(sl_x * sl_x + sl_y * sl_y)
+        ret.append(sl_d)
+    
+    return { 'x' : (ret[1], ret[3]),
+             'y' : (ret[0], ret[2]) }
+
+def bboxesIntersect(bbox_1, bbox_2):
+    """Determine if two bounding boxes intersect"""
+    bi_a1 = (bbox_1['w'], bbox_1['s'])
+    bi_a2 = (bbox_1['e'], bbox_1['n'])
+    bi_b1 = (bbox_2['w'], bbox_2['s'])
+    bi_b2 = (bbox_2['e'], bbox_2['n'])
+    cin = [False, False]
+    for i in (0, 1):
+        if (bi_a1[i] <= bi_b1[i] and bi_a2[i] >= bi_b2[i]) or \
+               (bi_a1[i] <= bi_b1[i] and bi_a2[i] >= bi_b2[i]) or \
+               (bi_a1[i] <= bi_b1[i] and bi_a1[i] >= bi_b2[i]) or \
+               (bi_a2[i] <= bi_b1[i] and bi_a2[i] >= bi_b2[i]):
+            cin[i] = True
+    
+    if cin[0] and cin[1]:
+        return True
+    
+    return False
+
+def main():
+    # Take into account those extra pixels we'll be a addin'
+    max_cols = int(options['maxcols']) - int(options['overlap'])
+    max_rows = int(options['maxrows']) - int(options['overlap'])
+    
+    # destination projection
+    if not options['destproj']:
+        dest_proj = grass.read_command('g.proj',
+                                       quiet = True,
+                                       flags = 'jf').rstrip('\n')
+        if not dest_proj:
+            grass.fatal('g.proj failed')
+    else:
+        dest_proj = options['destproj']
+    grass.debug("Getting destination projection -> '%s'" % dest_proj)
+        
+    # projection scale
+    if not options['destscale']:
+        ret = grass.parse_command('g.proj',
+                                 quiet = True,
+                                 flags = 'j')
+        if not ret:
+            grass.fatal('g.proj failed')
+        
+        dest_scale = ret['+to_meter'].strip()
+    else:
+        dest_scale = options['destscale']
+    grass.debug('Getting destination projection scale -> %s' % dest_scale)
+    
+    # set up the projections
+    srs_source = { 'proj' : options['sourceproj'],
+                   'scale' : float(options['sourcescale']) }
+    srs_dest   = { 'proj' : dest_proj,
+                   'scale' : float(dest_scale) }   
+    
+    if options['region']:
+        grass.run_command('g.region',
+                          quiet = True,
+                          region = options['region'])
+    dest_bbox = grass.region()
+    grass.debug('Getting destination region')
+    
+    # project the destination region into the source:
+    grass.verbose('Projecting destination region into source...')
+    dest_bbox_points = bboxToPoints(dest_bbox)
+    
+    dest_bbox_source_points = projectPoints(dest_bbox_points,
+                                            source = srs_dest,
+                                            dest = srs_source)
+    
+    source_bbox = pointsToBbox(dest_bbox_source_points)
+    
+    grass.verbose('Projecting source bounding box into destination...')
+    
+    source_bbox_points = bboxToPoints(source_bbox)
+    
+    source_bbox_dest_points = projectPoints(source_bbox_points,
+                                            source = srs_source,
+                                            dest = srs_dest)
+    
+    x_metric = 1 / dest_bbox['ewres']
+    y_metric = 1 / dest_bbox['nsres']
+    
+    grass.verbose('Computing length of sides of source bounding box...')
+    
+    source_bbox_dest_lengths = sideLengths(source_bbox_dest_points,
+                                           x_metric, y_metric)
+    
+    # Find the skewedness of the two directions.
+    # Define it to be greater than one
+    # In the direction (x or y) in which the world is least skewed (ie north south in lat long)
+    # Divide the world into strips. These strips are as big as possible contrained by max_
+    # In the other direction do the same thing.
+    # Theres some recomputation of the size of the world that's got to come in here somewhere.
+    
+    # For now, however, we are going to go ahead and request more data than is necessary.
+    # For small regions far from the critical areas of projections this makes very little difference
+    # in the amount of data gotten.
+    # We can make this efficient for big regions or regions near critical points later.
+    
+    bigger = []
+    bigger.append(max(source_bbox_dest_lengths['x']))
+    bigger.append(max(source_bbox_dest_lengths['y']))
+    maxdim = (max_cols, max_rows)
+    
+    # Compute the number and size of tiles to use in each direction
+    # I'm making fairly even sized tiles
+    # They differer from each other in height and width only by one cell
+    # I'm going to make the numbers all simpler and add this extra cell to
+    # every tile.
+    
+    grass.message('Computing tiling...')
+    tiles = [-1, -1]
+    tile_base_size = [-1, -1]
+    tiles_extra_1 = [-1, -1]
+    tile_size = [-1, -1]
+    tileset_size = [-1, -1]
+    tile_size_overlap = [-1, -1]
+    for i in range(len(bigger)):
+	# make these into integers.
+	# round up
+	bigger[i] = int(bigger[i] + 1)
+	tiles[i] = int((bigger[i] / maxdim[i]) + 1)
+        tile_size[i] = tile_base_size[i] = int(bigger[i] / tiles[i])
+	tiles_extra_1[i] = int(bigger[i] % tiles[i])
+	# This is adding the extra pixel (remainder) to all of the tiles:
+        if tiles_extra_1[i] > 0:
+            tile_size[i] = tile_base_size[i] + 1
+        tileset_size[i] = int(tile_size[i] * tiles[i])
+	# Add overlap to tiles (doesn't effect tileset_size
+        tile_size_overlap[i] = tile_size[i] + int(options['overlap'])
+    
+    grass.verbose("There will be %d by %d tiles each %d by %d cells" % \
+                  (tiles[0], tiles[1], tile_size[0], tile_size[1]))
+    
+    ximax = tiles[0]
+    yimax = tiles[1]
+    
+    min_x = source_bbox['w']
+    min_y = source_bbox['s']
+    max_x = source_bbox['e']
+    max_y = source_bbox['n']
+    span_x = (max_x - min_x)
+    span_y = (max_y - min_y)
+    
+    xi = 0
+    tile_bbox = { 'w' : -1, 's': -1, 'e' : -1, 'n' : -1 }
+    while xi < ximax:
+        tile_bbox['w'] = min_x + (xi * tile_size[0] / tileset_size[0]) * span_x
+        tile_bbox['e'] = min_x + ((xi + 1) * tile_size_overlap[0] / tileset_size[0]) * span_x
+	yi = 0
+        while yi < yimax:
+            tile_bbox['s'] = min_y + (yi * tile_size[1] / tileset_size[1]) * span_y
+            tile_bbox['n'] = min_y + ((yi + 1) * tile_size_overlap[1] / tileset_size[1]) * span_y
+            tile_bbox_points = bboxToPoints(tile_bbox)
+            tile_dest_bbox_points = projectPoints(tile_bbox_points,
+                                                  source = srs_source,
+                                                  dest = srs_dest)
+            tile_dest_bbox = pointsToBbox(tile_dest_bbox_points)
+            if bboxesIntersect(tile_dest_bbox, dest_bbox):
+                if flags['w']:
+                    print "bbox=%s,%s,%s,%s&width=%s&height=%s" % \
+                          (tile_bbox['w'], tile_bbox['s'], tile_bbox['e'], tile_bbox['n'],
+                           tile_size_overlap[0], tile_size_overlap[1])
+                elif flags['g']:
+                    print "w=%s;s=%s;e=%s;n=%s;cols=%s;rows=%s" % \
+                          (tile_bbox['w'], tile_bbox['s'], tile_bbox['e'], tile_bbox['n'],
+                           tile_size_overlap[0], tile_size_overlap[1])
+                else:
+                    fs = options['fs']
+                    print "%s%s%s%s%s%s%s%s%s%s%s" % \
+                          (tile_bbox['w'], fs, tile_bbox['s'], fs, tile_bbox['e'], fs, tile_bbox['n'], fs,
+                           tile_size_overlap[0], fs, tile_size_overlap[1])
+            yi += 1
+        xi += 1
+    
+if __name__ == "__main__":
+    cs2cs = 'cs2cs'
+    options, flags = grass.parser()
+    sys.exit(main())