Instantly share code, notes, and snippets.

# andrewharvey/math-working.tex Created Sep 17, 2012

What would you like to do?
Slippy map tiles to Static Map
 % 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}
 #!/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 # 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;