Slippy map tilenames

From OpenStreetMap Wiki
Jump to: navigation, search

This article describes the file naming conventions for the Slippy Map application.

  • Tiles are 256 × 256 pixel PNG files
  • Each zoom level is a directory, each column is a subdirectory, and each tile in that column is a file
  • Filename(url) format is /zoom/x/y.png

The slippy map expects tiles to be served up at URLs following this scheme, so all tile server URLs look pretty similar.

Tile servers

The first part of the URL specifies the tile server, and perhaps other parameters which might influence the style.

Generally several subdomains (server names) are provided to get around browser limitations on the number of simultaneous HTTP connections to each host. Browser-based applications can thus request multiple tiles from multiple subdomains faster than from one subdomain. For example, OSM, OpenCycleMap and CloudMade servers have three subdomains (a.tile, b.tile, c.tile), MapQuest has four (otile1, otile2, otile3, otile4), all pointing to the same CDN.

That all comes before the /zoom/x/y.png tail.

Here are some examples:

Name URL template zoomlevels
OSM 'standard' style http://[abc].tile.openstreetmap.org/zoom/x/y.png 0-19
OpenCycleMap http://[abc].tile.opencyclemap.org/cycle/zoom/x/y.png 0-18
OpenCycleMap Transport (experimental) http://[abc].tile2.opencyclemap.org/transport/zoom/x/y.png 0-18
MapQuest http://otile[1234].mqcdn.com/tiles/1.0.0/osm/zoom/x/y.jpg 0-19
MapQuest Open Aerial http://otile[1234].mqcdn.com/tiles/1.0.0/sat/zoom/x/y.jpg 0-11 globally, 12+ in the U.S.
Migurski's Terrain http://tile.stamen.com/terrain-background/zoom/x/y.jpg 4-18, US-only (for now)

Further tilesets are available from various '3rd party' sources.

Zoom levels

The zoom parameter is an integer between 0 (zoomed out) and 18 (zoomed in). 18 is normally the maximum, but some tile servers might go beyond that.

zoom level tile coverage number of tiles tile size in degrees
0 1 tile covers whole world 1 tile 360° x 170.1022°
1 2 × 2 tiles 4 tiles 180° x 85.0511°
2 4 × 4 tiles 16 tiles 90° x 42.5256°
n 2n × 2n tiles 22n tiles 360/2n° x 170.1022/2n°
12 4096 x 4096 tiles 16 777 216 0.0879° x 0.0415°
16 232 = 4 294 967 296 tiles
17 17 179 869 184 tiles
18 68 719 476 736 tiles
19 Maximum zoom for Mapnik layer 274 877 906 944 tiles

See Zoom levels for more details

X and Y

  • X goes from 0 (left edge is 180 °W) to 2zoom − 1 (right edge is 180 °E)
  • Y goes from 0 (top edge is 85.0511 °N) to 2zoom − 1 (bottom edge is 85.0511 °S) in a Mercator projection

For the curious, the num 85.0511 is the result of arctan(sinh(π)). By using this bound, the entire map becomes a (very large) square. See also the Osmarender bug.

Derivation of tile names

  • Reproject the coordinates to the Mercator projection (from EPSG:4326 to EPSG:3857):
    • x = lon
    • y = arsinh(tan(lat)) = log[tan(lat) + sec(lat)]
    (lat and lon are in radians)
  • Transform range of x and y to 0 – 1 and shift origin to top left corner:
    • x = [1 + (x / π)] / 2
    • y = [1 − (y / π)] / 2
  • Calculate the number of tiles across the map, n, using 2zoom
  • Multiply x and y by n. Round results down to give tilex and tiley.

Implementations

Pseudo-code

For those who like pseudo-code, here's some hints:

sec = 1/cos
arsinh(x) = log(x + (x^2 + 1)^0.5)
sec^2(x) = tan^2(x) + 1
→ arsinh(tan(x)) = log(tan(x) + sec(x))

Please note that "log" represents the natural logarithm (also known as ln or loge), not decimal logarithm (log10), as used on some calculators.

Lon./lat. to tile numbers

n = 2 ^ zoom
xtile = n * ((lon_deg + 180) / 360)
ytile = n * (1 - (log(tan(lat_rad) + sec(lat_rad)) / π)) / 2

Tile numbers to lon./lat.

n = 2 ^ zoom
lon_deg = xtile / n * 360.0 - 180.0
lat_rad = arctan(sinh(π * (1 - 2 * ytile / n)))
lat_deg = lat_rad * 180.0 / π

Mathematics

Idem with mathematic signs (lat and lon in degrees):

Latlon to tile.pngTile to latlon.png

Python

Lon./lat. to tile numbers

