110 lines
3.8 KiB
Python
110 lines
3.8 KiB
Python
import numpy as np
|
|
import rasterio
|
|
from rasterio.transform import from_origin
|
|
from rasterio.warp import transform_bounds
|
|
|
|
from pipeline.transform import noise_overlay_tiles
|
|
from pipeline.transform.noise_overlay_tiles import RasterInfo, _read_noise_tile
|
|
|
|
|
|
def _write_corridor_raster(path, nodata=-96.0):
|
|
"""A small EPSG:27700 raster: a column of 70 dB cells adjacent to genuine
|
|
0.0 (quiet) cells. Bilinear blending of the 0 cells would fabricate a halo
|
|
of intermediate dB values between 0 and 70."""
|
|
# 8x8 grid: leftmost two columns are 70 dB, the rest are genuine quiet 0.0.
|
|
data = np.zeros((8, 8), dtype=np.float32)
|
|
data[:, 0:2] = 70.0
|
|
# Place one true nodata cell to make sure it is also masked out.
|
|
data[0, 7] = nodata
|
|
|
|
# 10m cells anchored somewhere inside England's BNG extent.
|
|
left = 300_000.0
|
|
top = 300_080.0
|
|
transform = from_origin(left, top, 10.0, 10.0)
|
|
with rasterio.open(
|
|
path,
|
|
"w",
|
|
driver="GTiff",
|
|
height=data.shape[0],
|
|
width=data.shape[1],
|
|
count=1,
|
|
dtype=data.dtype,
|
|
crs="EPSG:27700",
|
|
transform=transform,
|
|
nodata=nodata,
|
|
) as dataset:
|
|
dataset.write(data, 1)
|
|
return path
|
|
|
|
|
|
def test_read_noise_tile_does_not_fabricate_halo(tmp_path):
|
|
raster_path = _write_corridor_raster(tmp_path / "corridor.tif")
|
|
|
|
with rasterio.open(raster_path) as dataset:
|
|
bounds_27700 = dataset.bounds
|
|
bounds_mercator = transform_bounds(
|
|
dataset.crs,
|
|
noise_overlay_tiles.WEB_MERCATOR_CRS,
|
|
*bounds_27700,
|
|
densify_pts=21,
|
|
)
|
|
|
|
info = RasterInfo(path=raster_path, bounds_mercator=bounds_mercator)
|
|
|
|
# Render at high resolution so any bilinear halo would surface as
|
|
# intermediate dB values along the corridor/quiet seam.
|
|
tile_size = 64
|
|
tile = _read_noise_tile([info], bounds_mercator, tile_size)
|
|
|
|
finite = tile[np.isfinite(tile)]
|
|
# Every finite cell must be the genuine corridor value (~70). There must be
|
|
# NO fabricated halo strictly between 0 and 70.
|
|
halo = finite[(finite > 0.0) & (finite < 70.0 - 1e-3)]
|
|
assert halo.size == 0, f"fabricated halo values present: {np.unique(halo)}"
|
|
# Sanity: the corridor itself must still be rendered.
|
|
assert finite.size > 0
|
|
assert np.all(finite >= 70.0 - 1e-3)
|
|
|
|
|
|
def test_read_noise_tile_preserves_peak_under_downsample(tmp_path):
|
|
# 8x8 EPSG:27700 raster: a single loud 75 dB cell in a 50 dB field.
|
|
# Downsampling into a smaller tile with bilinear would dilute the peak
|
|
# (arithmetic dB averaging); Resampling.max must keep the worst-case dB.
|
|
data = np.full((8, 8), 50.0, dtype=np.float32)
|
|
data[4, 4] = 75.0
|
|
transform = from_origin(300_000.0, 300_080.0, 10.0, 10.0)
|
|
raster_path = tmp_path / "peak.tif"
|
|
with rasterio.open(
|
|
raster_path,
|
|
"w",
|
|
driver="GTiff",
|
|
height=data.shape[0],
|
|
width=data.shape[1],
|
|
count=1,
|
|
dtype=data.dtype,
|
|
crs="EPSG:27700",
|
|
transform=transform,
|
|
nodata=-96.0,
|
|
) as dataset:
|
|
dataset.write(data, 1)
|
|
|
|
with rasterio.open(raster_path) as dataset:
|
|
bounds_mercator = transform_bounds(
|
|
dataset.crs,
|
|
noise_overlay_tiles.WEB_MERCATOR_CRS,
|
|
*dataset.bounds,
|
|
densify_pts=21,
|
|
)
|
|
|
|
info = RasterInfo(path=raster_path, bounds_mercator=bounds_mercator)
|
|
|
|
# Render the 8x8 source into a 4x4 tile: this downsamples, so bilinear
|
|
# would average the 75 dB peak away.
|
|
tile = _read_noise_tile([info], bounds_mercator, 4)
|
|
finite = tile[np.isfinite(tile)]
|
|
|
|
assert finite.size > 0
|
|
# The loud peak must survive the downsample (max, not arithmetic mean).
|
|
assert finite.max() >= 75.0 - 1e-3, f"peak diluted to {finite.max()}"
|
|
# Max resampling must never invent a value louder than the source.
|
|
assert finite.max() <= 75.0 + 1e-3
|