andrewharvey / math-working.tex

Created September 17, 2012


    
    View math-working.tex 
  
  
      math-working.tex
    
    
        Raw
      
    
    File suppressed. Click to show.
  
  
            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}

          
    View math-working.tex 
  
  
      tiles2staticmap.pl
    
    
        Raw
      
    
    File suppressed. Click to show.
  
  
            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
            164
          
          
            #!/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;