import math
def deg2num(lat_deg, lon_deg, zoom):
  lat_rad = math.radians(lat_deg)
  n = 2.0 ** zoom
  xtile = int((lon_deg + 180.0) / 360.0 * n)
  ytile = int((1.0 - math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * n)
  return (xtile, ytile)

Tile numbers to lon./lat.

import math
def num2deg(xtile, ytile, zoom):
  n = 2.0 ** zoom
  lon_deg = xtile / n * 360.0 - 180.0
  lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
  lat_deg = math.degrees(lat_rad)
  return (lat_deg, lon_deg)

This returns the NW-corner of the square. Use the function with xtile+1 and/or ytile+1 to get the other corners. With xtile+0.5 & ytile+0.5 it will return the center of the tile.

Another python implementation

Ruby

Lon./lat. to tile numbers

def get_tile_number(lat_deg, lng_deg, zoom)
  lat_rad = lat_deg/180 * Math::PI
  n = 2.0 ** zoom
  x = ((lng_deg + 180.0) / 360.0 * n).to_i
  y = ((1.0 - Math::log(Math::tan(lat_rad) + (1 / Math::cos(lat_rad))) / Math::PI) / 2.0 * n).to_i
 
  {:x => x, :y =>y}
end

Tile numbers to lon./lat.

def get_lat_lng_for_number(xtile, ytile, zoom)
  n = 2.0 ** zoom
  lon_deg = xtile / n * 360.0 - 180.0
  lat_rad = Math::atan(Math::sinh(Math::PI * (1 - 2 * ytile / n)))
  lat_deg = 180.0 * (lat_rad / Math::PI)
  {:lat_deg => lat_deg, :lng_deg => lon_deg}
end

Same as the Python implementation above, this returns the NW-corner of the square. Use the function with xtile+1 and/or ytile+1 to get the other corners. With xtile+0.5 & ytile+0.5 it will return the center of the tile.

Perl

Lon./lat. to tile numbers

use Math::Trig;
sub getTileNumber {
  my ($lat,$lon,$zoom) = @_;
  my $xtile = int( ($lon+180)/360 * 2**$zoom ) ;
  my $ytile = int( (1 - log(tan(deg2rad($lat)) + sec(deg2rad($lat)))/pi)/2 * 2**$zoom ) ;
  return ($xtile, $ytile);
}

Tile numbers to lon./lat.

use Math::Trig;
sub Project {
  my ($X,$Y, $Zoom) = @_;
  my $Unit = 1 / (2 ** $Zoom);
  my $relY1 = $Y * $Unit;
  my $relY2 = $relY1 + $Unit;
 
  # note: $LimitY = ProjectF(degrees(atan(sinh(pi)))) = log(sinh(pi)+cosh(pi)) = pi
  # note: degrees(atan(sinh(pi))) = 85.051128..
  #my $LimitY = ProjectF(85.0511);
 
  # so stay simple and more accurate
  my $LimitY = pi;
  my $RangeY = 2 * $LimitY;
  $relY1 = $LimitY - $RangeY * $relY1;
  $relY2 = $LimitY - $RangeY * $relY2;
  my $Lat1 = ProjectMercToLat($relY1);
  my $Lat2 = ProjectMercToLat($relY2);
  $Unit = 360 / (2 ** $Zoom);
  my $Long1 = -180 + $X * $Unit;
  return ($Lat2, $Long1, $Lat1, $Long1 + $Unit); # S,W,N,E
}
sub ProjectMercToLat($){
  my $MercY = shift;
  return rad2deg(atan(sinh($MercY)));
}
sub ProjectF
{
  my $Lat = shift;
  $Lat = deg2rad($Lat);
  my $Y = log(tan($Lat) + sec($Lat));
  return $Y;
}

Lon./lat. to bbox

use Math::Trig;
 
sub getTileNumber {
    my ($lat,$lon,$zoom) = @_;
    my $xtile = int( ($lon+180)/360 * 2**$zoom ) ;
    my $ytile = int( (1 - log(tan(deg2rad($lat)) + sec(deg2rad($lat)))/pi)/2 * 2**$zoom ) ;
    return ($xtile, $ytile);
}
 
sub getLonLat {
	my ($xtile, $ytile, $zoom) = @_;
	my $n = 2 ** $zoom;
	my $lon_deg = $xtile / $n * 360.0 - 180.0;
	my $lat_deg = rad2deg(atan(sinh(pi * (1 - 2 * $ytile / $n))));
	return ($lon_deg, $lat_deg);
}
 
# convert from permalink OSM format like:
# http://www.openstreetmap.org/?lat=43.731049999999996&lon=15.79375&zoom=13&layers=M
# to OSM "Export" iframe embedded bbox format like:
# http://www.openstreetmap.org/export/embed.html?bbox=15.7444,43.708,15.8431,43.7541&layer=mapnik
 
sub LonLat_to_bbox {
	my ($lat, $lon, $zoom) = @_;
 
	my $width = 425; my $height = 350;	# note: must modify this to match your embed map width/height in pixels
	my $tile_size = 256;
 
	my ($xtile, $ytile) = getTileNumber ($lat, $lon, $zoom);
 
	my $xtile_s = ($xtile * $tile_size - $width/2) / $tile_size;
	my $ytile_s = ($ytile * $tile_size - $height/2) / $tile_size;
	my $xtile_e = ($xtile * $tile_size + $width/2) / $tile_size;
	my $ytile_e = ($ytile * $tile_size + $height/2) / $tile_size;
 
	my ($lon_s, $lat_s) = getLonLat($xtile_s, $ytile_s, $zoom);
	my ($lon_e, $lat_e) = getLonLat($xtile_e, $ytile_e, $zoom);
 
	my $bbox = "$lon_s,$lat_s,$lon_e,$lat_e";
	return $bbox;
}

PHP

Lon./lat. to tile numbers

$xtile = floor((($lon + 180) / 360) * pow(2, $zoom));
$ytile = floor((1 - log(tan(deg2rad($lat)) + 1 / cos(deg2rad($lat))) / pi()) /2 * pow(2, $zoom));

Tile numbers to lon./lat.

$n = pow(2, $zoom);
$lon_deg = $xtile / $n * 360.0 - 180.0;
$lat_deg = rad2deg(atan(sinh(pi() * (1 - 2 * $ytile / $n))));

ECMAScript (JavaScript/ActionScript, etc.)

 function long2tile(lon,zoom) { return (Math.floor((lon+180)/360*Math.pow(2,zoom))); }
 function lat2tile(lat,zoom)  { return (Math.floor((1-Math.log(Math.tan(lat*Math.PI/180) + 1/Math.cos(lat*Math.PI/180))/Math.PI)/2 *Math.pow(2,zoom))); }

Inverse process:

 function tile2long(x,z) {
  return (x/Math.pow(2,z)*360-180);
 }
 function tile2lat(y,z) {
  var n=Math.PI-2*Math.PI*y/Math.pow(2,z);
  return (180/Math.PI*Math.atan(0.5*(Math.exp(n)-Math.exp(-n))));
 }

Example: Tilesname WebCalc V1.0

C/C++

int long2tilex(double lon, int z) 
{ 
	return (int)(floor((lon + 180.0) / 360.0 * pow(2.0, z))); 
}
 
int lat2tiley(double lat, int z)
{ 
	return (int)(floor((1.0 - log( tan(lat * M_PI/180.0) + 1.0 / cos(lat * M_PI/180.0)) / M_PI) / 2.0 * pow(2.0, z))); 
}
 
double tilex2long(int x, int z) 
{
	return x / pow(2.0, z) * 360.0 - 180;
}
 
double tiley2lat(int y, int z) 
{
	double n = M_PI - 2.0 * M_PI * y / pow(2.0, z);
	return 180.0 / M_PI * atan(0.5 * (exp(n) - exp(-n)));
}

Go

Example : [1] Doc : [2]

import (
	"math"
)
 
type Tile struct {
	Z    int
	X    int
	Y    int
	Lat  float64
	Long float64
}
 
type Conversion interface {
	deg2num(t *Tile) (x int, y int)
	num2deg(t *Tile) (lat float64, long float64)
}
 
func (*Tile) Deg2num(t *Tile) (x int, y int) {
	x = int(math.Floor((t.Long + 180.0) / 360.0 * (math.Exp2(float64(t.Z)))))
	y = int(math.Floor((1.0 - math.Log(math.Tan(t.Lat*math.Pi/180.0)+1.0/math.Cos(t.Lat*math.Pi/180.0))/math.Pi) / 2.0 * (math.Exp2(float64(t.Z)))))
	return
}
 
func (*Tile) Num2deg(t *Tile) (lat float64, long float64) {
	n := math.Pi - 2.0*math.Pi*float64(t.Y)/math.Exp2(float64(t.Z))
	lat = 180.0 / math.Pi * math.Atan(0.5*(math.Exp(n)-math.Exp(-n)))
	long = float64(t.X)/math.Exp2(float64(t.Z))*360.0 - 180.0
	return lat, long
}

Java

 public class slippytest {
 public static void main(String[] args) {
   int zoom = 10;
   double lat = 47.968056d;
   double lon = 7.909167d;
   System.out.println("http://tile.openstreetmap.org/" + getTileNumber(lat, lon, zoom) + ".png");
 }
 public static String getTileNumber(final double lat, final double lon, final int zoom) {
   int xtile = (int)Math.floor( (lon + 180) / 360 * (1<<zoom) ) ;
   int ytile = (int)Math.floor( (1 - Math.log(Math.tan(Math.toRadians(lat)) + 1 / Math.cos(Math.toRadians(lat))) / Math.PI) / 2 * (1<<zoom) ) ;
    if (xtile < 0)
     xtile=0;
    if (xtile >= (1<<zoom))
     xtile=((1<<zoom)-1);
    if (ytile < 0)
     ytile=0;
    if (ytile >= (1<<zoom))
     ytile=((1<<zoom)-1);
    return("" + zoom + "/" + xtile + "/" + ytile);
   }
 }

Tile bounding box

  class BoundingBox {
    double north;
    double south;
    double east;
    double west;   
  }
  BoundingBox tile2boundingBox(final int x, final int y, final int zoom) {
    BoundingBox bb = new BoundingBox();
    bb.north = tile2lat(y, zoom);
    bb.south = tile2lat(y + 1, zoom);
    bb.west = tile2lon(x, zoom);
    bb.east = tile2lon(x + 1, zoom);
    return bb;
  }
 
  static double tile2lon(int x, int z) {
     return x / Math.pow(2.0, z) * 360.0 - 180;
  }
 
  static double tile2lat(int y, int z) {
    double n = Math.PI - (2.0 * Math.PI * y) / Math.pow(2.0, z);
    return Math.toDegrees(Math.atan(Math.sinh(n)));
  }

VB.Net

Private Function CalcTileXY(ByVal lat As Single, ByVal lon As Single, ByVal zoom As Long) As Point
   CalcTileXY.X  = CLng(Math.Floor((lon + 180) / 360 * 2 ^ zoom))
   CalcTileXY.Y = CLng(Math.Floor((1 - Math.Log(Math.Tan(lat * Math.PI / 180) + 1 / Math.Cos(lat * Math.PI / 180)) / Math.PI) / 2 * 2 ^ zoom))
End Function

C#

 
public PointF WorldToTilePos(double lon, double lat, int zoom)
{
	PointF p = new Point();
	p.X = (float)((lon + 180.0) / 360.0 * (1 << zoom));
	p.Y = (float)((1.0 - Math.Log(Math.Tan(lat * Math.PI / 180.0) + 
		1.0 / Math.Cos(lat * Math.PI / 180.0)) / Math.PI) / 2.0 * (1 << zoom));
 
	return p;
}
 
public PointF TileToWorldPos(double tile_x, double tile_y, int zoom) 
{
	PointF p = new Point();
	double n = Math.PI - ((2.0 * Math.PI * tile_y) / Math.Pow(2.0, zoom));
 
	p.X = (float)((tile_x / Math.Pow(2.0, zoom) * 360.0) - 180.0);
	p.Y = (float)(180.0 / Math.PI * Math.Atan(Math.Sinh(n)));
 
	return p;
}

XSLT

Requires math extensions from exslt.org.

<xsl:transform
 xmlns:xsl="http://www.w3.org/1999/XSL/Transform"
 xmlns:m="http://exslt.org/math"
 extension-element-prefixes="m"
 version="1.0">
 
 <xsl:output method="text"/>
 <xsl:variable name="pi" select="3.14159265358979323846"/>
 
 <xsl:template name="tiley">
 	<xsl:param name="lat"/>
 	<xsl:param name="zoomfact"/>
 	<xsl:variable name="a" select="($lat * $pi) div 180.0"/>
 	<xsl:variable name="b" select="m:log(m:tan($a) + (1.0 div m:cos($a)))"/>
 	<xsl:variable name="c" select="(1.0 - ($b div $pi)) div 2.0"/>
	<xsl:value-of select="floor($c * $zoomfact)"/>
 </xsl:template>
 
 <xsl:template name="tilename">
 	<xsl:param name="lat"/>
 	<xsl:param name="lon"/>
 	<xsl:param name="zoom" select="10"/>
 	<xsl:variable name="zoomfact" select="m:power(2,$zoom)"/>
 	<xsl:variable name="x" select="floor((360.0 + ($lon * 2)) * $zoomfact div 720.0)"/>
 	<xsl:variable name="y">
 		<xsl:call-template name="tiley">
 			<xsl:with-param name="lat" select="$lat"/>
 			<xsl:with-param name="zoomfact" select="$zoomfact"/>
 		</xsl:call-template>
 	</xsl:variable>
 	<xsl:value-of select="concat($zoom,'/',$x,'/',$y)"/>
 </xsl:template>
 
 <xsl:template match="/">
 	<xsl:call-template name="tilename">
 		<xsl:with-param name="lat" select="49.867731999999997"/>
 		<xsl:with-param name="lon" select="8.6295369999999991"/>
 		<xsl:with-param name="zoom" select="14"/>
 	</xsl:call-template>
 </xsl:template>
</xsl:transform>

Haskell

-- https://github.com/j4/HSlippyMap
 
long2tilex lon z = floor((lon + 180.0) / 360.0 * (2.0 ** z))
 
lat2tiley lat z = floor((1.0 - log( tan(lat * pi/180.0) + 1.0 / cos(lat * pi/180.0)) / pi) / 2.0 * (2.0 ** z))
 
tilex2long x z = x / (2.0 ** z) * 360.0 - 180
 
tiley2lat y z = 180.0 / pi * atan(0.5 * (exp(n) - exp(-n)))
        where
                n = pi - 2.0 * pi * y / (2.0 ** z)
 
-- Example
main = do
        --print $ long2tilex 2.2712 17
        --print $ lat2tiley 48.8152 17
        --print $ tilex2long 66362 17
        --print $ tiley2lat 45115 17
        putStrLn "gps: (lat=48.8152,long=2.2712)"
        putStrLn $ "http://tile.openstreetmap.org/17/" ++ show x ++ "/" ++ show y ++ ".png"
        where
                z = 17
                x = long2tilex 2.2712 z
                y = lat2tiley 48.8152 z

Scala

import scala.math._
 
case class Tile(x: Int,y: Int, z: Short){
  def toLatLon = new LatLonPoint(
    toDegrees(atan(sinh(Pi * (1.0 - 2.0 * y.toDouble / (1<<z))))), 
    x.toDouble / (1<<z) * 360.0 - 180.0,
    z)
  def toURI = new java.net.URI("http://tile.openstreetmap.org/"+z+"/"+x+"/"+y+".png")
}
 
case class LatLonPoint(lat: Double, lon: Double, z: Short){
  def toTile = new Tile(
    ((lon + 180.0) / 360.0 * (1<<z)).toInt,
    ((1 - log(tan(toRadians(lat)) + 1 / cos(toRadians(lat))) / Pi) / 2.0 * (1<<z)).toInt, 
    z)
}
 
//Usage:
val point = LatLonPoint(51.51202,0.02435,17)
val tile = point.toTile
// ==> Tile(65544,43582,17)
val uri = tile.toURI
// ==> http://tile.openstreetmap.org/17/65544/43582.png

Revolution/Transcript


function osmTileRef iLat, iLong, iZoom --> part path
   local n, xTile, yTile
   put (2 ^ iZoom) into n
   put (iLong + 180) / 360 * n into xTile
   multiply iLat by (pi / 180) -- convert to radians
   put ((1 - ln(tan(iLat) + 1 / cos(iLat)) / pi) / 2) * n into yTile
   return "/" & iZoom & "/" & trunc(xTile) & "/" & trunc(yTile)
end osmTileRef

function osmTileCoords xTile, yTile, iZoom --> coordinates
   local twoPzoom, iLong, iLat, n
   put (2 ^ iZoom) into twoPzoom
   put xTile / twoPzoom * 360 - 180 into iLong
   put pi - 2 * pi * yTile / twoPzoom into n
   put "n1=" && n
   put 180 / pi * atan(0.5 * (exp(n) - exp(-n))) into iLat 
   return iLat & comma & iLong
end osmTileCoords

Mathematica

Deg2Num[lat_, lon_, zoom_] := 
 {IntegerPart[(2^(-3 + zoom)*(180 + lon))/45], IntegerPart[2^(-1 + zoom)*(1 - Log[Sec[Degree*lat] + Tan[Degree*lat]]/Pi)]}
Num2Deg[xtile_,ytile_,zoom_] := 
 {ArcTan[Sinh[Pi*(1 - 2*(ytile/2^zoom))]]/Degree, (xtile/2^zoom)*360 - 180} // N

Tcl

First of all, you need to use the package map::slippy from Tcllib:

package require map::slippy

Lat./lon. to tile number

map::slippy geo 2tile [list $zoom $lat $lon]

Tile number to lat/lon

map::slippy tile 2geo [list $zoom $row $col]

Pascal

(translated from the Pythoncode above to Pascal)

Coordinates to tile numbers

uses {...}, Math;
{...}
var
  zoom: Integer;
  lat_rad, lat_deg, lon_deg, n: Real;
begin
  lat_rad := DegToRad(lat_deg);
  n := Power(2, zoom);
  xtile := Trunc(((lon_deg + 180) / 360) * n);
  ytile := Trunc((1 - (ln(Tan(lat_rad) + (1 /Cos(lat_rad))) / Pi)) / 2 * n);
end;

Tile numbers to coordinates

uses {...}, Math;
{...}
var
  lat_rad, n: Real;
begin
  n := Power(2, zoom);
  lat_rad := Arctan (Sinh (Pi * (1 - 2 * ytile / n)));
  lat_deg := RadtoDeg (lat_rad);
  lon_deg := xtile / n * 360.0 - 180.0;
end;

R

Coordinates to tile numbers

deg2num<-function(lat_deg, lon_deg, zoom){
  lat_rad <- lat_deg * pi /180
  n <- 2.0 ^ zoom
  xtile <- floor((lon_deg + 180.0) / 360.0 * n)
  ytile = floor((1.0 - log(tan(lat_rad) + (1 / cos(lat_rad))) / pi) / 2.0 * n)
  return( c(xtile, ytile))
#  return(paste(paste("http://a.tile.openstreetmap.org", zoom, xtile, ytile, sep="/"),".png",sep=""))
}

Bourne shell with Awk

Tile numbers to lat./lon. / Coordinates to tile numbers / Sample of usage, with optional tms-format support

xtile2long()
{
 xtile=$1
 zoom=$2
 echo "${xtile} ${zoom}" | awk '{printf("%.9f", $1 / 2.0^$2 * 360.0 - 180)}'
} 
 
long2xtile()  
{ 
 long=$1
 zoom=$2
 echo "${long} ${zoom}" | awk '{ xtile = ($1 + 180.0) / 360 * 2.0^$2; 
  xtile+=xtile<0?-0.5:0.5;
  printf("%d", xtile ) }'
}
 
