Add support for local projections

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

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()