"""
Conversions between WGS84, Cartesian, Mercator and TMS coordinates and bounding boxes,
useful for converting geographic view ports as generated by Google Earth to TMS bounding boxes
"""
import math
import itertools
from abc import ABCMeta, abstractmethod
__author__ = "Martin Loetzsch, Gregor Sturm"
__licence__ = "Apache 2.0"
_earthradius = 6378137.0 # m
_tilesize = _initialresolution = _originshift = None
_originshift_180 = _180_originshift = None
_pi_2 = math.pi / 2.0
_pi_180 = math.pi / 180.0
_pi_360 = math.pi / 360.0
_180_pi = 180.0 / math.pi
[docs]def init_geometry(tilesize=256.0):
global _tilesize, _initialresolution, _originshift
global _originshift_180, _180_originshift
_tilesize = tilesize
_initialresolution = 2 * math.pi * _earthradius / _tilesize
_originshift = 2 * math.pi * _earthradius / 2.0
_originshift_180 = _originshift / 180.0
_180_originshift = 180.0 / _originshift
init_geometry()
[docs]def griditer(x, y, ncol, nrow=None, step=1):
"""
Iterate through a grid of tiles.
Args:
x (int): x start-coordinate
y (int): y start-coordinate
ncol (int): number of tile columns
nrow (int): number of tile rows. If not specified, this
defaults to ncol, s.t. a quadratic region is
generated
step (int): clear. Analogous to range().
Yields:
Tuple: all tuples (x, y) in the region delimited by
(x, y), (x + ncol, y + ncol).
"""
if nrow is None:
nrow = ncol
yield from itertools.product(range(x, x + ncol, step),
range(y, y + nrow, step))
[docs]def bboxiter(tile_bounds, tiles_per_row_per_region=1):
"""
Iterate through a grid of regions defined by a TileBB.
Args:
tile_bounds (GridBB):
tiles_per_row_per_region: Combine multiple tiles in one region.
E.g. if set to two, four tiles will be combined in one region.
See `kml` module description for more details. Leaving the
default value '1' simply yields all tiles in the bounding box.
Note:
If the number of regions would not be an integer due to specification
of the `tiles_per_row_per_region` parameter, the boundaries will
be rounded to the next smaller or next larger integer respectively.
Example:
We have the following bounding box with size 2x2 and set
tiles_per_row_per_region = 2, delimited by the coordinates (x, y)::
(5,5)--- ---
| |
| |
| |
--- ---(9,7)
Although this could be represented in one single region with two
tiles per row, it will create four regions::
(2,2)--- --- (5/2 = 2.5 -> 2, 5/2 = 2.5 -> 2)
| | |
--- ---
| | |
--- ---(5,4) (9/2 = 4.5 -> 5, 7/2 = 3.5 -> 4)
Yields:
Tuple: all tuples (x, y) in the region delimited by the TileBB
"""
x_lower = math.floor(tile_bounds.min.x / tiles_per_row_per_region)
y_lower = math.floor(tile_bounds.min.y / tiles_per_row_per_region)
ncol = math.ceil(tile_bounds.max.x / tiles_per_row_per_region) - x_lower
nrow = math.ceil(tile_bounds.max.y / tiles_per_row_per_region) - y_lower
yield from griditer(x_lower, y_lower, ncol, nrow)
[docs]class GeographicCoordinate:
"""
Represents a WGS84 Datum
Args:
lon: longitude in degrees
lat: latitude in degrees
height: height in meters above the surface of the earth spheroid
"""
def __init__(self, lon=None, lat=None, height=0.0):
self.lon = lon
self.lat = lat
self.height = height
[docs] def to_cartesian(self):
r = (_earthradius + self.height)
cosx = math.cos(_pi_180 * self.lon)
cosy = math.cos(_pi_180 * self.lat)
sinx = math.sin(_pi_180 * self.lon)
siny = math.sin(_pi_180 * self.lat)
return CartesianCoordinate(r * cosy * cosx, r * cosy * sinx, r * siny)
[docs] def to_mercator(self):
return MercatorCoordinate(
self.lon * _originshift_180,
math.log(math.tan((90 + self.lat) * _pi_360)) / _pi_180 * _originshift_180)
def __str__(self):
return "<lon: " + str(self.lon) + ", lat: " + str(self.lat) + ", height: " \
+ str(self.height) + ">"
[docs]class GeographicBB:
""" A bounding box defined by two geographic coordinates """
def __init__(self, min_lon=None, min_lat=None, max_lon=None, max_lat=None):
self.min = GeographicCoordinate(min_lon, min_lat)
self.max = GeographicCoordinate(max_lon, max_lat)
[docs] def intersects(self, other):
return self.min.lon <= other.max.lon and self.max.lon >= other.min.lon and \
self.min.lat <= other.max.lat and self.max.lat >= other.min.lat
[docs] def intersection(self, other):
intersects = self.intersects(other)
if intersects:
return GeographicBB(max(self.min.lon, other.min.lon),
max(self.min.lat, other.min.lat),
min(self.max.lon, other.max.lon),
min(self.max.lat, other.max.lat))
elif self.max.lon > 180.0:
return GeographicBB(-180, self.min.lat,
self.max.lon - 360.0, self.max.lat).intersection(other)
[docs] def center(self):
return GeographicCoordinate(self.min.lon + (self.max.lon - self.min.lon) / 2,
self.min.lat + (self.max.lat - self.min.lat) / 2)
[docs] def to_mercator(self):
return MercatorBB(self.min.to_mercator(), self.max.to_mercator())
def __str__(self):
return "<geographic min: " + str(self.min) + ", max: " + str(self.max) + ">"
[docs]class CartesianCoordinate:
""" Represents a coordinate in a geocentric Cartesian coordinate system """
def __init__(self, x, y, z):
self.x = x
self.y = y
self.z = z
def __sub__(self, other):
return CartesianCoordinate(self.x - other.x, self.y - other.y, self.z - other.z)
[docs] def length(self):
return math.sqrt(self.x ** 2 + self.y ** 2 + self.z ** 2)
def __str__(self):
return "<x: " + str(self.x) + ", y: " + str(self.y) + ", z: " + str(self.z) + ">"
[docs]class MercatorCoordinate:
""" Represents a coordinate in Spherical Mercator EPSG:900913 """
def __init__(self, x=None, y=None):
self.x = x
self.y = y
[docs] def to_tile(self, zoom):
res = _initialresolution / (2 ** zoom)
def transform(x):
return int(math.ceil(((x + _originshift) / res) / _tilesize) - 1)
return TileCoordinate(zoom, transform(self.x),
(2 ** zoom - 1) - transform(self.y))
[docs] def to_geographic(self):
return GeographicCoordinate(
None if self.x is None else self.x * _180_originshift,
None if self.y is None
else (2 * math.atan(math.exp(self.y * _180_originshift * _pi_180)) - _pi_2) \
* _180_pi)
def __str__(self):
return "<x: " + str(self.x) + ", y: " + str(self.y) + ">"
[docs]class MercatorBB:
""" A bounding box defined by two mercator coordinates """
def __init__(self, _min=MercatorCoordinate, _max=MercatorCoordinate):
self.min = _min
self.max = _max
[docs] def to_tile(self, zoom):
""" Converts EPSG:900913 to tile coordinates in given zoom level """
p1 = self.min.to_tile(zoom)
p2 = self.max.to_tile(zoom)
return GridBB(zoom, p1.x, p2.y, p2.x, p1.y)
def __str__(self):
return "<mercator min: " + str(self.min) + ", max: " + str(self.max) + ">"
[docs]class GridCoordinate(metaclass=ABCMeta):
"""
Represents a position in a worldwide grid.
"""
def __init__(self, zoom, x, y):
self.zoom = zoom
self.x = x
self.y = y
[docs] def zoom_in(self):
"""
Yields:
GridCoordinate: the four tiles of the next zoom level
"""
return
[docs] def encode_quad_tree(self):
quad_key = ""
tx = self.x
ty = self.y
for i in range(self.zoom, 0, -1):
digit = 0
mask = 1 << (i - 1)
if (tx & mask) != 0:
digit += 1
if (ty & mask) != 0:
digit += 2
quad_key += str(digit)
return quad_key
def __str__(self):
return "<zoom: " + str(self.zoom) + ", x: " + str(self.x) + ", y: " + str(self.y) + ">"
[docs]class TileCoordinate(GridCoordinate):
""" Represents a coordinate in a worldwide tile grid. (aka the google-system)"""
def __init__(self, zoom, x, y):
super(TileCoordinate, self).__init__(zoom, x, y)
[docs] def to_mercator(self):
res = _initialresolution / (2 ** self.zoom)
return MercatorCoordinate(
None if self.x is None else self.x * _tilesize * res - _originshift,
None if self.y is None else (2 ** self.zoom - self.y) * _tilesize * res - _originshift)
[docs] def to_geographic(self):
return self.to_mercator().to_geographic()
[docs] def geographic_bounds(self):
p1 = self.to_geographic()
p2 = TileCoordinate(self.zoom, self.x + 1, self.y + 1).to_geographic()
return GeographicBB(p1.lon, p2.lat, p2.lon, p1.lat)
[docs] def zoom_in(self):
for x, y in griditer(self.x * 2, self.y * 2, ncol=2):
yield TileCoordinate(self.zoom + 1, x, y)
[docs] def resolution(self):
"""
Get the tile resolution at the current position.
The scale in WG84 depends on
* the zoom level (obviously)
* the latitude
* the tile size
References:
* http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Resolution_and_Scale
* http://gis.stackexchange.com/questions/7430/what-ratio-scales-do-google-maps-zoom-levels-correspond-to
Returns:
float: meters per pixel
"""
geo_coords = self.to_geographic()
resolution = abs(_initialresolution * math.cos(geo_coords.lat * _pi_180) / (2**self.zoom))
return resolution
[docs]class RegionCoordinate(GridCoordinate):
"""
Represents a region containing multiple tiles in a worldwide grid of regions.
A region is a square containing multiple tiles, e.g.::
-- --
| 1| 2|
-- --
| 3| 4|
-- --
A region must contain at least one tile and each row must have
a power of two of tiles. The size of the region is specified
with log2(tiles per row per region), e.g.
* log_tiles_per_row = 0 means 2**0 = 1 tile
* log_tiles_per_row = 2 means 2**2 = 4 tiles per row, thus 16 tiles per region.
Args:
zoom (int): clear.
x (int): clear.
y (int): clear.
log_tiles_per_row (int): size of the region as log2(tiles per row per region).
needs to be at least 0 (-> 1 tile) and at most 4 (-> 256 tiles)
"""
def __init__(self, zoom, x, y, log_tiles_per_row=0):
assert log_tiles_per_row in range(0, 5)
super().__init__(zoom, x, y)
self.log_tiles_per_row = log_tiles_per_row
self.tiles_per_row = 2 ** log_tiles_per_row
self.root_tile = TileCoordinate(zoom, x * self.tiles_per_row, y * self.tiles_per_row)
[docs] def geographic_bounds(self):
p1 = self.root_tile.to_geographic()
p2 = TileCoordinate(self.root_tile.zoom,
self.root_tile.x + self.tiles_per_row,
self.root_tile.y + self.tiles_per_row).to_geographic()
return GeographicBB(p1.lon, p2.lat, p2.lon, p1.lat)
[docs] def get_tiles(self):
"""Get all TileCoordinates contained in the region"""
for x, y in griditer(self.root_tile.x, self.root_tile.y, ncol=self.tiles_per_row):
yield TileCoordinate(self.root_tile.zoom, x, y)
[docs] def zoom_in(self):
for x, y in griditer(self.x * 2, self.y * 2, ncol=2):
yield RegionCoordinate(self.zoom + 1, x, y, log_tiles_per_row=self.log_tiles_per_row)
def __str__(self):
return ("<zoom: " + str(self.zoom) + ", x: " + str(self.x) + ", y: " + str(self.y) +
", log_tiles_per_row: " + str(self.log_tiles_per_row) + ">")
[docs]class GridBB:
""" A bounding box defined by two grid coordinates """
def __init__(self, zoom, min_x, min_y, max_x, max_y):
self.zoom = zoom
self.min = GridCoordinate(zoom, min_x, min_y)
self.max = GridCoordinate(zoom, max_x, max_y)
[docs] def intersection(self, other):
intersects = self.min.x <= other.max.x and self.max.x >= other.min.x and \
self.min.y <= other.max.y and self.max.y >= other.min.y
if intersects:
return GridBB(self.zoom,
max(self.min.x, other.min.x),
max(self.min.y, other.min.y),
min(self.max.x, other.max.x),
min(self.max.y, other.max.y))
else:
return []
[docs] def intersections(self, other):
intersections = []
i = self.intersection(other)
if i: intersections.append(i)
if self.max.x >= 2 ** self.zoom:
i = GridBB(self.zoom, 0, self.min.y, self.max.x % (2 ** self.zoom), self.max.y).intersection(other)
if i: intersections.append(i)
return intersections
[docs] def is_inside(self, tile):
return self.min.y <= tile.y <= self.max.y and \
(self.min.x <= tile.x <= self.max.x or
(self.max.x > 2 ** tile.zoom - 1 and
0 <= tile.x <= self.max.x % (2 ** tile.zoom)))
def __str__(self):
return "<tile min: " + str(self.min) + ", max: " + str(self.max) + ">"