ytile2lat()
{
 ytile=$1;
 zoom=$2;
 tms=$3;
 if [ ! -z "${tms}" ]
 then
 #  from tms_numbering into osm_numbering
  ytile=`echo "${ytile}" ${zoom} | awk '{printf("%d\n",((2.0^$2)-1)-$1)}'`;
 fi
 lat=`echo "${ytile} ${zoom}" | awk -v PI=3.14159265358979323846 '{ 
       num_tiles = PI - 2.0 * PI * $1 / 2.0^$2;
       printf("%.9f", 180.0 / PI * atan2(0.5 * (exp(num_tiles) - exp(-num_tiles)),1)); }'`;
 echo "${lat}";
}
 
lat2ytile() 
{ 
 lat=$1;
 zoom=$2;
 tms=$3;
 ytile=`echo "${lat} ${zoom}" | awk -v PI=3.14159265358979323846 '{ 
   tan_x=sin($1 * PI / 180.0)/cos($1 * PI / 180.0);
   ytile = (1 - log(tan_x + 1/cos($1 * PI/ 180))/PI)/2 * 2.0^$2; 
   ytile+=ytile<0?-0.5:0.5;
   printf("%d", ytile ) }'`;
 if [ ! -z "${tms}" ]
 then
  #  from oms_numbering into tms_numbering
  ytile=`echo "${ytile}" ${zoom} | awk '{printf("%d\n",((2.0^$2)-1)-$1)}'`;
 fi
 echo "${ytile}";
}
# ------------------------------------
# Sample of use: 
# Position Brandenburg Gate, Berlin
# ------------------------------------
LONG=13.37771496361961;
LAT=52.51628011262304;
ZOOM=17;
TILE_X=70406;
TILE_Y=42987; 
TILE_Y_TMS=88084;
TMS=""; # when NOT empty: tms format assumed
# ------------------------------------
# assume input/output of y is in oms-format:
LONG=$( xtile2long ${TILE_X} ${ZOOM} );
LAT=$( ytile2lat ${TILE_Y} ${ZOOM} ${TMS} );
# Result should be longitude[13.375854492] latitude[52.517892228]
TILE_X=$( long2xtile ${LONG} ${ZOOM} );
TILE_Y=$( lat2ytile ${LAT} ${ZOOM} ${TMS} );
# Result should be x[70406] y_oms[42987] 
# ------------------------------------
# assume input/output of y is in tms-format:
TMS="tms";
TILE_Y_TMS=$( lat2ytile ${LAT} ${ZOOM} ${TMS} );
LAT_TMS=$( ytile2lat ${TILE_Y_TMS} ${ZOOM} ${TMS} );
echo "Result should be y_oms[${TILE_Y}] latitude[${LAT}] ; y_tms[${TILE_Y_TMS}] latitude_tms[${LAT_TMS}] "
# latitude and latitude_tms should have the same value ; y_oms and y_tms should have the given start values:
# Result should be y_oms[42987] latitude[52.517892228] ; y_tms[88084] latitude_tms[52.517892228]
# ------------------------------------

Tile bounding box and center

n=$(ytile2lat `expr ${TILE_Y}` ${ZOOM})
s=$(ytile2lat `expr ${TILE_Y} + 1` ${ZOOM})
e=$(xtile2long `expr ${TILE_X} + 1` ${ZOOM})
w=$(xtile2long `expr ${TILE_X}` ${ZOOM})
 
echo "bbox=$w,$s,$e,$n" 
echo "-I-> Result should be [bbox=13.375854492,52.516220864,13.378601074,52.517892228]";
 
center_lat=`echo "$s $n" | awk '{printf("%.8f", ($1 + $2) / 2.0)}'`
center_lon=`echo "$w $e" | awk '{printf("%.8f", ($1 + $2) / 2.0)}'`
 
echo "center=$center_lat,$center_lon"
echo "-I-> Result should be [center=52.51705655,13.37722778]";

tileinfo.sh

Based on the above functions, this script can show some basic information about a tile position:

Calling tileinfo.sh without parameters shows a Help-Text:

-H-> Usage: tileinfo zoom-level x-position y-position [tms]
 : if 'x/y-position' are double numbers	        : Convert to tiles using the zoom-level
 : if 'x/y-position' are integers numbers	: Convert to Latitude/Longitude using the zoom-level
 : if [tms] is set the tile y-position will be assumed as in the tms-format, otherwise as the osm-format
-H-> Sample: Lat/Long to Tile: tileinfo 17 13.37771496361961 52.51628011262304
-H-> Sample: Osm-Tile to Lat/Long: tileinfo 17 70406 42987
-H-> Sample: Tms-Tile to Lat/Long: tileinfo 17 70406 88084 tms
-H-> Sample: Test-Position: tileinfo -test [will show position of the Brandenburger Tor, Berlin Germany]
OR to create a MBTiles-Database 
 ---> [call from the main Directory of the tiles - where the sub-Directories are called '0,1,..,19']
-H-> MBTiles: Create from Osm Tiles-Diretory: tileinfo -mbtiles
-H-> MBTiles: Create from Tms Tiles-Diretory: tileinfo -mbtiles tms
ALSO possible
-H-> MBTiles: Create using spatialite instead of sqlite3: tileinfo -mbtiles -spatialite [tms] 
 ---> the created Database will be stored in the main Directory of the tiles.
 -----------------------------------------------------------------------------------------
 -> Problem: gdal2tiles creates the 'mercator' tiles (for osm) with tms-numbering
-> 0.png /1/0/0.png = North America [osm] ; South America [tms]
-> 1.png /1/0/1.png = South America [osm] ; North America [tms]
-> 0.png /1/1/0.png = Northern Asia [osm] ; Southern Asia [tms]
-> 1.png /1/1/1.png= Southern Asia [osm] ; Northern Asia [tms]
 -----------------------------------------------------------------------------------------
 ---> the .mbtiles Databases allways uses the tms-numbering for the y-position.
 ---> Sqlite is installed. Spatialite is installed.
