Coverage for gws-app/gws/lib/grid/__init__.py: 86%
124 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-05 13:35 +0200
1"""Map grids: tile pyramid math.
3A grid is a ``gws.MapGrid``: a CRS, a frame extent, a base resolution and a
4tile size. The origin is the north-west corner of the frame. Level ``z`` has
5the resolution ``baseResolution / 2**z``, so levels nest exactly (quad tree)
6and the ladder has no bottom: any resolution can be served by downsampling
7from the coarsest level whose resolution does not exceed it.
9Defaults: projected CRS use the web mercator square as frame and one tile at
10level 0; geographic CRS use ``(-180, -90, 180, 90)`` and two tiles (2x1) at
11level 0. Without an explicit base resolution, one tile at level 0 spans the
12frame height, so that the default resolutions are
13``156543.03392804097 / 2**z`` metres and ``0.703125 / 2**z`` degrees for a
14256px tile. Custom frames, base resolutions and tile sizes (tile source
15grids) are possible; a non-square custom frame needs an explicit base
16resolution. The frame is an indexing frame and does not have to be inside
17the area of use of the CRS.
19A custom grid can be snapped to the default grid of its CRS (``withSnap``):
20the base resolution becomes a ladder value, the nearest one if given,
21otherwise the coarsest whose tile spans the extent, and the extent is
22expanded outward to tile boundaries at that level. A snapped grid is a
23window on the default grid: its tile ``(x, y, z)`` is the default grid's
24``(x + ox * 2**z, y + oy * 2**z, z + k)`` for the integer offset ``(ox, oy)``
25of the window at level ``k``.
27Tiles are ``(x, y, z)``, ranges ``(min_x, min_y, max_x, max_y, z)``.
29Tile matrix sets (``gws.TileMatrixSet``, as in WMTS) describe foreign
30pyramids: each matrix has its own origin, size and resolution, and the
31resolutions need not form a ladder. ``matrix_set_for_grid`` expresses a grid
32as a matrix set, the ``matrix_*`` functions are the matrix counterparts of
33the grid functions.
35Example::
37 mg = gws.lib.grid.for_crs(gws.lib.crs.get(3857))
38 z = gws.lib.grid.level_for_resolution(mg, 10.0)
39 tr = gws.lib.grid.range_for_extent(mg, (1000000, 6000000, 1010000, 6010000), z)
40 for mt in gws.lib.grid.enum_tiles(tr):
41 extent = gws.lib.grid.extent_for_tile(mg, mt)
43A custom grid from a ``Config``-like set of options::
45 mg = gws.lib.grid.new(gws.lib.grid.Options(
46 crs=gws.lib.crs.get(25832),
47 extent=(280000, 5200000, 920000, 6100000),
48 tileSize=512,
49 ))
50"""
52import math
53from typing import Iterator, Optional
55import gws
56import gws.lib.crs
58RESOLUTION_TOLERANCE = 0.01
59"""Relative tolerance for a resolution to count as a level's resolution.
61Requests are often rounded (bboxes with a few decimals), so a resolution slightly above
62a level's must still map to that level rather than to the next finer one.
63"""
65DEFAULT_TILE_SIZE = 256
66"""Default tile size in pixels."""
69class Props(gws.Props):
70 origin: str
71 extent: gws.Extent
72 resolutions: list[float]
73 tileSize: int
78class Config(gws.Config):
79 """Tile grid. (changed in 8.5)"""
81 crs: Optional[gws.CrsName]
82 """Grid CRS."""
83 extent: Optional[gws.Extent]
84 """Grid frame extent, with the origin at its north-west corner."""
85 baseResolution: Optional[float]
86 """Resolution at level 0. (added in 8.5)"""
87 tileSize: Optional[int]
88 """Tile size in pixels."""
89 withSnap: bool = True
90 """Snap the extent and the base resolution to the default grid of the CRS. (added in 8.5)"""
93class Options(gws.Data):
94 """Options for creating a grid with ``new``."""
96 crs: gws.Crs
97 """Grid CRS."""
98 extent: Optional[gws.Extent]
99 """Frame extent. The default depends on the CRS."""
100 baseResolution: Optional[float]
101 """Resolution at level 0. By default, one tile spans the frame height."""
102 tileSize: Optional[int]
103 """Tile size in pixels, ``DEFAULT_TILE_SIZE`` by default."""
104 withSnap: Optional[bool]
105 """Snap a custom extent or base resolution to the default grid of the CRS, ``True`` by default."""
108def for_crs(crs: gws.Crs) -> gws.MapGrid:
109 """Create the default grid for a CRS.
111 Args:
112 crs: Grid CRS.
114 Returns:
115 A grid with the default frame, base resolution and tile size.
116 """
118 return new(Options(crs=crs))
121def new(opts: Options) -> gws.MapGrid:
122 """Create a grid.
124 Missing options are filled with the defaults for the CRS. If a custom extent or base resolution
125 is given and snapping is not disabled, the grid is snapped to the default grid of the CRS.
127 Args:
128 opts: Grid options.
130 Returns:
131 A new grid.
132 """
134 mg = gws.MapGrid()
135 mg.crs = opts.crs
136 mg.extent = opts.extent or (gws.lib.crs.WGS84.extent if mg.crs.isGeographic else gws.lib.crs.WEBMERCATOR_SQUARE)
137 mg.tileSize = opts.tileSize or DEFAULT_TILE_SIZE
138 mg.baseResolution = opts.baseResolution or (mg.extent[3] - mg.extent[1]) / mg.tileSize
139 with_snap = True if opts.withSnap is None else opts.withSnap
140 if with_snap and (opts.extent or opts.baseResolution):
141 _snap(mg, bool(opts.baseResolution))
142 return mg
145def _snap(mg: gws.MapGrid, has_base_resolution: bool):
146 """Snap the base resolution and the extent of a grid to the default grid of its CRS, in place."""
148 ref = for_crs(mg.crs)
150 if has_base_resolution:
151 k = round(math.log2(ref.baseResolution / mg.baseResolution))
152 else:
153 w = mg.extent[2] - mg.extent[0]
154 h = mg.extent[3] - mg.extent[1]
155 k = math.floor(math.log2(ref.baseResolution * mg.tileSize / max(w, h)) + 1e-9)
157 mg.baseResolution = resolution_for_level(ref, max(k, 0))
159 span = mg.baseResolution * mg.tileSize
160 eps = span * 1e-6
161 ox = ref.extent[0]
162 oy = ref.extent[3]
164 mg.extent = (
165 ox + math.floor((mg.extent[0] - ox + eps) / span) * span,
166 oy - math.ceil((oy - mg.extent[1] - eps) / span) * span,
167 ox + math.ceil((mg.extent[2] - ox - eps) / span) * span,
168 oy - math.floor((oy - mg.extent[3] + eps) / span) * span,
169 )
172def resolution_for_level(mg: gws.MapGrid, z: int) -> float:
173 """Return the resolution of a level.
175 Args:
176 mg: A grid.
177 z: Level.
179 Returns:
180 ``baseResolution / 2**z``.
181 """
183 return mg.baseResolution / (1 << z)
186def level_for_resolution(mg: gws.MapGrid, resolution: float) -> int:
187 """Return the coarsest level whose resolution does not exceed the given one.
189 A level whose resolution is within ``RESOLUTION_TOLERANCE`` of the given one also matches.
191 Args:
192 mg: A grid.
193 resolution: Resolution in CRS units per pixel.
195 Returns:
196 Level number.
198 Raises:
199 ``ValueError``: If the resolution is not positive or no level up to 99 matches.
200 """
202 if resolution <= 0:
203 raise ValueError(f'invalid resolution {resolution!r}')
204 for z in range(100):
205 r = resolution_for_level(mg, z)
206 if r <= resolution or math.isclose(r, resolution, rel_tol=RESOLUTION_TOLERANCE):
207 return z
208 raise ValueError(f'invalid resolution {resolution!r}')
211def props_for_resolutions(mg: gws.MapGrid, resolutions: list[float]) -> Props:
212 """Create client props for a grid.
214 Args:
215 mg: A grid.
216 resolutions: Resolutions the client uses. The finest one determines the finest level.
218 Returns:
219 Grid props with the resolutions of levels ``0`` up to the finest level needed.
220 """
222 zmax = level_for_resolution(mg, min(resolutions))
223 return Props(
224 origin=gws.Origin.nw,
225 extent=mg.extent,
226 resolutions=[resolution_for_level(mg, z) for z in range(zmax + 1)],
227 tileSize=mg.tileSize,
228 )
231def tile_count_for_level(mg: gws.MapGrid, z: int) -> tuple[int, int]:
232 """Return the number of tiles a level has across the frame.
234 Args:
235 mg: A grid.
236 z: Level.
238 Returns:
239 Number of columns and rows, at least 1 each.
240 """
242 span = resolution_for_level(mg, z) * mg.tileSize
243 return (
244 max(1, math.ceil((mg.extent[2] - mg.extent[0]) / span - 1e-6)),
245 max(1, math.ceil((mg.extent[3] - mg.extent[1]) / span - 1e-6)),
246 )
249def range_for_extent(mg: gws.MapGrid, extent: gws.Extent, z: int) -> gws.MapTileRange | None:
250 """Return the range of tiles covering an extent at a level.
252 The range is clipped to the grid frame.
254 Args:
255 mg: A grid.
256 extent: Extent in the grid CRS.
257 z: Level.
259 Returns:
260 A tile range, or ``None`` if the extent does not intersect the frame.
261 """
263 span = resolution_for_level(mg, z) * mg.tileSize
264 nx, ny = tile_count_for_level(mg, z)
265 eps = span * 1e-6
267 x0 = math.floor((extent[0] - mg.extent[0] + eps) / span)
268 x1 = math.floor((extent[2] - mg.extent[0] - eps) / span)
269 y0 = math.floor((mg.extent[3] - extent[3] + eps) / span)
270 y1 = math.floor((mg.extent[3] - extent[1] - eps) / span)
272 if x1 < x0 or y1 < y0 or x1 < 0 or y1 < 0 or x0 >= nx or y0 >= ny:
273 return None
274 return max(x0, 0), max(y0, 0), min(x1, nx - 1), min(y1, ny - 1), z
277def intersect_ranges(a: gws.MapTileRange, b: gws.MapTileRange) -> gws.MapTileRange | None:
278 """Return the intersection of two tile ranges of the same level.
280 Args:
281 a: First tile range.
282 b: Second tile range.
284 Returns:
285 The intersection, or ``None`` if the ranges do not intersect.
287 Raises:
288 ``ValueError``: If the ranges are at different levels.
289 """
291 if a[4] != b[4]:
292 raise ValueError(f'cannot intersect ranges of different levels: {a!r}, {b!r}')
293 x0, y0, x1, y1 = max(a[0], b[0]), max(a[1], b[1]), min(a[2], b[2]), min(a[3], b[3])
294 if x1 < x0 or y1 < y0:
295 return None
296 return x0, y0, x1, y1, a[4]
299def extent_for_range(mg: gws.MapGrid, tr: gws.MapTileRange) -> gws.Extent:
300 """Return the extent of a tile range.
302 Args:
303 mg: A grid.
304 tr: Tile range.
306 Returns:
307 Extent in the grid CRS.
308 """
310 x0, y0, x1, y1, z = tr
311 span = resolution_for_level(mg, z) * mg.tileSize
312 return (
313 mg.extent[0] + x0 * span,
314 mg.extent[3] - (y1 + 1) * span,
315 mg.extent[0] + (x1 + 1) * span,
316 mg.extent[3] - y0 * span,
317 )
320def extent_for_tile(mg: gws.MapGrid, tile: gws.MapTile) -> gws.Extent:
321 """Return the extent of a tile.
323 Args:
324 mg: A grid.
325 tile: Tile.
327 Returns:
328 Extent in the grid CRS.
329 """
331 x, y, z = tile
332 return extent_for_range(mg, (x, y, x, y, z))
335def matrix_set_for_grid(mg: gws.MapGrid, max_level: int, identifier: str = '') -> gws.TileMatrixSet:
336 """Express the levels ``0..max_level`` of a grid as a tile matrix set.
338 Matrix identifiers are the level numbers as strings. ``scale`` is not set.
340 Args:
341 mg: A grid.
342 max_level: Finest level to include.
343 identifier: Identifier of the matrix set.
345 Returns:
346 A tile matrix set.
347 """
349 tms = gws.TileMatrixSet(identifier=identifier, crs=mg.crs, matrices=[])
350 for z in range(max_level + 1):
351 nx, ny = tile_count_for_level(mg, z)
352 tms.matrices.append(
353 gws.TileMatrix(
354 identifier=str(z),
355 resolution=resolution_for_level(mg, z),
356 x=mg.extent[0],
357 y=mg.extent[3],
358 width=nx,
359 height=ny,
360 tileWidth=mg.tileSize,
361 tileHeight=mg.tileSize,
362 extent=mg.extent,
363 )
364 )
365 return tms
368def matrix_for_resolution(tms: gws.TileMatrixSet, resolution: float) -> gws.TileMatrix:
369 """Return the coarsest matrix that does not need upscaling.
371 A matrix whose resolution is within ``RESOLUTION_TOLERANCE`` of the given one also matches.
373 Args:
374 tms: Tile matrix set.
375 resolution: Wanted resolution in CRS units per pixel.
377 Returns:
378 The matching matrix, or the finest one if all are coarser.
379 """
381 # Coarsest matrix with tm.resolution <= resolution, i.e. never upscale (downscale up to 2x).
382 # Cross-CRS the wanted resolution rarely hits the source ladder, e.g. 3857 -> 25832
383 # at 51N needs 1.6x the target resolution, landing between two levels.
384 # Alternatives: nearest by ratio (upscale up to sqrt(2), coarser cartography, bigger labels)
385 # or a threshold as in MapProxy (allow upscale below a factor, default 1.15).
386 coarsest_first = sorted(tms.matrices, key=lambda m: -m.resolution)
387 for tm in coarsest_first:
388 if tm.resolution <= resolution or math.isclose(tm.resolution, resolution, rel_tol=RESOLUTION_TOLERANCE):
389 return tm
390 return coarsest_first[-1]
393def matrix_range_for_extent(tm: gws.TileMatrix, extent: gws.Extent) -> gws.MapTileRange | None:
394 """Return the range of matrix tiles covering an extent.
396 The range is clipped to the matrix.
398 Args:
399 tm: Tile matrix.
400 extent: Extent in the matrix CRS.
402 Returns:
403 A tile range with ``z`` set to 0, or ``None`` if the extent does not intersect the matrix.
404 """
406 tile_w = tm.resolution * tm.tileWidth
407 tile_h = tm.resolution * tm.tileHeight
409 # nudge the edges inward by a fraction of a tile, so that an edge lying exactly
410 # on a tile boundary does not pull in a neighbouring tile through float noise
411 eps = 1e-6
413 x0 = math.floor((extent[0] - tm.x + tile_w * eps) / tile_w)
414 x1 = math.floor((extent[2] - tm.x - tile_w * eps) / tile_w)
415 y0 = math.floor((tm.y - extent[3] + tile_h * eps) / tile_h)
416 y1 = math.floor((tm.y - extent[1] - tile_h * eps) / tile_h)
418 if x1 < x0 or y1 < y0 or x1 < 0 or y1 < 0 or x0 >= tm.width or y0 >= tm.height:
419 return None
420 return max(x0, 0), max(y0, 0), min(x1, int(tm.width) - 1), min(y1, int(tm.height) - 1), 0
423def matrix_extent_for_range(tm: gws.TileMatrix, tr: gws.MapTileRange) -> gws.Extent:
424 """Return the extent of a range of matrix tiles.
426 Args:
427 tm: Tile matrix.
428 tr: Tile range; ``z`` is ignored.
430 Returns:
431 Extent in the matrix CRS.
432 """
434 tile_w = tm.resolution * tm.tileWidth
435 tile_h = tm.resolution * tm.tileHeight
436 x0, y0, x1, y1, _ = tr
437 return (
438 tm.x + x0 * tile_w,
439 tm.y - (y1 + 1) * tile_h,
440 tm.x + (x1 + 1) * tile_w,
441 tm.y - y0 * tile_h,
442 )
445def enum_tiles(tr: gws.MapTileRange) -> Iterator[gws.MapTile]:
446 """Enumerate the tiles of a range, row by row.
448 Args:
449 tr: Tile range.
451 Yields:
452 Tiles ``(x, y, z)``.
453 """
455 x0, y0, x1, y1, z = tr
456 for y in range(y0, y1 + 1):
457 for x in range(x0, x1 + 1):
458 yield x, y, z
461def in_range(mt: gws.MapTile, tr: gws.MapTileRange) -> bool:
462 """Check if a tile lies within a range.
464 Args:
465 mt: Tile.
466 tr: Tile range.
468 Returns:
469 ``True`` if the tile is at the level of the range and inside it.
470 """
472 x, y, z = mt
473 return z == tr[4] and tr[0] <= x <= tr[2] and tr[1] <= y <= tr[3]