Slippy map tiles to Static Map

math-working.tex

% To the extent possible under law, the person who associated CC0
% with this work has waived all copyright and related or neighboring
% rights to this work.
% http://creativecommons.org/publicdomain/zero/1.0/

\documentclass[a4paper,10pt]{article}
\usepackage{verbatim} 
\usepackage{amsmath}
\usepackage{amssymb}
\setlength\parindent{0mm}
\usepackage{fullpage}
\usepackage[all]{xy}
\usepackage{graphicx}

\usepackage[pdftex,
            pdfauthor={Andrew Harvey},
            pdftitle={OSM Tiles to Static Map Working},
            pdfsubject={},
            pdfkeywords={},
            pdfstartview=FitB]{hyperref}

\begin{document}
This document describes the working of the accompying code base for creating a static map from OSM slippy map tiles.

\begin{center}
 \includegraphics{math-working-diagram.eps}
 % diagram.pdf: 176x213 pixel, 72dpi, 6.21x7.51 cm, bb=0 0 176 213
\end{center}

The blue outline is the extent of the entire world map. In this example we are at zoom 1 since there are a total of 4 tiles.

The dotted rectangle shows the extent of the static map requested within the world map.

Unless otherwise noted, $O$ shall be considered the origin and all units are in pixels.

$t$ is the width of one tile which is usually 256. We assume all tiles are square.

\section{Inputs}
\label{inputs:}

The paramaters given for the static map are,

\begin{itemize}
\item $I$ center point of static map in geodetic WGS 84
\item $w$ width of the static map in pixels
\item $h$ height of the static map in pixels
\item $z$ zoom level of tiles to use
\end{itemize}

\section{Outputs}
\label{outputs:}
In order to construct the final static map from the tiles we need to determine,
\begin{itemize}
\item $A'$, The coordinates of $A$ in pixels relative to the top left corner of the tile which $A$ falls within.
\item $T_A$, The slippy map tile reference of the tile which $A$ falls within.
\item $T_B$, The slippy map tile reference of the tile which $B$ falls within.
\end{itemize}

Once we have these we can simply retieve and position all the tiles formed by the area $T_A : T_B$ and offset it by $-A'$.

\section{Working}
\label{working:}
We shall first take $I$ and project it into the Popular Visualisation Mercator projection and denote it $C$.

$$C = \binom{\mbox{deg2rad}\left( I_x \right) }{asinh \left( \mbox{deg2rad}\left( I_y \right)\right)  }$$

We shall then transform $C$ into a ratio of its position within the entire world map and denote it $\omega$,

$$\omega = \begin{pmatrix}
1 + \dfrac{\left( \frac{C_x}{\pi}\right) }{2}\\
1 - \dfrac{\left( \frac{C_y}{\pi}\right) }{2}
\end{pmatrix}$$

Now,

$$l = t \times 2^z$$

So $I$ in absolute pixel coordinates, $P$ can be defined as,

$$P = l \times \omega$$

Hence,

$$A = P - \dfrac{w,h}{2}$$
$$B = P + \dfrac{w,h}{2}$$

We can then find the tile references of the tile containing A and B,

$$T_A = \mbox{floor} \left( \dfrac{A}{t} \right)$$
$$T_B = \mbox{floor} \left( \dfrac{B}{t} \right)$$

Finally we need $A'$ which is,

$$A' = A - \left( T_A \times 256 \right)$$

\end{document}

tiles2staticmap.pl

#!/usr/bin/perl -w

# This script will construct a "static map" from a tile server.
#
# Given a center point, zoom level and image dimensions it will retrieve
# all the tiles it needs to construct a single image of the area and will
# go ahead and construct that single static map image.
#
# Currently functions aren't abstracted, ideally the functional parts of
# this script would be abstracted away so it can easily be adapted to
# work as either a server side script (CGI), client side script (local
# program), and either retrive tiles from a HTTP tile server or via the
# local filesystem.

# Author: Andrew Harvey <andrew.harvey4@gmail.com>
# License: CC0 http://creativecommons.org/publicdomain/zero/1.0/
#
# To the extent possible under law, the person who associated CC0
# with this work has waived all copyright and related or neighboring
# rights to this work.
# http://creativecommons.org/publicdomain/zero/1.0/

use strict;

use Getopt::Long;
use POSIX; # for floor function
use Math::Trig qw(pi deg2rad asinh);

use GD;
use LWP;

use Geo::Proj4;

use Log::Log4perl qw(:easy);
Log::Log4perl->easy_init($INFO); #DEBUG, INFO, WARN, ERROR, FATAL

my $lat;
my $lon;
my $zoomLevel;

my $mapWidth;
my $mapHeight;

my $tileBase;
my $outputFile;

my $opts = GetOptions(
  "lat=f" => \$lat,
  "lon=f" => \$lon,
  "zoom=i" => \$zoomLevel,
  "width=i" => \$mapWidth,
  "height=i" => \$mapHeight,
  "tileBase=s" => \$tileBase,
  "output=s" => \$outputFile );

INFO "map request: lat $lat, lon $lon, zoom $zoomLevel of size ${mapWidth}x${mapHeight}px";

# width and height of tiles in pixels
my $tileSizeInPixels = 256;

