diff --git a/landez/proj.py b/landez/proj.py index 1b3d739..ae03b6b 100644 --- a/landez/proj.py +++ b/landez/proj.py @@ -124,3 +124,211 @@ def tileslist(self, bbox): continue l.append((z, x, y)) return l + + +class _OSRTransformer(object): + """ Utility class for converting coordinates based on OSR module """ + def __init__(self, proj_name): + self.lonlat = osr.SpatialReference() + self.lonlat.ImportFromProj4('+init=epsg:4326') + self.proj = osr.SpatialReference() + self.proj.ImportFromProj4('+init=' + proj_name.lower()) + self.from_ll_ct = osr.CoordinateTransformation(self.lonlat, self.proj) + self.to_ll_ct = osr.CoordinateTransformation(self.proj, self.lonlat) + + def from_lonlat(self, lon, lat): + return self.from_ll_ct.TransformPoint(lon, lat) + + def to_lonlat(self, x, y): + return self.to_ll_ct.TransformPoint(x, y) + + +class _PyProjTransformer(object): + """ Utility class for converting coordinates based on PyProj module """ + def __init__(self, proj_name): + self.lonlat = pyproj.Proj(init='EPSG:4326') + self.proj = pyproj.Proj(init=proj_name) + + def from_lonlat(self, lon, lat): + return pyproj.transform(self.lonlat, self.proj, lon, lat) + + def to_lonlat(self, x, y): + return pyproj.transform(self.proj, self.lonlat, x, y) + + +class _DumbTransformer(object): + """ + Fallback class for converting coordinates when no proj module is available + """ + def from_lonlat(self, lon, lat): + return lon, lat + + def to_lonlat(self, x, y): + return x, y + + +# Use the first projection module available +try: + import pyproj + ProjTransformer = _PyProjTransformer + HAS_PROJ = True +except ImportError: + try: + from osgeo import osr + ProjTransformer = _OSRTransformer + HAS_PROJ = True + except ImportError: + import logging + logger = logging.getLogger(__name__) + logger.warn('No projection module can be found') + ProjTransformer = _DumbTransformer + HAS_PROJ = False + + +class CustomTileSet(object): + """ Tileset with custom zoom levels and projections """ + + # NOTE 1: + # The MBTiles specs limits its scope to the global-mercator TMS profile. + # Using this class, you will break standard compliance of your MBTiles. + # Consider adding tileset parameters as metadata to enable overlaying of + # tiles. + + # NOTE 2: + # For the moment, the lower left corner of the extent is on tile (0, 0) at + # every zoom level. The TMS spec mentions an 'origin' parameter which + # makes possible to center the extent in the grid when the extent is not + # square. Should we implement this? + + def __init__(self, **kwargs): + """ + The tileset is defined by: + * tilesize: the pixel size of a tile (square tiles assumed). + Default: 256. + * proj: the coordinate system identifier (will be passed to the + underlying projection library). Default: 'EPSG:3857'. + * extent: maximal extent of the tileset expressed in map projection. + Format: (minx, miny, maxx, maxy). Mandatory. + * level_number: the number of zoom level to use. Default: 18. + * max_resolution: the pixel scale at the first zoom level. + Default value will be computed so that extent fits on 1 tile at + the first zoom level. + * resolutions: the pixel scale of each zoom level. If set this will + override level_number and max_resolution. Default to None. + """ + # Tile size + self.tilesize = kwargs.get('tilesize', 256) + + # Coordinate system + self.proj_name = kwargs.get('proj', 'EPSG:3857') + self.transformer = ProjTransformer(self.proj_name) + + # Tileset extent + if 'extent' in kwargs: + # XXX: should we accept and transform lon/lat extent? + if len(kwargs['extent']) != 4: + raise InvalidCoverageError(_("Wrong format of bounding box.")) + if kwargs['extent'][0] >= kwargs['extent'][2] or \ + kwargs['extent'][1] >= kwargs['extent'][3]: + raise InvalidCoverageError( + _("Bounding box format is (xmin, ymin, xmax, ymax)")) + self.extent = tuple(kwargs['extent']) + else: + # XXX: Is it possible to deduce extent automatically from the SRS? + raise TypeError('Mandatory parameter missing: extent') + + # Determine zoom levels + if 'resolutions' in kwargs: + self.resolutions = tuple(kwargs['resolutions']) + else: + if 'max_resolution' in kwargs: + max_res = float(kwargs['max_resolution']) + else: + map_size = max(self.extent[2] - self.extent[0], + self.extent[3] - self.extent[1]) + max_res = float(map_size) / (self.tilesize) + level_number = kwargs.get('level_number', 21) + self.resolutions = tuple([max_res / 2**n + for n in range(level_number)]) + + def __eq__(self, other): + """ Define == operator for TileSet object """ + return self.proj_name == other.proj_name and \ + self.resolutions == other.resolutions and \ + self.tilesize == other.tilesize and \ + self.extent == other.extent + + @property + def NAME(self): + """ A label for the map projection of this tileset """ + return self.proj_name + + def project_pixels(self, ll, zoom): + """ Return the pixel coordinates for the specified lon/lat position """ + coords = self.transformer.from_lonlat(ll[0], ll[1]) + res = self.resolutions[zoom] + x = int((coords[0] - self.extent[0]) / res) + y = int((coords[1] - self.extent[1]) / res) + return (x, y,) + + def unproject_pixels(self, px, zoom): + """ Return the lon/lat coordinates for the specified pixel position """ + res = self.resolutions[zoom] + x = (px[0] + 0.5) * res + self.extent[0] + y = (px[1] + 0.5) * res + self.extent[1] + coords = self.transformer.to_lonlat(x, y) + return coords + + def tile_at(self, zoom, position): + """ Return the tile coordinates for the specified lon/lat position """ + # XXX: Move this method to a BaseTileSet class? + x, y = self.project_pixels(position, zoom) + return (zoom, int(x/self.tilesize), int(y/self.tilesize)) + + def tile_bbox(self, (z, x, y)): + """ Return the lon/lat bbox of the specified tile """ + # XXX: Move this method to a BaseTileSet class? + # NOTE: Return the usual (lower-left, upper-right) instead of the + # (upper-left, lower-right) of GoogleProjection instances. + ll_p = (x * self.tilesize, y * self.tilesize) + ur_p = ((x+1) * self.tilesize - 1, (y+1) * self.tilesize - 1) + ll_g = self.unproject_pixels(ll_p, z) + ur_g = self.unproject_pixels(ur_p, z) + return ll_g + ur_g + + def project(self, (lng, lat)): + """ Convert coordinates from lon/lat to map projection """ + return self.transformer.from_lonlat(lng, lat) + + def unproject(self, (x, y)): + """ Convert coordinates from map projection to lon/lat """ + return self.transformer.to_lonlat(lng, lat) + + def tileslist(self, bbox, levels=None): + """ Return the subset of tiles within the specified lon/lat bbox """ + # XXX: Move this method to a BaseTileSet class? + # XXX: Could this be a generator? + if len(bbox) != 4: + raise InvalidCoverageError(_("Wrong format of bounding box.")) + xmin, ymin, xmax, ymax = bbox + if abs(xmin) > 180 or abs(xmax) > 180 or \ + abs(ymin) > 90 or abs(ymax) > 90: + raise InvalidCoverageError( + _("Some coordinates exceed [-180,+180], [-90, 90].")) + if xmin >= xmax or ymin >= ymax: + raise InvalidCoverageError( + _("Bounding box format is (xmin, ymin, xmax, ymax)")) + + l = [] + if levels is None: + levels = range(len(self.resolutions)) + for z in levels: + if z < 0 or z >= len(self.resolutions): + continue + ll_tile = self.tile_at(z, (xmin, ymin)) # lower left tile + ur_tile = self.tile_at(z, (xmax, ymax)) # upper right tile + + for x in range(ll_tile[1], ur_tile[1]+1): + for y in range(ll_tile[2], ur_tile[2]+1): + l.append((z, x, y)) + return l diff --git a/landez/tests.py b/landez/tests.py index 4306c82..315c551 100644 --- a/landez/tests.py +++ b/landez/tests.py @@ -4,7 +4,7 @@ from tiles import (TilesManager, MBTilesBuilder, ImageExporter, EmptyCoverageError, DownloadError) -from proj import InvalidCoverageError +from proj import InvalidCoverageError, ProjTransformer, HAS_PROJ, CustomTileSet from cache import Disk from sources import MBTilesReader @@ -202,6 +202,85 @@ def test_cache_folder(self): self.assertEqual(mb.cache.folder, '/tmp/landez/servercolortoalphaffffff') +class TestProjBindings(unittest.TestCase): + + def test_no_projection_required(self): + # There should not be any difference between lon/lat and EPSG:4326 + t = ProjTransformer('EPSG:4326') + self.assertEqual(t.to_lonlat(5.386362, 43.293414), (5.386362, 43.293414,)) + self.assertEqual(t.from_lonlat(5.386362, 43.293414), (5.386362, 43.293414,)) + + def test_common_projection(self): + if not HAS_PROJ: + self.skipTest('No projection library available') + p_lonlat = (5.386362, 43.293414,) + t = ProjTransformer('EPSG:3857') # "Google" Mercator + p_ref = (599607.075, 5356739.624,) + p = t.from_lonlat(*p_lonlat) + self.assertAlmostEqual(p[0], p_ref[0], places=3) + self.assertAlmostEqual(p[1], p_ref[1], places=3) + p = t.to_lonlat(*p_ref) + self.assertAlmostEqual(p[0], p_lonlat[0], places=3) + self.assertAlmostEqual(p[1], p_lonlat[1], places=3) + t = ProjTransformer('EPSG:2154') # Lambert 93 + p_ref = (893744.752, 6246732.852,) + p = t.from_lonlat(*p_lonlat) + self.assertAlmostEqual(p[0], p_ref[0], places=3) + self.assertAlmostEqual(p[1], p_ref[1], places=3) + p = t.to_lonlat(*p_ref) + self.assertAlmostEqual(p[0], p_lonlat[0], places=3) + self.assertAlmostEqual(p[1], p_lonlat[1], places=3) + + +class TestCustomTS(unittest.TestCase): + + def test_res_init(self): + ts1 = CustomTileSet(proj='EPSG:4326', + extent=(2.5, 40, 12.5, 50), + max_resolution=10./512, + level_number=3) + ts2 = CustomTileSet(proj='EPSG:4326', + extent=(2.5, 30, 12.5, 60), + resolutions=(10./512, 10./1024, 10./2048)) + self.assertEqual(ts1.resolutions, ts2.resolutions) + + def test_tile_at(self): + ts = CustomTileSet(proj='EPSG:4326', + extent=(2.5, 40, 12.5, 50), + level_number=3) + # FIXME: test breaks because of floating point imprecision + for z in range(len(ts.resolutions)): + res = ts.resolutions[z] + tile = ts.tile_at(z, (2.5 + 255*res, 40 + 255*res)) + self.assertEqual(tile, (0, 0, 0)) + tile = ts.tile_at(z, (2.5 + 256*res, 40 + 256*res)) + self.assertEqual(tile, (0, 1, 1)) + + def test_eq_operator(self): + ts1 = CustomTileSet(proj='EPSG:4326', + extent=(2.5, 40, 12.5, 50), + level_number=3) + ts2 = CustomTileSet(proj='EPSG:4326', + extent=(2.5, 40, 12.5, 50), + resolutions=(10./256, 10./512, 10./1024)) + ts3 = CustomTileSet(proj='EPSG:4326', + extent=(2.5, 30, 12.5, 60), + level_number=3) + self.assertTrue(ts1 == ts2) + self.assertFalse(ts1 == ts3) + + def test_tileslist(self): + ts = CustomTileSet(proj='EPSG:4326', + extent=(2.5, 40, 12.5, 50), + level_number=3) + tsize = 256 * ts.resolutions[2] + bbox= (7.5 - tsize/2, 45 - tsize/2, 7.5 + tsize/2, 45 + tsize/2) + tlist = ts.tileslist(bbox, range(3)) + self.assertEqual(tlist, [(0, 0, 0), (1, 0, 0), (1, 0, 1), (1, 1, 0), + (1, 1, 1), (2, 1, 1), (2, 1, 2), (2, 2, 1), (2, 2, 2)]) + tlist = ts.tileslist(bbox, [2]) + self.assertEqual(tlist, [(2, 1, 1), (2, 1, 2), (2, 2, 1), (2, 2, 2)]) + if __name__ == '__main__': logging.basicConfig(level=logging.DEBUG) unittest.main()