#!/bin/bash
#################################################
# Copyright (c) 2013, Mark Johnson
#
#    This program is free software: you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation, either version 3 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
#    See the GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program.  If not, see <http://www.gnu.org/licenses/>.
#******************************************************************************
# http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames
# http://osm.wff.ch/calc.htm
# http://wiki.osgeo.org/wiki/Tile_Map_Service_Specification
# https://github.com/mapbox/mbtiles-spec/blob/master/1.2/spec.md
# -----------------------------------------------------------------------
# Returns basic information about a tile-position and
# can create a MBTiles sqlite Database from an existing tiles
# -----------------------------------------------------------------------
# Version: 2013-12-24
# - inserting into mbtiles after reading every x directory
# - removed usage of spatialite
# - added support of soft-links
# - blank-images with tile_id as -R-G-B.rgb [will only be stored once]
# Version: 2013-10-27
# - added support for file-extentiong 'jpg','jpeg' togeather with 'png'
# Version: 2013-10-26
# - added sorting of sql statement, so that the z/x/y are sorted from top to bottom if SORT is installed
# - extensive testing with android application 'geopaparazzi'
# Version: 2013-10-25
# - corrected format= 'jpg' instead of 'png' if CONVERT is installed
# Version: 2013-09-05
# - corrected confusion in osm/tms numbering
# - removed rounding from lat/long to tile
# - added print of commands to download and georectify a queried tile
#  to 'WGS 84 / World Mercator' (3395)
# - convert png-tiles to jpg if 'convert' is installed [mbtiles]
# -----------------------------------------------------------------------
BENE="'Tutto bene!'";
BENENO="'No bene!'";
HABE_FERTIG="Ich habe fertig.";
BASE_NAME=`basename $0`;
BASE_NAME_CUT=`basename $0 | cut -d '.' -f1`;
MESSAGE_TYPE="-I->";
exit_rc=0;
tms="osm":
TILE_SIZE=256;
declare -a PIXEL_SIZE_ARRAY=( );
# All values based on the position: 13.37771496361961 52.51628011262304 [Brandenburg Gate, Berlin, Germany]
POSITION_BASE_ARRAY+=('Pixel Size = (156376.587086746090790,-156376.587086746090790)*256= 40032406.294 meters - gdal=yes'); # Zoom 00
POSITION_BASE_ARRAY+=('Pixel Size = (78188.293543373045395,-78188.293543373045395)*256= 20016203.147 meters - gdal=yes'); # Zoom 01
POSITION_BASE_ARRAY+=('Pixel Size = (39059.204695334061398,-39059.204695334061398)*256= 9999156.402 meters - gdal=yes'); # Zoom 02
POSITION_BASE_ARRAY+=('Pixel Size = (19513.214221145241027,-19513.214221145241027)*256= 4995382.840 meters - gdal=yes');  # Zoom 03
POSITION_BASE_ARRAY+=('Pixel Size = (9769.632134524423236,-9769.632134524423236)*256= 2501025.826 meters - gdal=yes');  # Zoom 04
POSITION_BASE_ARRAY+=('Pixel Size = (4885.857856551832811,-4885.857856551832811)*256= 1250779.611 meters - gdal=yes');  # Zoom 05
POSITION_BASE_ARRAY+=('Pixel Size = (2442.692472795480171,-2442.692472795480171)*256= 625329.273 meters - gdal=yes');  # Zoom 06
POSITION_BASE_ARRAY+=('Pixel Size = (1221.408542870289011,-1221.408542870289011)*256= 312680.586 meters - gdal=yes');  # Zoom 07
POSITION_BASE_ARRAY+=('Pixel Size = (610.719455956439333,-610.719455956439333)*256= 156344.180 meters - gdal=yes');  # Zoom 08
POSITION_BASE_ARRAY+=('Pixel Size = (305.363474400363145,-305.363474400363145)*256= 78173.049 meters - gdal=yes');  # Zoom 09
POSITION_BASE_ARRAY+=('Pixel Size = (152.682667445300353,-152.682667445300353)*256= 39086.762 meters - gdal=yes');  # Zoom 10
POSITION_BASE_ARRAY+=('Pixel Size = (76.341565429802998,-76.341565429802998)*256= 19543.440 meters - gdal=yes');  # Zoom 11
POSITION_BASE_ARRAY+=('Pixel Size = (38.170956013236221,-38.170956013236221)*256= 9771.764 meters - gdal=yes');  # Zoom 12
POSITION_BASE_ARRAY+=('Pixel Size = (19.085492344712762,-19.085492344712762)*256= 4885.886 meters - gdal=yes');  # Zoom 13
POSITION_BASE_ARRAY+=('Pixel Size = (9.542742608409494,-9.542742608409494)*256= 2442.942 meters - gdal=yes');  # Zoom 14
POSITION_BASE_ARRAY+=('Pixel Size = (4.771372073263216,-4.771372073263216)*256= 1221,471 meters - gdal=yes');  # Zoom 15
POSITION_BASE_ARRAY+=('Pixel Size = (2.385685756825894,-2.385685756825894)*256= 610.735 meters - gdal=yes');  # Zoom 16
POSITION_BASE_ARRAY+=('Pixel Size = (1.192842719290223,-1.192842719290223)*256= 305.367 meters - gdal=yes');  # Zoom 17
POSITION_BASE_ARRAY+=('Pixel Size = (0.596421364337474,-0.596421364337474)*256= 152.683 meters - gdal=yes');  # Zoom 18
POSITION_BASE_ARRAY+=('Pixel Size = (0.298210792187396,-0.298210792187396)*256= 76.341 meters - gdal=yes');  # Zoom 19
POSITION_BASE_ARRAY+=('Pixel Size = (0.149105396388009,-0.149105396388009)*256= 38.170 meters - gdal=yes');  # Zoom 20
POSITION_BASE_ARRAY+=('Pixel Size = (0.074552544930600,-0.074552544930600)*256= 19.085 meters - gdal=yes');  # Zoom 21
POSITION_BASE_ARRAY+=('Pixel Size = (0.037276174719086,-0.037276174719086)*256= 9,542 meters - gdal=yes');  # Zoom 22
POSITION_BASE_ARRAY+=('Pixel Size = (0.018638262558992,-0.018638262558992)*256= 4.771 meters - gdal=yes');  # Zoom 23
POSITION_BASE_ARRAY+=('Pixel Size = (0.009319239432482,-0.009319239432482)*256= 2.385 meters - gdal=no');  # Zoom 24
POSITION_BASE_ARRAY+=('Pixel Size = (0.004659619889283,-0.004659619889283)*256= 1.192 meters - gdal=no');  # Zoom 25
POSITION_BASE_ARRAY+=('Pixel Size = (0.002329523079167,-0.002329523079167)*256= 0.596 meters - gdal=no');  # Zoom 26
POSITION_BASE_ARRAY+=('Pixel Size = (0.001164761532319,-0.001164761532319)*256=  0.298 meters - gdal=no');  # Zoom 27
POSITION_BASE_ARRAY+=('Pixel Size = (0.000582380761390,-0.000582380761390)*256= 0.149 meters - gdal=no ');  # Zoom 28
POSITION_BASE_ARRAY+=('Pixel Size = (0.000291299252023,-0.000291299252023)*256= 0.074 meters - gdal=no');  # Zoom 29
POSITION_BASE_ARRAY+=('Pixel Size = (0.000145758555015,-0.000145758555015)*256= 0.037 meters - gdal=no');  # Zoom 30
POSITION_BASE_ARRAY+=('[gdalwarp] Failed to compute GCP transform: Transform is not solvable');  # Zoom 31
# 0: 152.363223666284540
# 1:  152.363338110592139
#-------------------------------------------------------------------------------
# Tests whether *entire string* is numerical.
# In other words, tests for integer variable.
#-----------------------------------------------------------------------------
is_digit()
{        # if ! is_digit "blombo"
  [ $# -eq 1 ] || return 1;
  case $1 in
   *[!0-9]*|"")
    return 1;
   ;;
   *)
    return 0;
   ;;
  esac
}
#-------------------------------------------------------------------------------
#-- TMS Global Mercator to osm Global Mercator
#-------------------------------------------------------------------------------
flip_tmsosm()
{ # from tms/osm_numbering into osm/tms_numbering  - [tms z=2 x=1 y=3 == osm z=2 x=1 y=0]
 TILE_Y_OSM=$((2**${2}-${1}-1));
 echo ${TILE_Y_OSM};
}
# y = Math.pow(2, z) - y - 1; $((2**${2}-${1}-1))
#-------------------------------------------------------------------------------
#-- osm Global Mercator to TMS Global Mercator
#-------------------------------------------------------------------------------
#-- TMS Global Geodetic
#-------------------------------------------------------------------------------
lat2_geodetic_ytile()
{
 lat=$1
 zoom=$2
 echo "${lat} ${zoom} ${TILE_SIZE} " | awk '{ res = (180 /  $3) / (2.0^$2);
  ypixel =(90 + $1) / res;
  ytile = ypixel / $3;
  ytile = (ytile == int(ytile)) ? ytile : int(ytile)+1;
  printf("%d", ytile ) }'
}
#-------------------------------------------------------------------------------
ytile2_geodetic_lat()
{
 ytile=$1
 zoom=$2
 echo "${ytile} ${zoom} ${TILE_SIZE} " | awk '{ res = (180 /  $3) / (2.0^$2);
  lat =($1 * $3 * res) - 90;
  printf(""%.9f"", lat ) }'
}
# ------------------------------------
long2_geodetic_xtile()
{
 long=$1
 zoom=$2
 echo "${long} ${zoom} ${TILE_SIZE} " | awk '{ res = (180 /  $3) / (2.0^$2);
  xpixel =(90 + $1) / res;
  xtile = xpixel / $3;
  xtile = (xtile == int(xtile)) ? xtile : int(xtile)+1;
  printf("%d", xtile-1 ) }'
}
#-------------------------------------------------------------------------------
xtile2_geodetic_long()
{
 xtile=$1
 zoom=$2
 echo "${xtile} ${zoom} ${TILE_SIZE} " | awk '{ res = (180 /  $3) / (2.0^$2);
  long =($1 * $3 * res) - 180;
  printf(""%.9f"", long ) }'
}
# ------------------------------------
USE_SPATIALITE="";
TMS="";
TILE_EXTENTION="png";
TILE_EXTENTION_CONVERT="png";
MBTILE_EXTENTION="mbtiles";
CMD_SQLITE=`which sqlite3`;
CMD_SPATIALITE=`which spatialite`;
CMD_SPATIALITE=""; # do not use
CMD_CONVERT=`which convert`;
CMD_IDENTIFY=`which identify`;
CMD_SORT=`which sort`;
CONVERT_COMPRESS_PARM="-compress zip";
IDENTIFY_FORMAT_PARM=" -format %k";
CONVERT_FORMAT_PARM="-format '%[fx:int(255*p{0,0}.r)]-%[fx:int(255*p{0,0}.g)]-%[fx:int(255*p{0,0}.b)]' info: ";
# identify -format %k 19-281166-352774.png
# 19-281166-352774.tms BLOB sz=417 JPEG image
# cp ../19/281166/352774.png .
# convert 19-281166-352774.png -format "%[fx:int(255*p{10,10}.r)]-%[fx:int(255*p{10,10}.g)]-%[fx:int(255*p{10,10}.b)]" info:
# convert -format "%[fx:int(255*p{10,10}.r)]-%[fx:int(255*p{10,10}.g)]-%[fx:int(255*p{10,10}.b)] " 19-281626-352336.png info:
CMD_GDAL_TRANSLATE=`which gdal_translate`;
# uncomment the next line if spatialite should be used
# USE_SPATIALITE="true";
# ------------------------------------
LONG=13.37771496361961;
LAT=52.51628011262304;
# tileinfo 17 13.37771496361961 52.51628011262304
ZOOM=17;
TILE_X=70406;
TILE_Y_OSM=42988;
TILE_Y_TMS=88083;
TILE_X_GEOD=${TILE_X};
TILE_Y_GEOD=${TILE_Y_OSM};
# ------------------------------------
MODUS="0";
INPUT_PARMS="";
# ------------------------------------
if [ $# -gt "0" ]
then
 if [ $# -eq "1" ]
 then
  if [ $1 == "-test" ]
  then
   MODUS=1;
   INPUT_PARMS="Test-Position Brandenburg Gate z/x,y[${ZOOM}/${LONG},${LAT}] modus[${MODUS}]";
  else
   if is_digit "$1"
   then
    MODUS=1;
    ZOOM=$1;
    INPUT_PARMS="Test-Position Brandenburg Gate z/x,y[${ZOOM}/${LONG},${LAT}] modus[${MODUS}]";
   else
    OIFS=$IFS;
    IFS='/';
    ARRAY_01=( $1 );
    IFS=$OIFS;
    if [  "${#ARRAY_01[@]}" -gt "2" ]
    then
     # /17/70406/88085
     POS_Z=1;
     POS_X=2;
     POS_Y=3;
     if [  "${#ARRAY_01[@]}" -eq "3" ]
     then
      # 17/70406/88085
      POS_Z=0;
      POS_X=1;
      POS_Y=2;
     fi
     TMS="osm";
     MODUS=2; # z/x/y to long/lat
     ZOOM=${ARRAY_01[${POS_Z}]};
     TILE_X=${ARRAY_01[${POS_X}]};
     TILE_Y=${ARRAY_01[${POS_Y}]};
     OIFS=$IFS;
     IFS='.';
     ARRAY_01=( ${TILE_Y} ); # remove possible '.jpg or .png'
     IFS=$OIFS;
     TILE_Y=${ARRAY_01[0]};
     TILE_Y_OSM=${TILE_Y};
     TILE_Y_TMS=$( flip_tmsosm ${TILE_Y} ${ZOOM} ${TMS} );
     INPUT_PARMS="Input: Osm-Tile z/x/y[${ZOOM}/${TILE_X}/${TILE_Y_OSM}] tms[${TILE_Y_TMS}] modus[${MODUS}] ";
    fi
   fi
  fi
 fi
 if [ $1 == "-mbtiles" ]
 then
  TMS="osm";
  MODUS=3;
  INPUT_PARMS="Create sqlite MBTiles Database from the OSM-Directory [${PWD##*/}] modus[${MODUS}]";
  if [ $# -gt "1" ]
  then
   if [ $2 == "tms" ]
   then
    TMS="tms";
    INPUT_PARMS="Create sqlite MBTiles Database from the TMS-Directory [${PWD##*/}]";
   fi
   if [ $2 == "-spatialite" ]
   then
    TILE_TYPE="OSM-";
    if [ $# -gt "2" ]
    then
     if [ $3 == "tms" ]
     then
      TMS="tms";
      TILE_TYPE="TMS-";
     fi
    fi
    if [ ! -z "${CMD_SPATIALITE}" ]
    then
     CMD_SQLITE="${CMD_SPATIALITE}";
     # CMD_SQLITE="${CMD_SPATIALITE} -silent";
     INPUT_PARMS="Create spatialite MBTiles Database from the ${TILE_TYPE} Directory [${PWD##*/}]";
    else
     INPUT_PARMS="Create sqlite MBTiles Database from the ${TILE_TYPE} Directory [${PWD##*/}] - spatialite not found";
    fi
   fi
  fi
 fi
 if [ $# -eq "2" ] && [ "${2}" == "tms" ]
 then
  OIFS=$IFS;
  IFS='/';
  ARRAY_01=( $1 );
  IFS=$OIFS;
  if [  "${#ARRAY_01[@]}" -gt "2" ]
  then
   # /17/70406/88085
   POS_Z=1;
   POS_X=2;
   POS_Y=3;
   if [  "${#ARRAY_01[@]}" -eq "3" ]
   then
    # 17/70406/88085
    POS_Z=0;
    POS_X=1;
    POS_Y=2;
   fi
   TMS="tms";
   MODUS=2; # z/x/y to long/lat
   ZOOM=${ARRAY_01[${POS_Z}]};
   TILE_X=${ARRAY_01[${POS_X}]};
   TILE_Y=${ARRAY_01[${POS_Y}]};
   OIFS=$IFS;
   IFS='.';
   ARRAY_01=( ${TILE_Y} ); # remove possible '.jpg or .png'
   IFS=$OIFS;
   TILE_Y=${ARRAY_01[0]};
   TILE_Y_TMS=${TILE_Y};
   TILE_Y_OSM=$( flip_tmsosm ${TILE_Y} ${ZOOM} ${TMS} );
   INPUT_PARMS="Input: Tms-Tile z/x/y[${ZOOM}/${TILE_X}/${TILE_Y_TMS}] osm[${TILE_Y_OSM}] modus[${MODUS}]";
  fi
 fi
 if [ "${MODUS}" -eq "3" ] && [ -z "${CMD_SQLITE}" ]
 then
  MODUS="";
 fi
 if [ "${MODUS}" -ne "3" ]
 then
  if [ $# -eq "3" ] || [ $# -eq "4" ]
  then
   OIFS=$IFS;
   IFS='.';
   ARRAY_02=( $2 );
   ARRAY_03=( $3 );
   IFS=$OIFS;
   if [  "${#ARRAY_02[@]}" -gt "1" ] && [  "${#ARRAY_03[@]}" -gt "1" ]
   then
    # This are decimal numbers
    MODUS=1; # long/lat to z/x/y
    ZOOM=$1;
    LONG=$2;
    LAT=$3;
    if [ "$4" == "tms" ]
    then
     TMS="tms";
    else
     TMS="osm";
    fi
    TILE_X="";
    TILE_Y="";
    TILE_Y_OSM="";
    TILE_Y_TMS="";
    TILE_Y_GEOD="";
    INPUT_PARMS="Input: Position z/x,y[${ZOOM}/${LONG},${LAT}] tms[${TMS}] modus[${MODUS}]";
   else
    MODUS=2; # z/x/y to long/lat
    ZOOM=$1;
    LONG="";
    TILE_X=$2;
    if [ -z "$4" ]
    then
     TMS="osm";
     TILE_Y=$3;
     TILE_Y_OSM=$3;
     TILE_Y_TMS=$( flip_tmsosm ${TILE_Y} ${ZOOM} ${TMS} );
     INPUT_PARMS="Input: Osm-Tile z/x/y[${ZOOM}/${TILE_X}/${TILE_Y_OSM}] tms[${TILE_Y_TMS}] modus[${MODUS}]";
    fi
    if [ "$4" == "tms" ]
    then
     TMS="tms";
     TILE_Y=$3;
     TILE_Y_TMS=$3;
     TILE_Y_OSM=$( flip_tmsosm ${TILE_Y} ${ZOOM} ${TMS} );
     INPUT_PARMS="Input: Tms-Tile z/x/y[${ZOOM}/${TILE_X}/${TILE_Y_TMS}] osm[${TILE_Y_OSM}] modus[${MODUS}]";
    fi
    if [ "$4" == "geod" ]
    then
     TILE_Y=$3;
     TILE_Y_GEOD=$3;
     TILE_Y_TMS=$( flip_tmsosm ${TILE_Y} ${ZOOM} ${TMS});
     INPUT_PARMS="Input: Geodetic-Tile z/x/y[${ZOOM}/${TILE_X}/${TILE_Y_GEOD}] tms[${TILE_Y_TMS}] modus[${MODUS}]";
    fi
   fi
  fi
 fi
fi
# ------------------------------------
greater_than()
{
 max=$(echo $1 $2 | awk '{if ($1 >= $2) print $1; else print $2}')
 echo ${max};
}
# ------------------------------------
less_than()
{
 min=$(echo $1 $2 | awk '{if ($1 <= $2) print $1; else print $2}')
 echo ${min};
}
# ------------------------------------
xtile2long()
{
 xtile=$1;
 zoom=$2;
 echo "${xtile} ${zoom}" | awk '{printf("%.9f", $1 / 2.0^$2 * 360.0 - 180)}'
}
# ------------------------------------
long2xtile()
{
 long=$1;
 zoom=$2;
 #   xtile+=xtile<0?-0.5:0.5; # this should be not be rounded up
 echo "${long} ${zoom}" | awk '{ xtile = ($1 + 180.0) / 360 * 2.0^$2;
  printf("%d", xtile ) }'
}
# ------------------------------------
ytile2lat()
{
 ytile=$1;
 zoom=$2;
 tms=$3;
 if [ "${tms}" == "tms" ]
 then
 #  from tms_numbering into osm_numbering
  ytile_tms=${ytile};
  ytile_osm=$( flip_tmsosm ${ytile} ${zoom} );
 else
  ytile_osm=${ytile};
  ytile_tms=$( flip_tmsosm ${ytile} ${zoom} );
 fi
 lat=`echo "${ytile_osm} ${zoom}" | awk -v PI=3.14159265358979323846 '{
       num_tiles = PI - 2.0 * PI * $1 / 2.0^$2;
       printf("%.9f", 180.0 / PI * atan2(0.5 * (exp(num_tiles) - exp(-num_tiles)),1)); }'`;
 echo "${lat}";
}
# ------------------------------------
lat2ytile()
{
 lat=$1;
 zoom=$2;
 tms=$3;
 #  ytile_osm+=ytile_osm<0?-0.5:0.5; # this should be not be rounded up
 ytile_osm=`echo "${lat} ${zoom}" | awk -v PI=3.14159265358979323846 '{
   tan_x=sin($1 * PI / 180.0)/cos($1 * PI / 180.0);
   ytile_osm = (1 - log(tan_x + 1/cos($1 * PI/ 180))/PI)/2 * 2.0^$2;
   printf("%d", ytile_osm ) }'`;
 if [ "${tms}" == "tms" ]
 then
  #  from osm_numbering into tms_numbering
   ytile=$( tile_osmtms ${ytile_osm} ${zoom} );
 else
  ytile=${ytile_osm};
 fi
 echo "${ytile}";
 return;
}
# ------------------------------------
#  Crude map scale
# Beware: True scale value will depend on your specific monitor size and current resolution.
# This is just a guide, map scale is only really meaningful for the printed page.
# This dpi to meters-on-the-ground "pixel factor" is based on someone's ancient CRT monitor, perhaps a 15" model at 800x600 resolution
calc_webtile_scale()
{
  lat=$1
  scale=$2
  echo "$lat $scale" | awk -v PI=3.14159265358979323846 \
    -v PIXELFACT=2817.947378 '{
       a = 6378137.0;
       printf("%.0f", (a * 2*PI * cos($1 * PI/180) * PIXELFACT) / (256 * 2^$2))
 
    }'
}
# ------------------------------------
as_bounding_box()
{
 zoom=$1;
 xtile=$2;
 ytile=$3;
 tms="";
 if [ $# -gt "3" ]
 then
  tms=$4;
 fi
 n=$(ytile2lat `expr ${ytile}` ${zoom} ${tms} );
 s=$(ytile2lat `expr ${ytile} + 1` ${zoom} ${tms} );
 e=$(xtile2long `expr ${xtile} + 1` ${zoom} ${tms} );
 w=$(xtile2long `expr ${xtile}` ${zoom} ${tms} );
 echo "bbox=$w,$s,$e,$n";
}
# ------------------------------------
as_center()
{
 zoom=$1;
 xtile=$2;
 ytile=$3;
 tms=$4;
 n=$(ytile2lat `expr ${ytile}` ${zoom} ${tms} );
 s=$(ytile2lat `expr ${ytile} + 1` ${zoom} ${tms} );
 e=$(xtile2long `expr ${xtile} + 1` ${zoom} ${tms} );
 w=$(xtile2long `expr ${xtile}` ${zoom} ${tms} );
 center_lat=`echo "$s $n" | awk '{printf("%.8f", ($1 + $2) / 2.0)}'`
 center_lon=`echo "$w $e" | awk '{printf("%.8f", ($1 + $2) / 2.0)}'`
 echo "center=$center_lon,$center_lat"
}
# ------------------------------------
find_image_type()
{ # find -L . -maxdepth 3 -mindepth 2   -type f -name *.png  | wc -l
 find_count=`find -L . -maxdepth 3 -mindepth 2 -type f  -name *.${TILE_EXTENTION}  | wc -l`;
 if [ "${find_count}" -le 0 ]
 then
  TILE_EXTENTION="jpg";
  find_count=`find -L . -maxdepth 3 -mindepth 2 -type f  -name *.${TILE_EXTENTION}  | wc -l`;
  if [ "${find_count}" -le 0 ]
  then
   TILE_EXTENTION="jpeg";
   find_count=`find -L . -maxdepth 3 -mindepth 2 -type f  -name *.${TILE_EXTENTION}  | wc -l`;
  fi
 fi
 echo "-I-> $0 find_image_type: MBTiles[${this_dir}] ; Searching for tiles as [*.${TILE_EXTENTION}] - count[${find_count}]";
 if [ "${find_count}" -le 0 ]
 then
  return 1;
 else
  return 0;
 fi
}
# ------------------------------------
loop_tile_dir()
{
 this_dir=$1;
 tms=$2;
 if [ -f "${this_dir}.sql" ]
 then
  rm ${this_dir}.sql;
 fi
 DIR_TYPE="Osm";
 if [ "${tms}" == "tms" ]
 then
  DIR_TYPE="Tms";
 fi
# ------------------------------------
 min_y=1000.0;
 max_y=-1000.0;
 max_x=$max_y;
 min_x=$min_y;
 max_z=0;
 min_z=30;
 z_tile_prev="";
 x_tile_prev="";
 tile_count=0;
 echo "-I-> $0 loop_tile_dir: MBTiles[${this_dir}] ; Searching for tiles as [*.${TILE_EXTENTION}] - parm_tms[${tms}]";
 if [ ! -z "${CMD_CONVERT}" ]  && [ "${TILE_EXTENTION}" = "png" ]
 then
  TILE_EXTENTION_CONVERT="jpg";
  CONVERT_COMPRESS_PARM="-compress zip"; # this is (I think) only for tif -seems to work with jpg
  echo -e "-I-> $0 loop_tile_dir: will convert .png to .jpg with 'convert' using: [${CONVERT_COMPRESS_PARM}] \n-I-> \t This will reduce the size of the Database by about 80%.";
 else
  if [ "${TILE_EXTENTION}" = "png" ]
  then
   echo -e "-W-> $0 'convert' is needed to convert from png to the jpg format. The Database could be 80% smaller with jpg files. [Install ImageMagik package]";
  fi
 fi
 if [ ! -z "${CMD_SORT}" ]
 then
  SORT_PARM="--key=24,90"; # sort from 24 [ ... '17-70382-43019.osm' to about here]
  echo -e "-I-> $0 loop_tile_dir: will sort the insert sql-commands so that the z/x/y values will be inserted from top to bottom with 'sort' using: [${SORT_PARM}].";
 else
  echo -e "-W-> $0 'sort' is needed to sort the sql, so that the z/x/y values will be inserted from top to bottom. "
 fi
 for I in `find -L . -maxdepth 3 -mindepth 2 -type f  -name *.${TILE_EXTENTION}`
 do
  ytile=`basename $I .${TILE_EXTENTION}`
  DIRNAME=`dirname $I | sed s/^.//`;
  z_tile=`echo $DIRNAME | cut -d "/" -f 2`;
  if  is_digit "$z_tile"
  then
   x_tile=`echo $DIRNAME | cut -d "/" -f 3`;
   tile_id="${z_tile}-${x_tile}-${ytile}";
   if [ "${tms}" == "tms" ]
   then
    tile_id="${tile_id}.tms";
    ytile_tms=${ytile};
    ytile_osm=$( flip_tmsosm ${ytile} ${z_tile} );
   fi
   if [ "${tms}" == "osm" ]
   then
    tile_id="${tile_id}.osm";
    ytile_osm=${ytile};
    ytile_tms=$( flip_tmsosm ${ytile} ${z_tile} );
   fi
   if  is_digit "$z_tile"  &&  is_digit "$x_tile"  && is_digit "$ytile_osm"  && is_digit "$ytile_tms"
   then
    ((tile_count++));
    if [ "${z_tile}" != "${z_tile_prev}" ]
    then
     echo "-I-> $0 loop_tile_dir: reading[${this_dir}/${z_tile}/${x_tile}] ; tiles found: [${tile_count}]";
     z_tile_prev=${z_tile};
    fi
    if [ "${x_tile}" != "${x_tile_prev}" ]
    then
     if [ ! -z "${x_tile_prev}" ]
     then
      echo "-I-> $0 loop_tile_dir: reading[${this_dir}/${z_tile}/${x_tile}] ; tiles found: [${tile_count}]";
      write_mbtiles "${this_dir}" "${min_y}" "${max_y}" "${min_x}" "${max_x}" "${min_z}" "${max_z}" "${x_tile_prev}";
      if [ ! -f "${this_dir}.sql" ]
      then
       echo "BEGIN TRANSACTION;"  > ${this_dir}.sql;
      fi
     fi
     x_tile_prev=${x_tile};
    fi
    # This will calculate the min/max values
    n=$(ytile2lat `expr ${ytile_osm}` ${z_tile} "osm" );
    s=$(ytile2lat `expr ${ytile_osm} + 1` ${z_tile} "osm" );
    e=$(xtile2long `expr ${x_tile} + 1` ${z_tile} "osm" );
    w=$(xtile2long `expr ${x_tile}` ${z_tile} "osm" );
    min_y=$(less_than "${s}" "${min_y}" );
    max_y=$(greater_than "${n}" "${max_y}" );
    min_x=$(less_than "${w}" "${min_x}" );
    max_x=$(greater_than "${e}" "${max_x}" );
    min_z=$(less_than "${z_tile}" "${min_z}" );
    max_z=$(greater_than "${z_tile}" "${max_z}" );
    BBOX="$min_x,$min_y,$max_x,$max_y - $min_z,$max_z";
    if [ ! -z "${CMD_IDENTIFY}" ]
    then
     RGB_COUNT=`${CMD_IDENTIFY} ${IDENTIFY_FORMAT_PARM} ${I}`;
     if [ "${RGB_COUNT}" -eq 1 ]
     then
      tile_id_save=${tile_id};
      # this image has one distinct colour
      RGB_COUNT=`${CMD_CONVERT} ${I} ${CONVERT_FORMAT_PARM}`;
      RGB_COUNT=`echo ${RGB_COUNT} | cut -d "'" -f 2`;
      OIFS=$IFS;
      IFS='-';
      ARRAY_01=( ${RGB_COUNT} );
      IFS=$OIFS;
      R_VALUE=${ARRAY_01[0]};
      G_VALUE=${ARRAY_01[1]};
      B_VALUE=${ARRAY_01[2]};
      R_VALUE=$(printf "%02x" "${R_VALUE}");
      G_VALUE=$(printf "%02x" "${G_VALUE}");
      B_VALUE=$(printf "%02x" "${B_VALUE}");
      # store as 'e6-bc-95.rgb'
      tile_id="${R_VALUE}-${G_VALUE}-${B_VALUE}.rgb";
      echo "-I-> blank rgb: old[${tile_id_save}] new[${tile_id}] [${RGB_COUNT}] r[${R_VALUE}] g[${G_VALUE}] b[${B_VALUE}]";
     fi
    fi
    # make this sortable
    sql_insert_map="INSERT OR REPLACE INTO map (tile_id,zoom_level,tile_column,tile_row,grid_id) VALUES('${tile_id}',${z_tile},${x_tile},${ytile_tms},'');";
    input_file=${I};
    if [ ! -z "${CMD_CONVERT}" ] && [ "${TILE_EXTENTION}" = "png" ]
    then
     input_file="convert.${TILE_EXTENTION_CONVERT}";
     OUTPUTPARMS="${CONVERT_COMPRESS_PARM} $I ${input_file}";
     ${CMD_CONVERT} ${OUTPUTPARMS};
     local_rc=$?;
     if [ "${local_rc}" -ne 0 ]
     then
      input_file=${I};
     fi
    fi
    # sql_insert_image="INSERT OR REPLACE INTO images VALUES(X'`hexdump -ve '1/1 "%.2x"' ${input_file}`', '${tile_id}');"
    # make this sortable
    sql_insert_image="INSERT OR REPLACE INTO images (tile_id,tile_data) VALUES('${tile_id}',X'`hexdump -ve '1/1 "%.2x"' ${input_file}`');"
    if [ ! -z "${CMD_CONVERT}" ]  && [ "${TILE_EXTENTION}" = "png" ]
    then
     rm ${input_file};
    fi
    # echo "-I-> sql_insert_map[${sql_insert_map}] ;";
    # echo "-I-> sql_insert_image[${sql_insert_image}] ;";
    # echo "-I-> I[${I}] ; DIRNAME[${DIRNAME}] z_tile[${z_tile}] x_tile[${x_tile}] ytile_osm[${ytile_osm}] ytile_tms[${ytile_tms}] tile_id[${tile_id}] BBOX[${BBOX}";
    if [ ! -f "${this_dir}.sql" ]
    then
     echo "-I-> $0 loop_tile_dir: creating[${this_dir}.sql] ; ";
    fi
    echo ${sql_insert_map}  >> ${this_dir}.sql;
    echo ${sql_insert_image}>> ${this_dir}.sql;
   else
    tile_text="${z_tile}-${x_tile}-${ytile_tms}-${ytile_osm}";
    echo "-W-> $0 loop_tile_dir: non-numeric values found[${tile_text}] ; format should be : [z-x-y_tms-y_osm]";
   fi
  else
   echo "-W-> $0 loop_tile_dir: non-numeric directory found[${z_tile}] skipping this [${DIRNAME}] ";
  fi
 done
 x_tile_prev="final";
 write_mbtiles "${this_dir}" "${min_y}" "${max_y}" "${min_x}" "${max_x}" "${min_z}" "${max_z}" "${x_tile_prev}";
 echo "-I-> $0 loop_tile_dir: reading[compleated] ; tiles found: [${tile_count}]";
}
# ------------------------------------
duplicate_images_mbtiles()
{
 this_dir=$1;
 echo "
   SELECT tile_id,tile_data FROM images WHERE tile_data IN (SELECT tile_data FROM images GROUP BY tile_data HAVING (COUNT(hex(tile_data))>1));
   " > ${this_dir}.sql;
  RC_RESULT=`cat ${this_dir}.sql | ${CMD_SQLITE} ${this_dir}.${MBTILE_EXTENTION}`;
  echo ${RC_RESULT};
}
# ------------------------------------
write_mbtiles()
{
 this_dir=$1;
 min_y=$2;
 max_y=$3;
 min_x=$4;
 max_x=$5;
 min_z=$6;
 max_z=$7;
 x_tile_prev=$8;
 if [ -f "${this_dir}.sql" ]
 then
  if [ ! -z "${CMD_SORT}" ]
  then
   # echo "-I-> moving: [${this_dir}.sql to ${this_dir}.sql.sort]";
   mv ${this_dir}.sql ${this_dir}.sql.sort
   # echo "-I-> sorting: ${this_dir}.sql [${CMD_SORT} -o ${this_dir}.sql ${SORT_PARM} ${this_dir}.sql.sort]";
   ${CMD_SORT} -o ${this_dir}.sql ${SORT_PARM} ${this_dir}.sql.sort
   local_rc=$?;
   if [ "${local_rc}" -eq 0 ] &&  [ -f "${this_dir}.sql" ]
   then
    # echo "-I-> removing: [${this_dir}.sql.sort]";
    rm ${this_dir}.sql.sort
   else
    mv ${this_dir}.sql.sort ${this_dir}.sql
   fi
  fi
  # ------------------------------------
  center_lat=`echo "${min_y} ${max_y}" | awk '{printf("%.8f", ($1 + $2) / 2.0)}'`
  center_lon=`echo "${min_x} ${max_x}" | awk '{printf("%.8f", ($1 + $2) / 2.0)}'`
  if [ "${x_tile_prev}" = "final" ]
  then
   echo "
   INSERT OR REPLACE INTO metadata VALUES('minzoom','${min_z}');
   INSERT OR REPLACE INTO metadata VALUES('maxzoom','${max_z}');
   INSERT OR REPLACE INTO metadata VALUES('center','${center_lon},${center_lat},${max_z}');
   INSERT OR REPLACE INTO metadata VALUES('bounds','${min_x},${min_y},${max_x},${max_y}');
   END TRANSACTION;
   ANALYZE;
   VACUUM;
   " >> ${this_dir}.sql;
  else
   echo "END TRANSACTION; "  >> ${this_dir}.sql;
  fi
  if [ ! -f "${this_dir}.${MBTILE_EXTENTION}" ]
  then
   # Needed for the Android Application 'OruxMaps': 'CREATE TABLE android_metadata (locale text);'
   echo "
   BEGIN TRANSACTION;
   CREATE TABLE android_metadata (locale text);
   CREATE TABLE grid_key (grid_id TEXT,key_name TEXT);
   CREATE TABLE grid_utfgrid (grid_id TEXT,grid_utfgrid BLOB);
   CREATE TABLE keymap (key_name TEXT,key_json TEXT);
   CREATE TABLE images (tile_data blob,tile_id text);
   CREATE TABLE map (zoom_level INTEGER,tile_column INTEGER,tile_row INTEGER,tile_id TEXT,grid_id TEXT);
   CREATE TABLE metadata (name text,value text);
   CREATE VIEW tiles AS SELECT map.zoom_level AS zoom_level,map.tile_column AS tile_column,map.tile_row AS tile_row,images.tile_data AS tile_data FROM map JOIN images ON images.tile_id = map.tile_id ORDER BY zoom_level,tile_column,tile_row;
   CREATE VIEW grids AS SELECT map.zoom_level AS zoom_level,map.tile_column AS tile_column,map.tile_row AS tile_row,grid_utfgrid.grid_utfgrid AS grid FROM map JOIN grid_utfgrid ON grid_utfgrid.grid_id = map.grid_id;
   CREATE VIEW grid_data AS SELECT map.zoom_level AS zoom_level,map.tile_column AS tile_column,map.tile_row AS tile_row,keymap.key_name AS key_name,keymap.key_json AS key_json FROM map JOIN grid_key ON map.grid_id = grid_key.grid_id JOIN keymap ON grid_key.key_name = keymap.key_name;
   CREATE UNIQUE INDEX grid_key_lookup ON grid_key (grid_id,key_name);
   CREATE UNIQUE INDEX grid_utfgrid_lookup ON grid_utfgrid (grid_id);
   CREATE UNIQUE INDEX keymap_lookup ON keymap (key_name);
   CREATE UNIQUE INDEX images_id ON images (tile_id);
   CREATE UNIQUE INDEX map_index ON map (zoom_level, tile_column, tile_row);
   CREATE UNIQUE INDEX name ON metadata (name);
   INSERT INTO metadata VALUES('name','${this_dir}');
   INSERT INTO metadata VALUES('type','baselayer');
   INSERT INTO metadata VALUES('version','1.1');
   INSERT INTO metadata VALUES('description', '${this_dir} import from ${DIR_TYPE}-Tiles');
   INSERT INTO metadata VALUES('format','${TILE_EXTENTION_CONVERT}');
   INSERT INTO metadata VALUES('minzoom','${min_z}');
   INSERT INTO metadata VALUES('maxzoom','${max_z}');
   INSERT INTO metadata VALUES('center','${center_lon},${center_lat},${max_z}');
   INSERT INTO metadata VALUES('bounds','${min_x},${min_y},${max_x},${max_y}');
   INSERT INTO metadata VALUES('tile_row_type','tms');
   " >> "${this_dir}_create.sql"
   echo "-I-> Creating: ${this_dir}.${MBTILE_EXTENTION} ; bbox=${min_x},${min_y},${max_x},${max_y} min/max_zoom=${min_z},${max_z} ";
   cat ${this_dir}_create.sql ${this_dir}.sql | ${CMD_SQLITE} ${this_dir}.${MBTILE_EXTENTION}
   rm ${this_dir}_create.sql ${this_dir}.sql ;
  else
   echo "-I-> Inserting into: ${this_dir}.${MBTILE_EXTENTION} ; bbox=${min_x},${min_y},${max_x},${max_y} min/max_zoom=${min_z},${max_z} ";
   cat ${this_dir}.sql | ${CMD_SQLITE} ${this_dir}.${MBTILE_EXTENTION}
   local_rc=$?;
   if [ "${local_rc}" -eq 0 ] &&  [ -f "${this_dir}.${MBTILE_EXTENTION}" ]
   then
    # echo "-I-> Removing: ${this_dir}.sql ";
    rm  ${this_dir}.sql ;
   fi
  fi
 fi
}
# ------------------------------------
as_ewkt()
{
 zoom=$1;
 xtile=$2;
 ytile=$3;
 tms=$4;
 n=$(ytile2lat `expr ${ytile}` ${zoom} ${tms} );
 s=$(ytile2lat `expr ${ytile} + 1` ${zoom} ${tms} );
 e=$(xtile2long `expr ${xtile} + 1` ${zoom} ${tms} );
 w=$(xtile2long `expr ${xtile}` ${zoom} ${tms} );
 center_lat=`echo "${s} ${n}" | awk '{printf("%.8f", ($1 + $2) / 2.0)}'`
 center_lon=`echo "${w} ${e}" | awk '{printf("%.8f", ($1 + $2) / 2.0)}'`
 # SRID=4326;POINT(-83.3203125000000 30.9587685707799)
 # SRID=4326;POLYGON((-91.3842773437500 35.2231850497018, -75.2563476562500 35.2231850497018, -75.2563476562500 26.4951567924440, -91.3842773437500 26.4951567924440, -91.3842773437500 35.2231850497018))
 echo "SRID=4326;POINT(${center_lon} ${center_lat})\nSRID=4326;POLYGON((${w} ${n},${e} ${n},${e} ${s},${w} ${s},${w} ${n}))"
}
# ------------------------------------
as_gdalwarp()
{
 zoom=$1;
 xtile=$2;
 ytile=$3;
 tms=$4;
 if [ "${tms}" == "tms" ]
 then
  ytile_tms=${ytile};
  ytile_name="${ytile}_tms";
  ytile_osm=$( flip_tmsosm ${ytile} ${zoom} );
  ytile_name_to="${ytile_osm}_osm";
 else
  ytile_osm=${ytile};
  ytile_name="${ytile}_osm";
  ytile_tms=$( flip_tmsosm ${ytile} ${zoom} );
  ytile_name_to="${ytile_tms}_tms";
 fi
 OUTPUT_WGET="wget 'http://mt3.google.com/vt/lyrs=s,h&z=${zoom}&x=${xtile}&y=${ytile_osm}' -O ${zoom}_${xtile}_${ytile_name}.png\n";
 OUTPUT_LN="ln ${ytile}.png ${zoom}_${xtile}_${ytile_name}.png\n";
 n=$(ytile2lat `expr ${ytile}` ${zoom} ${tms} );
 s=$(ytile2lat `expr ${ytile} + 1` ${zoom} ${tms} );
 e=$(xtile2long `expr ${xtile} + 1` ${zoom} ${tms} );
 w=$(xtile2long `expr ${xtile}` ${zoom} ${tms} );
 OUTPUT_TRANSLATE="gdal_translate -gcp 0 $TILE_SIZE $w $s -gcp $TILE_SIZE 0 $e $n  -gcp 0 0 $w $n  -gcp $TILE_SIZE $TILE_SIZE  $e $s ${zoom}_${xtile}_${ytile_name}.png ${zoom}_${xtile}_${ytile_name}.gcp.tif\n";
 OUTPUT_WARP="gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3395  ${zoom}_${xtile}_${ytile_name}.gcp.tif ${zoom}_${xtile}_${ytile_name}.3395.tif \n";
 OUTPUT_INFO="gdalinfo ${zoom}_${xtile}_${ytile_name}.3395.tif  | grep Pixel \n";
 OUTPUT_CP="cp ${zoom}_${xtile}_${ytile_name}.3395.tif  ${zoom}_${xtile}_${ytile_name_to}.3395.tif \n";
 if [ "${zoom}" -gt "30" ]
 then
  zoom=30;
 fi
 POSITION_BASE="${1}_${xtile}_${ytile_name}.3395.tif: ${POSITION_BASE_ARRAY[$zoom]}";
 # Pixel Size = (1226.152784436568254,-1226.152784436568254) wsg84(0.011014717869251,-0.006857672434668)
 echo "${OUTPUT_WGET}${OUTPUT_LN}${OUTPUT_TRANSLATE}${OUTPUT_WARP}${OUTPUT_INFO}${OUTPUT_CP}rm ${zoom}_${xtile}_${ytile_name}.gcp.tif ${zoom}_${xtile}_${ytile_name}.png\n# ${POSITION_BASE}";
 # tileinfo 20 13.377914429 2.516220864
}
# ------------------------------------
as_position()
{
 zoom=$1;
 xtile=$2;
 ytile=$3;
 tms=$4;
 CENTER=$(as_center ${zoom} ${xtile} ${ytile} ${tms} );
 BBOX=$(as_bounding_box ${zoom} ${xtile} ${ytile} ${tms} );
 echo "${CENTER}\n${BBOX}"
}
# ------------------------------------
as_tile()
{
 zoom=$1;
 xtile=$2;
 ytile=$3;
 tms=$4;
 TILE_TMS="";
 TILE_OSM="";
 if [ -z "${tms}" ] ||  [ "${tms}" == "osm" ]
 then
  TILE_OSM="tile_osm : ${zoom}/${xtile}/${ytile}";
  ytile=$( flip_tmsosm ${ytile} ${zoom} );
  TILE_TMS="tile_tms : ${zoom}/${xtile}/${ytile}";
 fi
 if [ "${tms}" == "tms" ]
 then
  TILE_TMS="tile_tms : ${zoom}/${xtile}/${ytile}";
  ytile=$( flip_tmsosm ${ytile} ${zoom} );
  TILE_OSM="tile_osm : ${zoom}/${xtile}/${ytile}";
 fi
 echo -e "${TILE_OSM}\n${TILE_TMS}";
}
# ------------------------------------
On_Help()
{
 T="     ";
 POS_TEXT="${T}${T}";
 COMMENT_TEXT="";
 if [ ! -z "${CMD_SPATIALITE}" ]
 then
  COMMENT_TEXT="${COMMENT_TEXT} Spatialite is installed.";
 else
  COMMENT_TEXT="${COMMENT_TEXT} Spatialite is not installed.";
 fi
 if [ ! -z "${CMD_SQLITE}" ]
 then
  COMMENT_TEXT="Sqlite is installed.";
 else
  COMMENT_TEXT="Sqlite is not installed.";
 fi
 #------------------------------------------------------------------------------
 RC_TEXT="$RC_TEXT-H-> Usage: ${BASE_NAME} zoom-level x-position y-position [tms]";
 RC_TEXT="$RC_TEXT\n${POS_TEXT}: if 'x/y-position' are double numbers\t: Convert to tiles using the zoom-level";
 RC_TEXT="$RC_TEXT\n${POS_TEXT}: if 'x/y-position' are integers numbers\t: Convert to Latitude/Longitude using the zoom-level";
 RC_TEXT="$RC_TEXT\n${POS_TEXT}: if [tms] is set the tile y-position will be assumed as in the tms-format, otherwise as the osm-format";
 RC_TEXT="$RC_TEXT\n-H->${POS_TEXT}Sample:${T}Lat/Long to Tile:${T}${T}${BASE_NAME} ${ZOOM} ${LONG} ${LAT}";
 RC_TEXT="$RC_TEXT\n-H->${POS_TEXT}Sample:${T}Osm-Tile to Lat/Long:${T}${BASE_NAME} ${ZOOM} ${TILE_X} ${TILE_Y_OSM}";
 RC_TEXT="$RC_TEXT\n-H->${POS_TEXT}Sample:${T}Tms-Tile to Lat/Long:${T}${BASE_NAME} ${ZOOM} ${TILE_X} ${TILE_Y_TMS} tms";
 RC_TEXT="$RC_TEXT\n-H->${POS_TEXT}or:${T}Osm-Tile to Lat/Long:${T}${BASE_NAME} /${ZOOM}/${TILE_X}/${TILE_Y_OSM}";
 RC_TEXT="$RC_TEXT\n-H->${POS_TEXT}or:${T}Tms-Tile to Lat/Long:${T}${BASE_NAME} ${ZOOM}/${TILE_X}/${TILE_Y_TMS} tms";
 RC_TEXT="$RC_TEXT\n-H->${POS_TEXT}or:${T}Tms-Tile to Lat/Long:${T}${BASE_NAME} /${ZOOM}/${TILE_X}/${TILE_Y_TMS}.extention tms";
 RC_TEXT="$RC_TEXT\n-H->${POS_TEXT}Sample:${T}Test-Position:${T}${BASE_NAME} -test [will show position of the Brandenburg Gate, Berlin Germany at zoom 17]";
 RC_TEXT="$RC_TEXT\n-H->${POS_TEXT}Sample:${T}Zoom to Lat/Long:${T}${BASE_NAME} ${ZOOM} [will show position of the Brandenburg Gate, Berlin Germany at the given zoom-level]";
 RC_TEXT="$RC_TEXT\n-H->${POS_TEXT}Note:${T}Sample gdal-commands will be shown that can be used to georeference a tile to : 'WGS 84 / World Mercator' srid:3395";
 RC_TEXT="$RC_TEXT\n${POS_TEXT}---> creating [${ZOOM}_${TILE_X}_${TILE_Y_OSM}_osm.3395.tif, ${ZOOM}_${TILE_X}_${TILE_Y_TMS}_tms.3395.tif]";
 RC_TEXT="$RC_TEXT\nOR to create a MBTiles-Database ";
 RC_TEXT="$RC_TEXT\n${POS_TEXT}---> [call from the main Directory of the tiles - where the sub-Directories are called '0,1,..,19']";
 RC_TEXT="$RC_TEXT\n-H->${POS_TEXT}MBTiles:${T}Create from Osm Tiles-Diretory:${T}${BASE_NAME} -mbtiles";
 RC_TEXT="$RC_TEXT\n-H->${POS_TEXT}MBTiles:${T}Create from Tms Tiles-Diretory:${T}${BASE_NAME} -mbtiles tms";
 RC_TEXT="$RC_TEXT\nALSO possible\n-H->${POS_TEXT}MBTiles:${T}Create using spatialite instead of sqlite3:${T}${BASE_NAME} -mbtiles -spatialite [tms] ";
 RC_TEXT="$RC_TEXT\n${POS_TEXT}---> the created Database will be stored in the main Directory of the tiles.";
 RC_TEXT="$RC_TEXT\n${POS_TEXT}-----------------------------------------------------------------------------------------";
 RC_TEXT="$RC_TEXT\n${POS_TEXT}-> Problem: gdal2tiles creates the 'mercator' tiles (for osm) with tms-numbering";
 RC_TEXT="$RC_TEXT\n${POS_TEXT}-> http://tile.openstreetmap.org/1/0/0.png = North America [osm]  ; South America [tms]";
 RC_TEXT="$RC_TEXT\n${POS_TEXT}-> http://tile.openstreetmap.org/1/0/1.png =  South America [osm]  ; North America [tms]";
 RC_TEXT="$RC_TEXT\n${POS_TEXT}-> http://tile.openstreetmap.org/1/1/0.png =  Northern Asia [osm]  ; Southern Asia [tms]";
 RC_TEXT="$RC_TEXT\n${POS_TEXT}-> http://tile.openstreetmap.org/1/1/1.png =  Southern Asia [osm]   ; Northern Asia [tms]";
 RC_TEXT="$RC_TEXT\n${POS_TEXT}-----------------------------------------------------------------------------------------";
 RC_TEXT="$RC_TEXT\n${POS_TEXT}---> the .mbtiles Databases allways uses the tms-numbering for the y-position. .png files will be converted to jpg";
 RC_TEXT="$RC_TEXT\n${POS_TEXT}---> ${COMMENT_TEXT}";
 #------------------------------------------------------------------------------
 echo ${RC_TEXT};
}
# ------------------------------------
if [ "${MODUS}"  -gt "0" ]
then
 echo "${MESSAGE_TYPE} ${INPUT_PARMS}";
 if [ "${MODUS}" -eq "3" ]
 then
  find_image_type;
  local_rc=$?;
  if [ "${local_rc}" -eq 0 ]
  then
   loop_tile_dir "${PWD##*/}" "${TMS}";
  fi
 else
  if [ "${MODUS}" -eq "1" ]
  then
   TILE_X=$(long2xtile ${LONG} ${ZOOM});
   # echo -e "called long2xtile [${LONG} ${ZOOM}] result[${TILE_X}]";
   TILE_Y=$(lat2ytile ${LAT} ${ZOOM} ${TMS});
   echo -e "called lat2ytile [${LAT} ${ZOOM} ${TMS}] result[${TILE_Y}]";
   # long2_geodetic_xtile()
   # lat2_geodetic_ytile()
   TILE_X_GEOD=$( long2_geodetic_xtile ${LONG} ${ZOOM} )
   TILE_Y_GEOD=$( lat2_geodetic_ytile ${LAT} ${ZOOM} )
  fi
  if [ "${MODUS}" -eq "2" ]
  then
   LONG=$(xtile2long ${TILE_X} ${ZOOM});
   # echo -e "called xtile2long [${TILE_X} ${ZOOM}] result[${LONG}]";
   LAT=$(ytile2lat ${TILE_Y} ${ZOOM} ${TMS});
  fi
  TILE=$(as_tile ${ZOOM} ${TILE_X} ${TILE_Y} ${TMS});
  echo -e "${TILE}";
  POSITION=$(as_position ${ZOOM} ${TILE_X} ${TILE_Y} ${TMS});
  echo -e "${POSITION}";
  EWKT=$(as_ewkt ${ZOOM} ${TILE_X} ${TILE_Y} ${TMS});
  echo -e "${EWKT}";
  GDAL=$(as_gdalwarp ${ZOOM} ${TILE_X} ${TILE_Y} ${TMS});
  echo -e "${GDAL}";
 fi
else
 # ------------------------------------
 RC_TEXT=$( On_Help );
 echo -e "${RC_TEXT}";
 MESSAGE_TYPE="-H->";
fi
#---------------------------------------------------
if [ "$exit_rc" -eq "0" ]
then
 RC_TEXT=$BENE;
else
 RC_TEXT="$BENENO";
 MESSAGE_TYPE="-E->";
fi
echo "${MESSAGE_TYPE} ${BASE_NAME_CUT} rc=$exit_rc [${RC_TEXT}] - ${HABE_FERTIG}";
exit $exit_rc;
#---------------------------------------------------
tileinfo 17 70406 42988
-I-> Input: Osm-Tile z/x/y[17/70406/42988] tms[88083]
tile_osm : 17/70406/42988
tile_tms : 17/70406/88083
center=13.37722778,52.51538515
bbox=13.375854492,52.514549436,13.378601074,52.516220864
SRID=4326;POINT(13.37722778 52.51538515)
SRID=4326;POLYGON((13.375854492 52.516220864,13.378601074 52.516220864,13.378601074 52.514549436,13.375854492 52.514549436,13.375854492 52.516220864))
wget 'http://mt3.google.com/vt/lyrs=s,h&z=17&x=70406&y=42988' -O 17_70406_42988_osm.png
ln 42988.png 17_70406_42988_osm.png
gdal_translate -gcp 0 256 13.375854492 52.514549436 -gcp 256 0 13.378601074 52.516220864  -gcp 0 0 13.375854492 52.516220864  -gcp 256 256  13.378601074 52.514549436 17_70406_42988_osm.png 17_70406_42988_osm.gcp.tif
gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3395  17_70406_42988_osm.gcp.tif 17_70406_42988_osm.3395.tif
gdalinfo 17_70406_42988_osm.3395.tif  | grep Pixel
cp 17_70406_42988_osm.3395.tif  17_70406_88083_tms.3395.tif
rm 17_70406_42988_osm.gcp.tif 17_70406_42988_osm.png
# 17_70406_42988.3395.tif: Pixel Size = (1.192842756645181,-1.192842756645181)
#---------------------------------------------------
tileinfo 17 70406 88083 tms
-I-> Input: Tms-Tile z/x/y[17/70406/88083] osm[42988]
tile_osm : 17/70406/42988
tile_tms : 17/70406/88083
center=13.37722778,52.51705655
bbox=13.375854492,52.517892228,13.378601074,52.516220864
SRID=4326;POINT(13.37722778 52.51705655)
SRID=4326;POLYGON((13.375854492 52.516220864,13.378601074 52.516220864,13.378601074 52.517892228,13.375854492 52.517892228,13.375854492 52.516220864))
wget 'http://mt3.google.com/vt/lyrs=s,h&z=17&x=70406&y=42988' -O 17_70406_88083_tms.png
ln 88083.png 17_70406_88083_tms.png
gdal_translate -gcp 0 256 13.375854492 52.517892228 -gcp 256 0 13.378601074 52.516220864  -gcp 0 0 13.375854492 52.516220864  -gcp 256 256  13.378601074 52.517892228 17_70406_88083_tms.png 17_70406_88083_tms.gcp.tif
gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3395  17_70406_88083_tms.gcp.tif 17_70406_88083_tms.3395.tif
gdalinfo 17_70406_88083_tms.3395.tif  | grep Pixel
cp 17_70406_88083_tms.3395.tif  17_70406_42988_osm.3395.tif
rm 17_70406_42988_osm.gcp.tif 17_70406_42988_osm.png
# 17_70406_88083.3395.tif: Pixel Size = (1.192842756645181,-1.192842756645181)
 
#---------------------------------------------------
tileinfo 17 13.37771496361961 52.51628011262304
-I-> Input: Position z/x,y[17/13.37771496361961,52.51628011262304] tms[osm] modus[1]
called long2xtile [13.37771496361961 17] result[70406]
tile_osm : 17/70406/42988
tile_tms : 17/70406/88083
center=13.37722778,52.51538515
bbox=13.375854492,52.514549436,13.378601074,52.516220864
SRID=4326;POINT(13.37722778 52.51538515)
SRID=4326;POLYGON((13.375854492 52.516220864,13.378601074 52.516220864,13.378601074 52.514549436,13.375854492 52.514549436,13.375854492 52.516220864))
wget 'http://mt3.google.com/vt/lyrs=s,h&z=17&x=70406&y=42988' -O 17_70406_42988_osm.png
ln 42988.png 17_70406_42988_osm.png
gdal_translate -gcp 0 256 13.375854492 52.514549436 -gcp 256 0 13.378601074 52.516220864  -gcp 0 0 13.375854492 52.516220864  -gcp 256 256  13.378601074 52.514549436 17_70406_42988_osm.png 17_70406_42988_osm.gcp.tif
gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3395  17_70406_42988_osm.gcp.tif 17_70406_42988_osm.3395.tif
gdalinfo 17_70406_42988_osm.3395.tif  | grep Pixel
cp 17_70406_42988_osm.3395.tif  17_70406_88083_tms.3395.tif
rm 17_70406_42988_osm.gcp.tif 17_70406_42988_osm.png
# 17_70406_42988.3395.tif: Pixel Size = (1.192842756645181,-1.192842756645181)
#---------------------------------------------------

creating a mbtiles-file with tileinfo.sh

Creating a mbtiles-file with tileinfo.sh is relativly easy:

  • assuming tileinfo.sh has been copied somewhere were it can be found like '/usr/local/bin' and is excutable
  • assuming a tile directory exists with sub-directories called '1,2,3 .. 18,19'
  • assuming you know how the y-numbering was done ['osm' - North to South or 'tms' - South to North]
    • open a terminal in this directory, with 'ls' you should see something like this: '12 17 18'
  • file extensions supported: jpg,jpeg,png
    • whereby png will be converted to jpg if imagemagik convert is installed

---

Senario 1:

  • directory: 'Berlin_Area'
    • 'osm' was used to create these tiles
tileinfo.sh -mbtiles

will, when finished, have a Berlin_Area.mbtiles file in the Berlin_Area directory

---

Senario 2:

  • directory: 'Schoeneberg_Verwaltungsbezirk_1931'
    • 'tms' which was created with 'gdal2tiles.py' using --profile=mercator
tileinfo.sh -mbtiles tms

will, when finished, have a Schoeneberg_Verwaltungsbezirk_1931.mbtiles file in the Schoeneberg_Verwaltungsbezirk_1931 directory

---

Files created with this script have been tested with the android application 'geopaparazzi'

  • care must be taken to avoid mbtiles files that are larger than 2.1 GB
    • many of the sdcard's on android cannot store such files

Crude map scale

Beware: True scale value will depend on your specific monitor size and current resolution. This is just a guide, map scale is only really meaningful for the printed page. This dpi to meters-on-the-ground "pixel factor" is based on someone's ancient CRT monitor, perhaps a 15" model at 800x600 resolution.

calc_webtile_scale()
{
    lat=$1
    scale=$2
    echo "$lat $scale" | awk -v PI=3.14159265358979323846 \
      -v PIXELFACT=2817.947378 '{
       a = 6378137.0;
       printf("%.0f", (a * 2*PI * cos($1 * PI/180) * PIXELFACT) / (256 * 2^$2))
 
    }'
}
 
mapscale=$(calc_webtile_scale $center_lat $zoom)

Octave

Lon./lat. to tile numbers

% convert the degrees to radians
rho = pi/180;
lon_rad = lon_deg * rho;
lat_rad = lat_deg * rho;
 
n = 2 ^ zoom
xtile = n * ((lon_deg + 180) / 360)
ytile = n * (1 - (log(tan(lat_rad) + sec(lat_rad)) / pi)) / 2

Tile numbers to lon./lat.

n=2^zoom
lon_deg = xtile / n * 360.0 - 180.0
lat_rad = arctan(sinh(pi * (1 - 2 * ytile / n)))
lat_deg = lat_rad * 180.0 / pi


Emacs-lisp

(defun longitude2tile (lon zoom) (* (expt 2 zoom) (/ (+ lon 180) 360)))
 
(defun tile2longitude (x zoom) (- (/ (* x 360) (expt 2 zoom)) 180))
 
(defun latitude2tile (lat zoom) (* (expt 2 zoom) (/ (- 1 (/ (log (+ (tan (/ (* lat pi) 180)) (/ 1 (cos (/ (* lat pi) 180))))) pi)) 2)))
 
(defun sinh (value) (/ (- (exp value) (exp (- value))) 2))
(defun tile2latitude (y zoom) (/ (* 180 (atan (sinh (* pi (- 1 (* 2 (/ y (expt 2 zoom)))))))) pi))

Erlang

-module(slippymap).
-export([deg2num/3]).
-export([num2deg/3]).
 
deg2num(Lat,Lon,Zoom)->
    X=math:pow(2, Zoom) * ((Lon + 180) / 360),
    Sec=1/math:cos(deg2rad(Lat)),
    R = math:log(math:tan(deg2rad(Lat)) + Sec)/math:pi(),
    Y=math:pow(2, Zoom) * (1 - R) / 2,
    {round(X),round(Y)}.
 
num2deg(X,Y,Zoom)->
    N=math:pow(2, Zoom),
    Lon=X/N*360-180,
    Lat_rad=math:atan(math:sinh(math:pi()*(1-2*Y/N))),
    Lat=Lat_rad*180/math:pi(),
    {Lon,Lat}.
 
deg2rad(C)->
    C*math:pi()/180.

Lua

function deg2num(lon, lat, zoom)
    local n = 2 ^ zoom
    local lon_deg = tonumber(lon)
    local lat_rad = math.rad(lat)
    local xtile = math.floor(n * ((lon_deg + 180) / 360))
    local ytile = math.floor(n * (1 - (math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi)) / 2)
    return xtile, ytile
end
 
function num2deg(x, y, z)
    local n = 2 ^ z
    local lon_deg = x / n * 360.0 - 180.0
    local lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * y / n)))
    local lat_deg = lat_rad * 180.0 / math.pi
    return lon_deg, lat_deg
end

PostgreSQL

CREATE OR REPLACE FUNCTION lon2tile(lon double precision, zoom integer)
  RETURNS integer AS
$BODY$
BEGIN
    return floor( (lon + 180) / 360 * (1 << zoom) )::integer;
END;
$BODY$
  LANGUAGE plpgsql IMMUTABLE;
 
 
CREATE OR REPLACE FUNCTION lat2tile(lat double precision, zoom integer)
  RETURNS integer AS
$BODY$
BEGIN
    return floor( (1.0 - ln(tan(radians(lat)) + 1.0 / cos(radians(lat))) / pi()) / 2.0 * (1 << zoom) )::integer;
END;
$BODY$
  LANGUAGE plpgsql IMMUTABLE;
 
 
CREATE OR REPLACE FUNCTION tile2lat(y integer, zoom integer)
  RETURNS double precision AS
$BODY$
DECLARE
 n float;
 sinh float;
 E float = 2.7182818284;
BEGIN
    n = pi() - (2.0 * pi() * y) / power(2.0, zoom);
    sinh = (1 - power(E, -2*n)) / (2 * power(E, -n));
    return degrees(atan(sinh));
END;
$BODY$
  LANGUAGE plpgsql IMMUTABLE;
 
 
CREATE OR REPLACE FUNCTION tile2lon(x integer, zoom integer)
  RETURNS double precision AS
$BODY$
BEGIN
 return x * 1.0 / (1 << zoom) * 360.0 - 180.0;
END;
$BODY$
  LANGUAGE plpgsql IMMUTABLE;

Subtiles

If you're looking at tile x,y and want to zoom in, the subtiles are (in the next zoom-level's coordinate system):