# number of vectical or horizontal tiles for the current zoom level
my $numTilesAlongSingleAxis = 2 ** $zoomLevel;

# number of pixels in vertical and horizontal direction for the whole world
my $worldSizeInPixels = $tileSizeInPixels * $numTilesAlongSingleAxis;

DEBUG "At zoom $zoomLevel there are ${numTilesAlongSingleAxis}x${numTilesAlongSingleAxis} tiles";

# project to the Popular Visualisation Mercator projection
my $toPopularVisMercator = Geo::Proj4->new ('+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs +over');
my ($centerInMercX, $centerInMercY) = $toPopularVisMercator->forward($lat, $lon);

my ($projExtentX, $projExtentY) = $toPopularVisMercator->forward(-85, 180);

$projExtentY = -$projExtentX; # FIXME why is this really needed?

DEBUG "The requested map center in Popular Visualisation Mercator projection is easting $centerInMercX of $projExtentX, northing $centerInMercY of $projExtentY";

# transform range of x and y to 0-1 and shift origin to top left corner
my $centerRatioX = (1 + ($centerInMercX / $projExtentX)) / 2;
my $centerRatioY = (1 - ($centerInMercY / -$projExtentY)) / 2;

DEBUG "Center has a ratio of the entire world that is $centerRatioX, $centerRatioY";

# get absolute pixel of centre point
my $centerAbsoluteX = $centerRatioX * $worldSizeInPixels;
my $centerAbsoluteY = $centerRatioY * $worldSizeInPixels;

DEBUG "Center has a worldwide absolute pixel that is $centerAbsoluteX of $worldSizeInPixels, $centerAbsoluteY of $worldSizeInPixels";

my $topLeftPixelX = $centerAbsoluteX - ($mapWidth / 2);
my $topLeftPixelY = $centerAbsoluteY - ($mapHeight / 2);

my $bottomRightPixelX = $centerAbsoluteX + ($mapWidth / 2) - 1;
my $bottomRightPixelY = $centerAbsoluteY + ($mapHeight / 2) - 1;

my $tileRefAX = floor($topLeftPixelX / $tileSizeInPixels);
my $tileRefAY = floor($topLeftPixelY / $tileSizeInPixels);

my $tileRefBX = floor($bottomRightPixelX / $tileSizeInPixels);
my $tileRefBY = floor($bottomRightPixelY / $tileSizeInPixels);

DEBUG "Tile reference of A is $tileRefAX, $tileRefAY; B is $tileRefBX, $tileRefBY";

my $offsetX = $topLeftPixelX - ($tileRefAX * $tileSizeInPixels);
my $offsetY = $topLeftPixelY - ($tileRefAY * $tileSizeInPixels);

DEBUG "Will need to retrieve tiles ${tileRefAX}/${tileRefAY} through ${tileRefBX}/${tileRefBY}";

INFO "Will retrieve a total of " . 
  ($tileRefBX - $tileRefAX + 1) . 
  "x" . 
  ($tileRefBY - $tileRefAY + 1) . 
  " = " . 
  (($tileRefBX - $tileRefAX + 1) * ($tileRefBY - $tileRefAY + 1)) . 
  " tiles.";

DEBUG "We will then need too offset this tile set by -$offsetX, -$offsetY.";

# now construct the final static map from the tiles
# tell GD to always use 24bit color
GD::Image->trueColor(1);

my $img = GD::Image->new($mapWidth, $mapHeight);

# handle transparency properly
$img->alphaBlending(0);
$img->saveAlpha(1);

my $ua = LWP::UserAgent->new();

# get all the tiles we need to cover the area of this static map
for (my $tx = $tileRefAX; $tx <= $tileRefBX; $tx++) {
  for (my $ty = $tileRefAY; $ty <= $tileRefBY; $ty++) {
    if (($tx >= 0) && ($ty >= 0) && ($tx < $numTilesAlongSingleAxis) && ($ty < $numTilesAlongSingleAxis)) {
      my $tileURL = $tileBase;
      $tileURL =~ s/{x}/$tx/g;
      $tileURL =~ s/{y}/$ty/g;
      $tileURL =~ s/{z}/$zoomLevel/g;

      INFO "GET $tileURL";

      my $getResponse = $ua->get($tileURL);
      die "GET $tileURL failed with " . $getResponse->status_line . "\n" unless ($getResponse->is_success);
      my $tile = GD::Image->new($getResponse->content);
      die "Unexpected tile size of " . $tile->width . "x" . $tile->height . "\n" unless (($tile->width == $tileSizeInPixels) && ($tile->height == $tileSizeInPixels));

      my $dx = (($tx - $tileRefAX) * $tileSizeInPixels) - $offsetX;
      my $dy = (($ty - $tileRefAY) * $tileSizeInPixels) - $offsetY;

      DEBUG "Pasting tile $tx/$ty at ${dx}px, ${dy}px";

      $img->copy($tile, $dx, $dy, 0, 0, $tileSizeInPixels, $tileSizeInPixels);
    } # else tile is outside valid range so don't fill it in in the final image
  }
}


# write out the static map
binmode STDOUT;
open my $output_fh, ">$outputFile";
print $output_fh $img->png();
close $output_fh;