public
Last active

Slippy map tiles to Static Map

  • Download Gist
math-working.tex
TeX
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96
% 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
Perl
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163
#!/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;

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.