2x, 2y 2x + 1, 2y
2x, 2y + 1 2x + 1, 2y + 1

Similarly, zoom out by halving x and y (in the previous zoom level)

Resolution and Scale

Exact length of the equator (according to wikipedia) is 40075.016686 km in WGS-84. A horizontal tile size at zoom 0 would be 156543.034 meters:

6378137.0 * 2 * pi / 256 = 156543.034

Which gives us a formula to calculate resolution at any given zoom:

resolution = 156543.034 meters/pixel * cos(latitude) / (2 ^ zoomlevel)

Some applications need to know a map scale, that is, how 1 cm on a screen translates to 1 cm of a map.

scale = 1 : (screen_dpi * 39.37 in/m * resolution)

And here is the table to rid you of those calculations. All values are shown for equator, and you have to multiply them by cos(latitude) to adjust to a given latitude. For example, divide those by 2 for latitude 60 (Oslo, Helsinki, Saint-Petersburg).

zoom resolution, m/px scale 96 dpi 1 screen cm is scale 120 dpi
0 156543.03 1 : 554 678 932 5547 km 1 : 739 571 909
1 78271.52 1 : 277 339 466 2773 km 1 : 369 785 954
2 39135.76 1 : 138 669 733 1337 km 1 : 184 892 977
3 19567.88 1 : 69 334 866 693 km 1 : 92 446 488
4 9783.94 1 : 34 667 433 347 km 1 : 46 223 244
5 4891.97 1 : 17 333 716 173 km 1 : 23 111 622
6 2445.98 1 : 8 666 858 86.7 km 1 : 11 555 811
7 1222.99 1 : 4 333 429 43.3 km 1 : 5 777 905
8 611.50 1 : 2 166 714 21.7 km 1 : 2 888 952
9 305.75 1 : 1 083 357 10.8 km 1 : 1 444 476
10 152.87 1 : 541 678 5.4 km 1 : 722 238
11 76.437 1 : 270 839 2.7 km 1 : 361 119
12 38.219 1 : 135 419 1.4 km 1 : 180 559
13 19.109 1 : 67 709 677 m 1 : 90 279
14 9.5546 1 : 33 854 339 m 1 : 45 139
15 4.7773 1 : 16 927 169 m 1 : 22 569
16 2.3887 1 : 8 463 84.6 m 1 : 11 284
17 1.1943 1 : 4 231 42.3 m 1 : 5 642
18 0.5972 1 : 2 115 21.2 m 1 : 2 821

See also Zoom levels

Tools

References

(note: Slippy tiles and Google map tiles count tile 0,0 down from the top-left of the tile grid; the TMS spec specifies tiles count up from 0,0 in the lower-left!)