Coverage for gws-app/gws/lib/grid/__init__.py: 87%
144 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-10-23 20:06 +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(
46 gws.lib.grid.Options(
47 crs=gws.lib.crs.get(25832),
48 extent=(280000, 5200000, 920000, 6100000),
49 tileSize=512,
50 )
51 )
52"""
54import math
55from typing import Optional
56from collections.abc import Iterator
58import gws
59import gws.lib.crs
61RESOLUTION_TOLERANCE = 0.01
62"""Relative tolerance for a resolution to count as a level's resolution.
64Requests are often rounded (bboxes with a few decimals), so a resolution slightly above
65a level's must still map to that level rather than to the next finer one.
66"""
68UPSCALE_BIAS = 2
69"""Preference for upscaling a source matrix over downscaling one, as a factor on a log scale.
711 picks the matrix with the smaller scale factor, values above 1 nudge towards upscaling the
72coarser matrix, values below 1 towards downscaling the finer one. Downscaling makes labels and
73thin lines unreadable, upscaling only makes them soft. With a 2x ladder, a matrix is upscaled by
74at most ``2 ** (UPSCALE_BIAS / (UPSCALE_BIAS + 1))``, about 1.59 for the default.
75"""
77DEFAULT_TILE_SIZE = 256
78"""Default tile size in pixels."""
81class Props(gws.Props):
82 origin: str
83 extent: gws.Extent
84 resolutions: list[float]
85 tileSize: int
88class Config(gws.Config):
89 """Tile grid. (changed in 8.5)"""
91 crs: Optional[gws.CrsName]
92 """Grid CRS."""
93 extent: Optional[gws.Extent]
94 """Grid frame extent, with the origin at its north-west corner."""
95 baseResolution: Optional[float]
96 """Resolution at level 0. (added in 8.5)"""
97 tileSize: Optional[int]
98 """Tile size in pixels."""
99 withSnap: bool = True
100 """Snap the extent and the base resolution to the default grid of the CRS. (added in 8.5)"""
103class Options(gws.Data):
104 """Options for creating a grid with ``new``."""
106 crs: gws.Crs
107 """Grid CRS."""
108 extent: Optional[gws.Extent]
109 """Frame extent. The default depends on the CRS."""
110 baseResolution: Optional[float]
111 """Resolution at level 0. By default, one tile spans the frame height."""
112 tileSize: Optional[int]
113 """Tile size in pixels, ``DEFAULT_TILE_SIZE`` by default."""
114 withSnap: Optional[bool]
115 """Snap a custom extent or base resolution to the default grid of the CRS, ``True`` by default."""
118def for_crs(crs: gws.Crs) -> gws.MapGrid:
119 """Create the default grid for a CRS.
121 Args:
122 crs: Grid CRS.
124 Returns:
125 A grid with the default frame, base resolution and tile size.
126 """
128 return new(Options(crs=crs))
131def new(opts: Options) -> gws.MapGrid:
132 """Create a grid.
134 Missing options are filled with the defaults for the CRS. If a custom extent or base resolution
135 is given and snapping is not disabled, the grid is snapped to the default grid of the CRS.
137 Args:
138 opts: Grid options.
140 Returns:
141 A new grid.
142 """
144 mg = gws.MapGrid()
145 mg.crs = opts.crs
146 mg.extent = opts.extent or (gws.lib.crs.WGS84.extent if mg.crs.isGeographic else gws.lib.crs.WEBMERCATOR_SQUARE)
147 mg.tileSize = opts.tileSize or DEFAULT_TILE_SIZE
148 mg.baseResolution = opts.baseResolution or (mg.extent[3] - mg.extent[1]) / mg.tileSize
149 with_snap = True if opts.withSnap is None else opts.withSnap
150 if with_snap and (opts.extent or opts.baseResolution):
151 _snap(mg, bool(opts.baseResolution))
152 return mg
155def _snap(mg: gws.MapGrid, has_base_resolution: bool):
156 """Snap the base resolution and the extent of a grid to the default grid of its CRS, in place."""
158 ref = for_crs(mg.crs)
160 if has_base_resolution:
161 k = round(math.log2(ref.baseResolution / mg.baseResolution))
162 else:
163 w = mg.extent[2] - mg.extent[0]
164 h = mg.extent[3] - mg.extent[1]
165 k = math.floor(math.log2(ref.baseResolution * mg.tileSize / max(w, h)) + 1e-9)
167 mg.baseResolution = resolution_for_level(ref, max(k, 0))
169 span = mg.baseResolution * mg.tileSize
170 eps = span * 1e-6
171 ox = ref.extent[0]
172 oy = ref.extent[3]
174 mg.extent = (
175 ox + math.floor((mg.extent[0] - ox + eps) / span) * span,
176 oy - math.ceil((oy - mg.extent[1] - eps) / span) * span,
177 ox + math.ceil((mg.extent[2] - ox - eps) / span) * span,
178 oy - math.floor((oy - mg.extent[3] + eps) / span) * span,
179 )
182def resolution_for_level(mg: gws.MapGrid, z: int) -> float:
183 """Return the resolution of a level.
185 Args:
186 mg: A grid.
187 z: Level.
189 Returns:
190 ``baseResolution / 2**z``.
191 """
193 return mg.baseResolution / (1 << z)
196def level_for_resolution(mg: gws.MapGrid, resolution: float) -> int:
197 """Return the coarsest level whose resolution does not exceed the given one.
199 A level whose resolution is within ``RESOLUTION_TOLERANCE`` of the given one also matches.
201 Args:
202 mg: A grid.
203 resolution: Resolution in CRS units per pixel.
205 Returns:
206 Level number.
208 Raises:
209 ``ValueError``: If the resolution is not positive or no level up to 99 matches.
210 """
212 if resolution <= 0:
213 raise ValueError(f'invalid resolution {resolution!r}')
214 for z in range(100):
215 r = resolution_for_level(mg, z)
216 if r <= resolution or math.isclose(r, resolution, rel_tol=RESOLUTION_TOLERANCE):
217 return z
218 raise ValueError(f'invalid resolution {resolution!r}')
221def props_for_resolutions(mg: gws.MapGrid, resolutions: list[float]) -> Props:
222 """Create client props for a grid.
224 Args:
225 mg: A grid.
226 resolutions: Resolutions the client uses. The finest one determines the finest level.
228 Returns:
229 Grid props with the resolutions of levels ``0`` up to the finest level needed.
230 """
232 zmax = level_for_resolution(mg, min(resolutions))
233 return Props(
234 origin=gws.Origin.nw,
235 extent=mg.extent,
236 resolutions=[resolution_for_level(mg, z) for z in range(zmax + 1)],
237 tileSize=mg.tileSize,
238 )
241def tile_count_for_level(mg: gws.MapGrid, z: int) -> tuple[int, int]:
242 """Return the number of tiles a level has across the frame.
244 Args:
245 mg: A grid.
246 z: Level.
248 Returns:
249 Number of columns and rows, at least 1 each.
250 """
252 span = resolution_for_level(mg, z) * mg.tileSize
253 return (
254 max(1, math.ceil((mg.extent[2] - mg.extent[0]) / span - 1e-6)),
255 max(1, math.ceil((mg.extent[3] - mg.extent[1]) / span - 1e-6)),
256 )
259def range_for_extent(mg: gws.MapGrid, extent: gws.Extent, z: int) -> gws.MapTileRange | None:
260 """Return the range of tiles covering an extent at a level.
262 The range is clipped to the grid frame.
264 Args:
265 mg: A grid.
266 extent: Extent in the grid CRS.
267 z: Level.
269 Returns:
270 A tile range, or ``None`` if the extent does not intersect the frame.
271 """
273 span = resolution_for_level(mg, z) * mg.tileSize
274 nx, ny = tile_count_for_level(mg, z)
275 eps = span * 1e-6
277 x0 = math.floor((extent[0] - mg.extent[0] + eps) / span)
278 x1 = math.floor((extent[2] - mg.extent[0] - eps) / span)
279 y0 = math.floor((mg.extent[3] - extent[3] + eps) / span)
280 y1 = math.floor((mg.extent[3] - extent[1] - eps) / span)
282 if x1 < x0 or y1 < y0 or x1 < 0 or y1 < 0 or x0 >= nx or y0 >= ny:
283 return None
284 return max(x0, 0), max(y0, 0), min(x1, nx - 1), min(y1, ny - 1), z
287def intersect_ranges(a: gws.MapTileRange, b: gws.MapTileRange) -> gws.MapTileRange | None:
288 """Return the intersection of two tile ranges of the same level.
290 Args:
291 a: First tile range.
292 b: Second tile range.
294 Returns:
295 The intersection, or ``None`` if the ranges do not intersect.
297 Raises:
298 ``ValueError``: If the ranges are at different levels.
299 """
301 if a[4] != b[4]:
302 raise ValueError(f'cannot intersect ranges of different levels: {a!r}, {b!r}')
303 x0, y0, x1, y1 = max(a[0], b[0]), max(a[1], b[1]), min(a[2], b[2]), min(a[3], b[3])
304 if x1 < x0 or y1 < y0:
305 return None
306 return x0, y0, x1, y1, a[4]
309def extent_for_range(mg: gws.MapGrid, tr: gws.MapTileRange) -> gws.Extent:
310 """Return the extent of a tile range.
312 Args:
313 mg: A grid.
314 tr: Tile range.
316 Returns:
317 Extent in the grid CRS.
318 """
320 x0, y0, x1, y1, z = tr
321 span = resolution_for_level(mg, z) * mg.tileSize
322 return (
323 mg.extent[0] + x0 * span,
324 mg.extent[3] - (y1 + 1) * span,
325 mg.extent[0] + (x1 + 1) * span,
326 mg.extent[3] - y0 * span,
327 )
330def extent_for_tile(mg: gws.MapGrid, tile: gws.MapTile) -> gws.Extent:
331 """Return the extent of a tile.
333 Args:
334 mg: A grid.
335 tile: Tile.
337 Returns:
338 Extent in the grid CRS.
339 """
341 x, y, z = tile
342 return extent_for_range(mg, (x, y, x, y, z))
345def matrix_set_for_grid(mg: gws.MapGrid, max_level: int, identifier: str = '') -> gws.TileMatrixSet:
346 """Express the levels ``0..max_level`` of a grid as a tile matrix set.
348 Matrix identifiers are the level numbers as strings. ``scale`` is not set.
350 Args:
351 mg: A grid.
352 max_level: Finest level to include.
353 identifier: Identifier of the matrix set.
355 Returns:
356 A tile matrix set.
357 """
359 tms = gws.TileMatrixSet(identifier=identifier, crs=mg.crs, matrices=[])
360 for z in range(max_level + 1):
361 nx, ny = tile_count_for_level(mg, z)
362 tms.matrices.append(
363 gws.TileMatrix(
364 identifier=str(z),
365 resolution=resolution_for_level(mg, z),
366 x=mg.extent[0],
367 y=mg.extent[3],
368 width=nx,
369 height=ny,
370 tileWidth=mg.tileSize,
371 tileHeight=mg.tileSize,
372 extent=mg.extent,
373 )
374 )
375 return tms
378def best_resolution(resolution: float, resolutions: list[float]) -> float:
379 """Return the resolution to render a resolution from.
381 Of the two resolutions around the given one, the coarser one needs upscaling and the
382 finer one downscaling. The scale factors are compared on a log scale, weighted by ``UPSCALE_BIAS``.
384 Args:
385 resolution: Wanted resolution.
386 resolutions: Available resolutions, in any order.
388 Returns:
389 The matching resolution, or the finest one if all are coarser.
390 """
392 by_res = sorted(resolutions, reverse=True)
393 prev = None
395 for r in by_res:
396 if math.isclose(r, resolution, rel_tol=RESOLUTION_TOLERANCE):
397 return r
398 if r > resolution:
399 prev = r
400 continue
401 if prev is None:
402 return r
403 up_scale_factor = math.log(prev / resolution)
404 down_scale_factor = math.log(resolution / r)
405 if up_scale_factor <= UPSCALE_BIAS * down_scale_factor:
406 return prev
407 return r
409 return by_res[-1]
412def matrix_for_resolution(tms: gws.TileMatrixSet, resolution: float) -> gws.TileMatrix:
413 """Return the matrix to render a resolution from, see ``best_resolution``.
415 Args:
416 tms: Tile matrix set.
417 resolution: Wanted resolution in CRS units per pixel.
419 Returns:
420 The matching matrix, or the finest one if all are coarser.
421 """
423 r = best_resolution(resolution, [tm.resolution for tm in tms.matrices])
424 for tm in tms.matrices:
425 if tm.resolution == r:
426 return tm
427 raise ValueError(f'no matrix found for resolution {r}')
430def matrix_range_for_extent(tm: gws.TileMatrix, extent: gws.Extent) -> gws.MapTileRange | None:
431 """Return the range of matrix tiles covering an extent.
433 The range is clipped to the matrix.
435 Args:
436 tm: Tile matrix.
437 extent: Extent in the matrix CRS.
439 Returns:
440 A tile range with ``z`` set to 0, or ``None`` if the extent does not intersect the matrix.
441 """
443 tile_w = tm.resolution * tm.tileWidth
444 tile_h = tm.resolution * tm.tileHeight
446 # nudge the edges inward by a fraction of a tile, so that an edge lying exactly
447 # on a tile boundary does not pull in a neighbouring tile through float noise
448 eps = 1e-6
450 x0 = math.floor((extent[0] - tm.x + tile_w * eps) / tile_w)
451 x1 = math.floor((extent[2] - tm.x - tile_w * eps) / tile_w)
452 y0 = math.floor((tm.y - extent[3] + tile_h * eps) / tile_h)
453 y1 = math.floor((tm.y - extent[1] - tile_h * eps) / tile_h)
455 if x1 < x0 or y1 < y0 or x1 < 0 or y1 < 0 or x0 >= tm.width or y0 >= tm.height:
456 return None
457 return max(x0, 0), max(y0, 0), min(x1, int(tm.width) - 1), min(y1, int(tm.height) - 1), 0
460def matrix_extent_for_range(tm: gws.TileMatrix, tr: gws.MapTileRange) -> gws.Extent:
461 """Return the extent of a range of matrix tiles.
463 Args:
464 tm: Tile matrix.
465 tr: Tile range; ``z`` is ignored.
467 Returns:
468 Extent in the matrix CRS.
469 """
471 tile_w = tm.resolution * tm.tileWidth
472 tile_h = tm.resolution * tm.tileHeight
473 x0, y0, x1, y1, _ = tr
474 return (
475 tm.x + x0 * tile_w,
476 tm.y - (y1 + 1) * tile_h,
477 tm.x + (x1 + 1) * tile_w,
478 tm.y - y0 * tile_h,
479 )
482def enum_tiles(tr: gws.MapTileRange) -> Iterator[gws.MapTile]:
483 """Enumerate the tiles of a range, row by row.
485 Args:
486 tr: Tile range.
488 Yields:
489 Tiles ``(x, y, z)``.
490 """
492 x0, y0, x1, y1, z = tr
493 for y in range(y0, y1 + 1):
494 for x in range(x0, x1 + 1):
495 yield x, y, z
498def in_range(mt: gws.MapTile, tr: gws.MapTileRange) -> bool:
499 """Check if a tile lies within a range.
501 Args:
502 mt: Tile.
503 tr: Tile range.
505 Returns:
506 ``True`` if the tile is at the level of the range and inside it.
507 """
509 x, y, z = mt
510 return z == tr[4] and tr[0] <= x <= tr[2] and tr[1] <= y <= tr[3]