From 9141e2eaac3c121b75a5b925c516f9ebe0f3aa05 Mon Sep 17 00:00:00 2001 From: Balthasar Indermuehle Date: Fri, 24 Jul 2026 13:58:46 +1000 Subject: [PATCH] Add Copernicus DEM (GLO-90/GLO-30) terrain backend Add two SrtmConf servers, copernicus_glo90 and copernicus_glo30, that read the Copernicus DEM Cloud-Optimised GeoTIFFs from the AWS Open Data buckets. These give a reliably working auto-download path (the SRTM servers are effectively retired), global pole-to-pole coverage, and void-free terrain over water. The GeoTIFF reader reconciles the two ways Copernicus tiles differ from SRTM .hgt: pixel-centre (area) registration vs grid-node, and the latitude- dependent longitude spacing above |50 deg| (tiles are not square there). The public API (srtm_height_data, srtm_height_profile, PathProp) is unchanged. rasterio is used behind an import guard, matching the pattern in gis.py. Add two SrtmConf failure-mode options: on_missing ('zeros'|'raise') and void_fill ('zero'|'nan'|'interp'), both defaulting to the historic behaviour. Fix the SRTM void mask: the canonical void sentinel -32768 was not masked (only -32767/+32767 were), so genuine voids leaked through and were linearly blended with valid neighbours into deep spurious pits. Include tests, docs with the required Copernicus attribution, and a changelog entry. --- CHANGES.rst | 31 ++ docs/pathprof/working_with_srtm.rst | 77 +++++ pycraf/pathprof/srtm.py | 460 +++++++++++++++++++++++++--- pycraf/pathprof/tests/test_srtm.py | 381 +++++++++++++++++++++++ 4 files changed, 900 insertions(+), 49 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index f693e79a..2216d465 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -1,3 +1,34 @@ +unreleased +======================= + +New Features +------------ +pycraf.pathprof +^^^^^^^^^^^^^^ +- Add support for the Copernicus DEM (GLO-90 and GLO-30) as a terrain + source, via two new `pycraf.pathprof.SrtmConf` servers, + ``'copernicus_glo90'`` and ``'copernicus_glo30'``. These are hosted as + Cloud-Optimised GeoTIFFs on the AWS Open Data buckets (no authentication), + are global (pole-to-pole) and void-free over water. The reader handles the + pixel-centre registration and latitude-dependent longitude spacing of the + Copernicus tiles, so the public API is unchanged. Reading the GeoTIFF tiles + needs the optional ``rasterio`` package. +- Add a `pycraf.pathprof.SrtmConf` option ``on_missing`` (``'zeros'`` or + ``'raise'``) to control whether a tile that is missing on disk becomes zero + terrain (with a warning, the historic behaviour) or raises a + ``TileNotAvailableOnDiskError``. +- Add a `pycraf.pathprof.SrtmConf` option ``void_fill`` (``'zero'``, + ``'nan'`` or ``'interp'``) to control how void pixels are handled. + +Bugfixes +-------- +pycraf.pathprof +^^^^^^^^^^^^^^ +- Mask the canonical SRTM void sentinel ``-32768`` when reading ``.hgt`` + tiles. Previously only ``-32767``/``+32767`` were masked, so genuine + ``-32768`` voids leaked through and were linearly blended with valid + neighbours, producing spurious deep pits in height profiles. + 2.1.0 (2025-02-02) ======================= This is a maintenance release. pycraf now uses numpy version 2.0.0 or higher. diff --git a/docs/pathprof/working_with_srtm.rst b/docs/pathprof/working_with_srtm.rst index 83c748a3..0a9374fd 100644 --- a/docs/pathprof/working_with_srtm.rst +++ b/docs/pathprof/working_with_srtm.rst @@ -145,10 +145,87 @@ We refer to `~scipy.interpolate.RectBivariateSpline` description for further information. +Copernicus DEM (GLO-90 / GLO-30) +================================ + +Besides the SRTM *.hgt* tiles, `~pycraf.pathprof` can use the `Copernicus DEM +`_ +(GLO-90 and GLO-30) as a terrain source. Compared to SRTM this has a number of +advantages: + +- It is hosted as `Cloud-Optimised GeoTIFFs + `_ on the AWS Open Data + buckets (no authentication needed) and thus - unlike the retired SRTM + download servers - provides a reliably working automatic-download path. +- It is global (pole-to-pole), whereas SRTM only covers 60 deg N to 56 deg S. +- It is void-free over water (edited/flattened water bodies), avoiding the + data-gap artifacts that SRTM tiles can contain. + +To use it, select one of the two Copernicus servers:: + + >>> from pycraf.pathprof import SrtmConf + >>> SrtmConf.set( # doctest: +SKIP + ... srtm_dir='/path/to/copernicus', + ... download='missing', + ... server='copernicus_glo90', # or 'copernicus_glo30' + ... ) + +The public API (`~pycraf.pathprof.srtm_height_data`, +`~pycraf.pathprof.srtm_height_profile`, and the P.452 `~pycraf.pathprof.PathProp` +engine) is unchanged; internally the reader accounts for the Copernicus tiles' +pixel-centre (area) registration and their latitude-dependent longitude spacing +(the tiles are not square above :math:`|\mathrm{lat}| = 50` deg). + +.. note:: + + Reading the GeoTIFF tiles requires the optional `rasterio + `_ package. If a Copernicus server is + selected without it, a clear error is raised. + +.. note:: + + The Copernicus DEM uses the EGM2008 geoid as its vertical datum (SRTM uses + EGM96). The difference is at the metre level and is not converted, since + propagation profiles only depend on *relative* terrain heights. + +**Attribution.** When using or redistributing Copernicus DEM data, the +following statement must be reproduced: + + Produced using Copernicus WorldDEM-90 © DLR e.V. 2010-2014 and + © Airbus Defence and Space GmbH 2014-2018 provided under COPERNICUS by + the European Union and ESA; all rights reserved. + + +Handling missing tiles and voids +================================ + +Two `~pycraf.pathprof.SrtmConf` options control what happens when data is +missing or invalid. By default, if a tile that *should* exist on the chosen +server is not found on disk (and is not downloaded), its terrain is set to +zero and a `~pycraf.pathprof.TileNotAvailableOnDiskWarning` is emitted. To turn +this into a hard error instead (so that computations never silently run over +zero terrain), use:: + + >>> SrtmConf.set(on_missing='raise') # doctest: +IGNORE_OUTPUT + +Void pixels (SRTM data gaps; the Copernicus DEM is void-free) are, by default, +replaced with zero. They can instead be kept as ``NaN`` (to be detected +downstream) or filled from the nearest valid pixels:: + + >>> SrtmConf.set(void_fill='interp') # doctest: +IGNORE_OUTPUT + +.. note:: + + Restore the historic defaults with + ``SrtmConf.set(on_missing='zeros', void_fill='zero')``. + + Download links ============== - `NASA v2.1 `__ - `NASA v1.0 `__ - `viewfinderpanoramas.org `__ +- `Copernicus GLO-90 (AWS Open Data) `__ +- `Copernicus GLO-30 (AWS Open Data) `__ diff --git a/pycraf/pathprof/srtm.py b/pycraf/pathprof/srtm.py index e1dbea63..b3cf273c 100644 --- a/pycraf/pathprof/srtm.py +++ b/pycraf/pathprof/srtm.py @@ -18,6 +18,46 @@ V4.1 and viewfinderpanoramas forbid commercial use (without explicit permission). + +Copernicus DEM +-------------- + +In addition to the SRTM '.hgt' tiles, `pycraf` can use the Copernicus DEM +(GLO-90 and GLO-30) as a terrain source. These are hosted as Cloud-Optimised +GeoTIFFs on the AWS Open Data buckets (no authentication required) and, unlike +the retired SRTM auto-download servers, provide a reliable download path. The +Copernicus DEM is global (pole-to-pole, whereas SRTM only covers 60N to 56S) +and is void-free over water. To use it:: + + from pycraf.pathprof import SrtmConf + SrtmConf.set( + srtm_dir='/path/to/copernicus', download='missing', + server='copernicus_glo90', + ) + +Reading the GeoTIFF tiles requires the (optional) `rasterio` package. + +Model type: the Copernicus DEM is a Digital Surface Model (DSM), i.e. it +"represents the surface of the Earth including buildings, infrastructure and +vegetation" (source: Copernicus DEM readme, +https://copernicus-dem-30m.s3.amazonaws.com/readme.html; see also the +Copernicus DEM Product Handbook). It is derived from the TanDEM-X radar +mission (X-band). It is therefore NOT a bare-earth terrain model; canopy and +building heights are included. This matches SRTM, which is also a radar DSM +(C-band), so the surface-vs-terrain character is unchanged when switching +between the two. Note that P.452 additionally models clutter separately (via +the `zone_t` / `zone_r` options), so assigning clutter over a DSM in built-up +or forested terminals can double-count vegetation/building height. + +Attribution (required when using or redistributing Copernicus DEM data): + + Produced using Copernicus WorldDEM-90 (c) DLR e.V. 2010-2014 and + (c) Airbus Defence and Space GmbH 2014-2018 provided under COPERNICUS by + the European Union and ESA; all rights reserved. + +Note, the Copernicus DEM uses the EGM2008 geoid as vertical datum (SRTM uses +EGM96); the difference is at the metre level and is not converted here, as +propagation profiles only depend on relative terrain heights. ''' @@ -59,6 +99,30 @@ VIEWPANO_TILES = np.load(_VIEWPANO_NAME) +# Copernicus DEM (GLO-90 / GLO-30) on the AWS Open Data buckets. +# The "code" is the internal arc-second*10/3 identifier that appears in the +# tile names ('30' -> GLO-90 = 3 arcsec, '10' -> GLO-30 = 1 arcsec); it is +# *not* the resolution in metres. `hgt_res` is the nominal (equatorial) +# latitude pixel spacing in metres, used to pick the height-profile sampling. +COPERNICUS_SERVERS = { + 'copernicus_glo90': { + 'base_url': 'https://copernicus-dem-90m.s3.amazonaws.com/', + 'code': '30', + 'hgt_res': 90., + }, + 'copernicus_glo30': { + 'base_url': 'https://copernicus-dem-30m.s3.amazonaws.com/', + 'code': '10', + 'hgt_res': 30., + }, + } + +# cache of the authoritative tile inventories (server -> set of tile names), +# lazily populated from a local copy of the bucket "tileList.txt" or, if +# downloading is enabled, from the bucket root +_COPERNICUS_TILE_LISTS = {} + + class TileNotAvailableOnServerError(Exception): pass @@ -109,9 +173,31 @@ class SrtmConf(utils.MultiState): The default behavior is to not download anything (`download='never'`). There is even an option, to always force download (`download='always'`). - The default download server will be `server='nasa_v2.1'`. One could - also use the (very old) data (`server='nasa_v1.0'`) or inofficial - tiles from viewfinderpanorama (`server='viewpano'`). + The default download server is `server='viewpano'` (inofficial tiles + from viewfinderpanorama). Alternatively, one can use the Copernicus DEM + (`server='copernicus_glo90'` or `server='copernicus_glo30'`), which - + unlike the retired SRTM servers - offers a reliably working download + path (AWS Open Data, no authentication), is global (pole-to-pole) and + void-free over water. The Copernicus tiles are Cloud-Optimised GeoTIFFs + and require the optional `rasterio` package (a clear error is raised if + a Copernicus server is selected without it). See the module + documentation for the required attribution statement. + + Two further options control the behaviour when data is missing or + invalid. `on_missing` decides what happens if a tile that *should* + exist on the chosen server is not found on disk (and cannot be + downloaded): `on_missing='zeros'` (default) sets the terrain of that + tile to zero and emits a `TileNotAvailableOnDiskWarning` (the historic + behaviour), whereas `on_missing='raise'` raises a + `TileNotAvailableOnDiskError` instead - useful to avoid silently + computing over zero-terrain. `void_fill` controls how void pixels + (SRTM data gaps; the Copernicus DEM is void-free) are treated: + `void_fill='zero'` (default) replaces voids with zero, + `void_fill='nan'` keeps them as `NaN` (so they can be detected + downstream), and `void_fill='interp'` fills them from the nearest + valid pixels. To change these use:: + + SrtmConf.set(on_missing='raise', void_fill='interp') Of course, one can set several of these options simultaneously:: @@ -156,7 +242,7 @@ class SrtmConf(utils.MultiState): _attributes = ( 'srtm_dir', 'download', 'server', 'interp', 'spline_opts', - 'tile_size', 'hgt_res' + 'on_missing', 'void_fill', 'tile_size', 'hgt_res' ) srtm_dir = os.environ.get('SRTMDATA', '.') @@ -164,6 +250,8 @@ class SrtmConf(utils.MultiState): server = 'viewpano' interp = 'linear' spline_opts = (3, 0) + on_missing = 'zeros' + void_fill = 'zero' tile_size = 1201 hgt_res = 90. # m; basic SRTM resolution (refers to 3 arcsec resolution) @@ -195,10 +283,28 @@ def validate(cls, **kwargs): 'are supported for "download" option.' ) if k == 'server': - if v not in ['viewpano']: + if v not in [ + 'viewpano', + 'copernicus_glo90', 'copernicus_glo30', + ]: raise ValueError( - 'Only the value "viewpano" is currently ' - 'supported for "server" option.' + 'Only the values "viewpano", "copernicus_glo90", ' + 'and "copernicus_glo30" are currently supported for ' + 'the "server" option.' + ) + + if k == 'on_missing': + if v not in ['zeros', 'raise']: + raise ValueError( + 'Only the values "zeros" and "raise" are supported ' + 'for the "on_missing" option.' + ) + + if k == 'void_fill': + if v not in ['zero', 'nan', 'interp']: + raise ValueError( + 'Only the values "zero", "nan", and "interp" are ' + 'supported for the "void_fill" option.' ) if k == 'interp': @@ -244,6 +350,12 @@ def hook(cls, **kwargs): # check if srtm_dir changed and clear cache if kwargs['srtm_dir'] != cls.srtm_dir: get_tile_interpolator.cache_clear() + # the Copernicus tile inventory is per directory + _COPERNICUS_TILE_LISTS.clear() + + if 'server' in kwargs and kwargs['server'] != cls.server: + # the cached Copernicus inventory is per server + _COPERNICUS_TILE_LISTS.clear() if 'download' in kwargs: # check if 'download' strategy was changed and clear cache @@ -259,22 +371,35 @@ def hook(cls, **kwargs): if kwargs['server'] != cls.server: get_tile_interpolator.cache_clear() + if 'on_missing' in kwargs: + # changes whether a missing tile becomes zeros or an error, + # i.e. the cached tile data would differ + if kwargs['on_missing'] != cls.on_missing: + get_tile_interpolator.cache_clear() + + if 'void_fill' in kwargs: + # changes how void pixels are filled in the cached interpolator + if kwargs['void_fill'] != cls.void_fill: + get_tile_interpolator.cache_clear() + @classmethod def __repr__(cls): return ( ''.format( + 'interp: {}, spline_opts: {}, on_missing: {}, void_fill: {}>' + ''.format( cls.srtm_dir, cls.download, cls.server, - cls.interp, cls.spline_opts + cls.interp, cls.spline_opts, cls.on_missing, cls.void_fill )) @classmethod def __str__(cls): return ( 'SrtmConf\n directory: {}\n download: {}\n server: {}\n' - ' interp: {}\n spline_opts: {}'.format( + ' interp: {}\n spline_opts: {}\n on_missing: {}\n' + ' void_fill: {}'.format( cls.srtm_dir, cls.download, cls.server, - cls.interp, cls.spline_opts + cls.interp, cls.spline_opts, cls.on_missing, cls.void_fill )) @@ -289,6 +414,58 @@ def _hgt_filename(ilon, ilat): ) +def _copernicus_tilename(ilon, ilat, server=None): + # construct the Copernicus DEM tile (base) name for the tile whose + # south-west corner is at the integer degree (ilon, ilat) + + if server is None: + server = SrtmConf.server + code = COPERNICUS_SERVERS[server]['code'] + + return ( + 'Copernicus_DSM_COG_{code}_{ns:1s}{ilat:02d}_00_{ew:1s}{ilon:03d}' + '_00_DEM'.format( + code=code, + ns='N' if ilat >= 0 else 'S', ilat=abs(ilat), + ew='E' if ilon >= 0 else 'W', ilon=abs(ilon), + ) + ) + + +def _copernicus_tile_set(server=None): + # return the authoritative set of tile names for a Copernicus server, + # or None if the inventory is not available (and cannot be fetched) + + if server is None: + server = SrtmConf.server + + if server in _COPERNICUS_TILE_LISTS: + return _COPERNICUS_TILE_LISTS[server] + + srtm_dir = SrtmConf.srtm_dir + list_name = os.path.join(srtm_dir, server + '_tileList.txt') + + tiles = None + if os.path.exists(list_name): + with open(list_name, 'r') as f: + tiles = set(line.strip() for line in f if line.strip()) + elif SrtmConf.download in ['missing', 'always']: + # fetch the authoritative list from the bucket root and cache it + base_url = COPERNICUS_SERVERS[server]['base_url'] + tmp_path = download_file(base_url + 'tileList.txt') + try: + os.makedirs(srtm_dir, exist_ok=True) + shutil.copyfile(tmp_path, list_name) + except OSError: + # srtm_dir not writable - keep the in-memory copy only + pass + with open(tmp_path, 'r') as f: + tiles = set(line.strip() for line in f if line.strip()) + + _COPERNICUS_TILE_LISTS[server] = tiles + return tiles + + def _check_availability(ilon, ilat): # check availability of a tile on download servers # returns continent name (for NASA server) or zip file name (Pano) @@ -296,7 +473,23 @@ def _check_availability(ilon, ilat): server = SrtmConf.server tile_name = _hgt_filename(ilon, ilat) - if server.startswith('nasa_v'): + if server.startswith('copernicus'): + + cop_name = _copernicus_tilename(ilon, ilat) + tiles = _copernicus_tile_set() + + # if the inventory is unknown (download='never' and no cached list), + # we cannot rule the tile out - defer the decision to the disk lookup + if tiles is not None and cop_name not in tiles: + raise TileNotAvailableOnServerError( + 'No tile found for ({}d, {}d) in list of available ' + 'tiles.'.format( + ilon, ilat + )) + + return cop_name + + elif server.startswith('nasa_v'): for continent, tiles in NASA_TILES.items(): if tile_name in tiles: @@ -354,6 +547,21 @@ def _download(ilon, ilat): srtm_dir = SrtmConf.srtm_dir server = SrtmConf.server + if server.startswith('copernicus'): + + # Copernicus tiles are single Cloud-Optimised GeoTIFFs, stored as + # "/.tif" in the bucket; we keep them flat on disk + base_url = COPERNICUS_SERVERS[server]['base_url'] + cop_name = _copernicus_tilename(ilon, ilat) + full_url = base_url + cop_name + '/' + cop_name + '.tif' + tile_path = os.path.join(srtm_dir, cop_name + '.tif') + + tmp_path = download_file(full_url) + os.makedirs(srtm_dir, exist_ok=True) + shutil.move(tmp_path, tile_path) + + return + tile_name = _hgt_filename(ilon, ilat) tile_path = os.path.join(srtm_dir, tile_name) @@ -476,67 +684,221 @@ def get_hgt_file(ilon, ilat): return hgt_file -def get_tile_data(ilon, ilat): - # angles in deg +def get_copernicus_file(ilon, ilat): + # locate (and, if requested, download) the Copernicus GeoTIFF tile whose + # south-west corner is at the integer degree (ilon, ilat) + + _check_availability(ilon, ilat) + + srtm_dir = SrtmConf.srtm_dir + tif_name = _copernicus_tilename(ilon, ilat) + '.tif' + tif_file = _get_hgt_diskpath(tif_name) + + download = SrtmConf.download + if download == 'always' or (tif_file is None and download == 'missing'): + + _download(ilon, ilat) + + tif_file = _get_hgt_diskpath(tif_name) + if tif_file is None: + raise TileNotAvailableOnDiskError( + 'No Copernicus tile found for ({}d, {}d), was looking for {}\n' + 'in directory: {}'.format( + ilon, ilat, tif_name, srtm_dir + )) + + return tif_file + + +# metres per degree of latitude (mean, spherical Earth); only used to turn the +# tile's latitude pixel spacing into an approximate resolution for choosing the +# height-profile sampling step +_M_PER_DEG_LAT = 111120. + + +def _read_copernicus_cog(tif_file): + # read a Copernicus DEM Cloud-Optimised GeoTIFF and return coordinate and + # height arrays in pycraf's tile convention: + # lons -> shape (nlon, 1), ascending west -> east + # lats -> shape (1, nlat), ascending south -> north + # tile -> shape (nlat, nlon), tile[lat_idx, lon_idx], voids set to NaN + # Copernicus tiles use pixel-centre (area) registration and their + # longitude spacing widens above |lat| 50 deg, so the actual geotransform + # is read from the file rather than assumed. try: - hgt_file = get_hgt_file(ilon, ilat) - # need to run check after get_hgt_file, because download could happen - _check_consistent_tile_sizes(SrtmConf.srtm_dir) - tile = np.fromfile(hgt_file, dtype='>i2') - tile_size = int(np.sqrt(tile.size) + 0.5) - hgt_res = 90. * 1200 / (tile_size - 1) - SrtmConf.set(tile_size=tile_size, _do_validate=False) - SrtmConf.set(hgt_res=hgt_res, _do_validate=False) - tile = tile.reshape((tile_size, tile_size))[::-1] - - bad_mask = (tile == 32767) | (tile == -32767) - tile = tile.astype(np.float32) - tile[bad_mask] = np.nan + import rasterio + except ImportError as e: + raise ImportError( + 'The "rasterio" package is required to read Copernicus DEM ' + '(GeoTIFF) tiles. Install it (e.g. "pip install rasterio") or ' + 'select a different "server" in pycraf.pathprof.SrtmConf.' + ) from e + + with rasterio.open(tif_file) as ds: + tile = ds.read(1).astype(np.float32) # (nlat, nlon), row 0 = north + transf = ds.transform + nodata = ds.nodata + + nlat, nlon = tile.shape + # pixel-centre coordinates from the affine transform + lon0 = transf.c + 0.5 * transf.a # centre of column 0 + lat0 = transf.f + 0.5 * transf.e # centre of row 0 (northern-most) + lon_axis = lon0 + np.arange(nlon) * transf.a # west -> east + lat_axis = lat0 + np.arange(nlat) * transf.e # north -> south + + # flip rows to go south -> north, matching the '.hgt' convention + tile = tile[::-1] + lat_axis = lat_axis[::-1] + + # NoData handling: Copernicus is void-free over water, but coastal tiles + # can carry NaN or negative sentinels; mask conservatively (sea -> 0 m is + # applied later via "void_fill", consistent with the SRTM convention) + bad_mask = ~np.isfinite(tile) | (tile < -500.) + if nodata is not None and np.isfinite(nodata): + bad_mask |= (tile == np.float32(nodata)) + tile[bad_mask] = np.nan + + hgt_res = abs(transf.e) * _M_PER_DEG_LAT + SrtmConf.set(tile_size=nlat, _do_validate=False) + SrtmConf.set(hgt_res=hgt_res, _do_validate=False) + + lons = lon_axis[:, np.newaxis] + lats = lat_axis[np.newaxis, :] + return lons, lats, tile - except TileNotAvailableOnServerError: - # always use very small tile size for zero tiles - # (just enough to make spline interpolation work) - tile_size = 5 - tile = np.zeros((tile_size, tile_size), dtype=np.float32) - except TileNotAvailableOnDiskError: - # also set to zero, but raise a warning - tile_size = 5 - tile = np.zeros((tile_size, tile_size), dtype=np.float32) - - tile_name = _hgt_filename(ilon, ilat) - srtm_dir = SrtmConf.srtm_dir - warnings.warn( - ''' -No hgt-file found for ({}d, {}d) - was looking for file {} +def _missing_tile_warning(ilon, ilat, tile_name): + # emit the historic "tile not on disk -> zeros" warning + + srtm_dir = SrtmConf.srtm_dir + warnings.warn( + ''' +No tile found for ({}d, {}d) - was looking for file {} in directory: {} Will set terrain heights in this area to zero. Note, you can have pycraf download missing tiles automatically - just use "pycraf.pathprof.SrtmConf" -(see its documentation).'''.format(ilon, ilat, tile_name, srtm_dir), - category=TileNotAvailableOnDiskWarning, - stacklevel=1, - ) +(see its documentation). To turn this into an error instead, set +"SrtmConf.set(on_missing='raise')".'''.format( + ilon, ilat, tile_name, srtm_dir), + category=TileNotAvailableOnDiskWarning, + stacklevel=1, + ) + +def _zero_tile(ilon, ilat): + # a minimal (5x5) zero tile, just big enough for spline interpolation + tile_size = 5 + tile = np.zeros((tile_size, tile_size), dtype=np.float32) dx = dy = 1. / (tile_size - 1) x, y = np.ogrid[0:tile_size, 0:tile_size] lons, lats = x * dx + ilon, y * dy + ilat return lons, lats, tile +def _get_srtm_tile_data(ilon, ilat): + + hgt_file = get_hgt_file(ilon, ilat) + # need to run check after get_hgt_file, because download could happen + _check_consistent_tile_sizes(SrtmConf.srtm_dir) + tile = np.fromfile(hgt_file, dtype='>i2') + tile_size = int(np.sqrt(tile.size) + 0.5) + hgt_res = 90. * 1200 / (tile_size - 1) + SrtmConf.set(tile_size=tile_size, _do_validate=False) + SrtmConf.set(hgt_res=hgt_res, _do_validate=False) + tile = tile.reshape((tile_size, tile_size))[::-1] + + # void/NoData sentinels: the canonical SRTM void is -32768 (0x8000); + # -32767 and +32767 are also seen in some products. (The historic code + # masked only -32767/+32767, so genuine -32768 voids leaked through and + # were linearly blended with valid neighbours, producing spurious pits.) + bad_mask = (tile == -32768) | (tile == -32767) | (tile == 32767) + tile = tile.astype(np.float32) + tile[bad_mask] = np.nan + + dx = dy = 1. / (tile_size - 1) + x, y = np.ogrid[0:tile_size, 0:tile_size] + lons, lats = x * dx + ilon, y * dy + ilat + return lons, lats, tile + + +def get_tile_data(ilon, ilat): + # angles in deg + + server = SrtmConf.server + + try: + if server.startswith('copernicus'): + tif_file = get_copernicus_file(ilon, ilat) + return _read_copernicus_cog(tif_file) + else: + return _get_srtm_tile_data(ilon, ilat) + + except TileNotAvailableOnServerError: + # tile is genuinely absent from the server (ocean, polar cap): use a + # zero tile silently, as before + return _zero_tile(ilon, ilat) + + except TileNotAvailableOnDiskError: + # tile should exist but is not on disk (and wasn't downloaded); either + # raise or fall back to zeros (+ warning), depending on "on_missing" + if SrtmConf.on_missing == 'raise': + raise + + if server.startswith('copernicus'): + tile_name = _copernicus_tilename(ilon, ilat) + '.tif' + else: + tile_name = _hgt_filename(ilon, ilat) + _missing_tile_warning(ilon, ilat, tile_name) + return _zero_tile(ilon, ilat) + + +def _fill_voids(tile, void_fill): + # resolve NaN void pixels in a tile according to the "void_fill" policy; + # returns a float array free of NaNs unless void_fill == 'nan' + + if void_fill == 'nan': + return tile + + bad_mask = ~np.isfinite(tile) + if not bad_mask.any(): + return tile + + if void_fill == 'interp': + if bad_mask.all(): + return np.nan_to_num(tile) + from scipy import ndimage + # fill each void with the value of the nearest valid pixel + idx = ndimage.distance_transform_edt( + bad_mask, return_distances=False, return_indices=True + ) + return tile[tuple(idx)] + + # void_fill == 'zero' + return np.nan_to_num(tile) + + # cannot use SrtmConf inside to query interp and spline_opts, because -# caching might cause problems +# caching might cause problems (changing them does not clear this cache). +# "void_fill" is safe to query here, because the SrtmConf.hook clears this +# cache whenever "void_fill" (or the tile source) changes. @lru_cache(maxsize=36, typed=False) def get_tile_interpolator(ilon, ilat, interp, spline_opts): # angles in deg lons, lats, tile = get_tile_data(ilon, ilat) - # have to treat NaNs in some way; set to zero for now - tile = np.nan_to_num(tile) + # resolve voids (NaNs); default replaces them with zero + tile = _fill_voids(tile, SrtmConf.void_fill) if interp in ['nearest', 'linear']: + # bounds_error=False + fill_value=None extrapolates for query points + # that fall just outside the pixel-centre grid. This is needed for the + # Copernicus (area-registered) tiles, whose pixel centres do not reach + # the southern/eastern tile edge, and is a no-op for the (node- + # registered) SRTM tiles, which always cover the assigned degree cell. _tile_interpolator = RegularGridInterpolator( (lons[:, 0], lats[0]), tile.T, method=interp, + bounds_error=False, fill_value=None, ) elif interp == 'spline': kx = ky = spline_opts[0] diff --git a/pycraf/pathprof/tests/test_srtm.py b/pycraf/pathprof/tests/test_srtm.py index 63d47dc9..47848a39 100644 --- a/pycraf/pathprof/tests/test_srtm.py +++ b/pycraf/pathprof/tests/test_srtm.py @@ -2,6 +2,7 @@ # -*- coding: utf-8 -*- import os +import importlib import pytest import numpy as np from numpy.testing import assert_equal, assert_allclose @@ -13,6 +14,12 @@ TOL_KWARGS = {'atol': 1.e-4, 'rtol': 1.e-4} +# skip Copernicus (GeoTIFF) reading tests, if rasterio is not installed +skip_rasterio = pytest.mark.skipif( + importlib.util.find_spec('rasterio') is None, + reason='"rasterio" package not installed' + ) + class TestSrtmConf: @@ -570,3 +577,377 @@ def test_srtm_height_data_broadcasting(srtm_temp_dir): [[433.44000244, 416.20001221, 704.52001953, 826.08001709], [358.72000122, 395.55999756, 263.83999634, 469.39999390]] ]) * apu.m) + + +# --------------------------------------------------------------------------- +# Failure-mode hardening (missing tiles, void masking) +# --------------------------------------------------------------------------- + +def _write_hgt(path, tile_size, fill, void_positions=(), void_value=-32768): + '''Write a minimal big-endian int16 SRTM ".hgt" tile.''' + + arr = np.full((tile_size, tile_size), fill, dtype='>i2') + for pos in void_positions: + arr[pos] = void_value + arr.astype('>i2').tofile(path) + + +@pytest.mark.parametrize('void_value', [-32768, -32767, 32767]) +def test_srtm_void_is_masked(srtm_temp_dir, void_value): + # regression test: the canonical SRTM void sentinel is -32768; make sure + # all of -32768/-32767/+32767 are masked (the historic code only handled + # -32767/+32767, so genuine -32768 voids leaked through as huge values) + + tdir = os.path.join(srtm_temp_dir, 'voidtest_{}'.format(void_value)) + os.makedirs(tdir, exist_ok=True) + # N50E006 (ilon=6, ilat=50) is a valid viewpano tile + _write_hgt( + os.path.join(tdir, 'N50E006.hgt'), 6, 100, + void_positions=[(2, 3)], void_value=void_value, + ) + + with srtm.SrtmConf.set(srtm_dir=tdir, server='viewpano'): + srtm.get_tile_interpolator.cache_clear() + _, _, tile = srtm.get_tile_data(6, 50) + + assert np.isnan(tile).sum() == 1 + # no sentinel value survived into the (float) tile + assert np.nanmin(tile) == 100 and np.nanmax(tile) == 100 + + +@pytest.mark.parametrize('on_missing', ['zeros', 'raise']) +def test_on_missing_behaviour(srtm_temp_dir, on_missing): + # a tile that is "available on the server" but not on disk, with + # download='never', either warns+zeros or raises, depending on on_missing + + tdir = os.path.join(srtm_temp_dir, 'missing_{}'.format(on_missing)) + os.makedirs(tdir, exist_ok=True) + + # 15E, 47N is not a viewpano tile -> would be a *server* miss (silent + # zeros); use 6E, 50N which *is* a viewpano tile -> a *disk* miss + with srtm.SrtmConf.set( + srtm_dir=tdir, server='viewpano', + download='never', on_missing=on_missing, + ): + srtm.get_tile_interpolator.cache_clear() + + if on_missing == 'raise': + with pytest.raises(srtm.TileNotAvailableOnDiskError): + srtm.get_tile_data(6, 50) + else: + with pytest.warns(srtm.TileNotAvailableOnDiskWarning): + _, _, tile = srtm.get_tile_data(6, 50) + assert_allclose(tile, np.zeros((5, 5), dtype=np.float32)) + + +# --------------------------------------------------------------------------- +# Copernicus DEM backend +# --------------------------------------------------------------------------- + +def test_copernicus_tilename(): + + cases = [ + (6, 45, 'copernicus_glo90', + 'Copernicus_DSM_COG_30_N45_00_E006_00_DEM'), + (109, -31, 'copernicus_glo90', + 'Copernicus_DSM_COG_30_S31_00_E109_00_DEM'), + (-97, 51, 'copernicus_glo90', + 'Copernicus_DSM_COG_30_N51_00_W097_00_DEM'), + (-46, -22, 'copernicus_glo90', + 'Copernicus_DSM_COG_30_S22_00_W046_00_DEM'), + (6, 45, 'copernicus_glo30', + 'Copernicus_DSM_COG_10_N45_00_E006_00_DEM'), + ] + + for ilon, ilat, server, name in cases: + assert srtm._copernicus_tilename(ilon, ilat, server) == name + + +def _write_cop_tile( + path, ilon, ilat, nlon, nlat, dx_asec, dy_asec, values, nodata=None + ): + '''Write a synthetic Copernicus-style COG (pixel-centre registration). + + The pixel *centre* of the north-west pixel is placed exactly at + (ilon, ilat + 1), matching the real Copernicus tiles (named by their + south-west corner). + ''' + + import rasterio + + dx = dx_asec / 3600. + dy = dy_asec / 3600. + west_edge = ilon - 0.5 * dx + north_edge = (ilat + 1) + 0.5 * dy + transform = rasterio.Affine(dx, 0.0, west_edge, 0.0, -dy, north_edge) + + kwargs = dict( + driver='GTiff', height=nlat, width=nlon, count=1, + dtype='float32', crs='EPSG:4326', transform=transform, + ) + if nodata is not None: + kwargs['nodata'] = nodata + + with rasterio.open(path, 'w', **kwargs) as dst: + dst.write(values.astype('float32'), 1) + + +@skip_rasterio +def test_copernicus_registration(srtm_temp_dir): + + tdir = os.path.join(srtm_temp_dir, 'cop_reg') + os.makedirs(tdir, exist_ok=True) + + # coarse full-degree tile at SW corner (6, 45); values encode (row, col) + nlon = nlat = 12 + dasec = 300. # 5 arcmin -> 12 px per degree + rows, cols = np.mgrid[0:nlat, 0:nlon] + values = 1000. + 10. * rows + cols # row 0 = north + _write_cop_tile( + os.path.join(tdir, 'Copernicus_DSM_COG_30_N45_00_E006_00_DEM.tif'), + 6, 45, nlon, nlat, dasec, dasec, values, + ) + + with srtm.SrtmConf.set( + srtm_dir=tdir, server='copernicus_glo90', + download='never', interp='nearest', + ): + srtm.get_tile_interpolator.cache_clear() + + lons, lats, tile = srtm.get_tile_data(6, 45) + + # pixel-centre registration: NW pixel centre at (6, 46) + assert tile.shape == (nlat, nlon) + assert_allclose(lons[0, 0], 6.0) # west-most centre + assert_allclose(lats[0, -1], 46.0) # north-most centre + # hgt_res inferred from the tile (5 arcmin ~ 9.3 km here) + assert srtm.SrtmConf.hgt_res > 0 + + # query an exact pixel centre (stored row 3, col 5 -> value 1035) + dx = dasec / 3600. + lon_q = 6.0 + 5 * dx + lat_q = 46.0 - 3 * dx + h = srtm.srtm_height_data([lon_q] * apu.deg, [lat_q] * apu.deg) + assert_quantity_allclose(h, [1035.] * apu.m) + + +@skip_rasterio +def test_copernicus_variable_spacing(srtm_temp_dir): + # above |lat| 50 deg the longitude spacing widens; tiles are not square + # and the reader must use the actual geotransform + + tdir = os.path.join(srtm_temp_dir, 'cop_var') + os.makedirs(tdir, exist_ok=True) + + nlon, nlat = 4, 12 + dx_asec, dy_asec = 900., 300. # 9 arcsec would be a GLO-90 70-75 deg tile + values = np.arange(nlon * nlat, dtype='float32').reshape(nlat, nlon) + _write_cop_tile( + os.path.join(tdir, 'Copernicus_DSM_COG_30_N70_00_E020_00_DEM.tif'), + 20, 70, nlon, nlat, dx_asec, dy_asec, values, + ) + + with srtm.SrtmConf.set( + srtm_dir=tdir, server='copernicus_glo90', download='never', + ): + srtm.get_tile_interpolator.cache_clear() + + lons, lats, tile = srtm.get_tile_data(20, 70) + + assert tile.shape == (nlat, nlon) + assert lons.shape == (nlon, 1) + assert lats.shape == (1, nlat) + # longitude spacing != latitude spacing + assert_allclose((lons[1, 0] - lons[0, 0]) * 3600., dx_asec) + assert_allclose((lats[0, 1] - lats[0, 0]) * 3600., dy_asec) + + +@skip_rasterio +@pytest.mark.parametrize('void_fill', ['zero', 'nan', 'interp']) +def test_copernicus_void_fill(srtm_temp_dir, void_fill): + + tdir = os.path.join(srtm_temp_dir, 'cop_void_{}'.format(void_fill)) + os.makedirs(tdir, exist_ok=True) + + nlon = nlat = 12 + values = np.full((nlat, nlon), 200., dtype='float32') + values[5, 5] = np.nan # NaN void + values[6, 6] = -1000. # sea/void sentinel (< -500) + _write_cop_tile( + os.path.join(tdir, 'Copernicus_DSM_COG_30_N45_00_E006_00_DEM.tif'), + 6, 45, nlon, nlat, 300., 300., values, + ) + + with srtm.SrtmConf.set( + srtm_dir=tdir, server='copernicus_glo90', + download='never', interp='nearest', void_fill=void_fill, + ): + srtm.get_tile_interpolator.cache_clear() + + _, _, tile = srtm.get_tile_data(6, 45) + # both the NaN and the < -500 sentinel are masked to NaN on read + assert np.isnan(tile).sum() == 2 + + dx = 300. / 3600. + # stored row 5, col 5 -> the NaN void + lon_q = 6.0 + 5 * dx + lat_q = 46.0 - 5 * dx + h = srtm.srtm_height_data( + [lon_q] * apu.deg, [lat_q] * apu.deg + ).to_value(apu.m) + + if void_fill == 'zero': + assert_allclose(h, [0.]) + elif void_fill == 'nan': + assert np.isnan(h).all() + elif void_fill == 'interp': + # filled from the nearest valid pixel (200 m) + assert_allclose(h, [200.]) + + +@skip_rasterio +def test_copernicus_check_availability(srtm_temp_dir): + # with a cached tileList, tiles not in the list are treated as ocean + # (silent zeros); tiles in the list but missing on disk hit on_missing + + tdir = os.path.join(srtm_temp_dir, 'cop_avail') + os.makedirs(tdir, exist_ok=True) + with open( + os.path.join(tdir, 'copernicus_glo90_tileList.txt'), 'w' + ) as f: + f.write('Copernicus_DSM_COG_30_N45_00_E006_00_DEM\n') + + with srtm.SrtmConf.set( + srtm_dir=tdir, server='copernicus_glo90', download='never', + ): + srtm.get_tile_interpolator.cache_clear() + + # not in the list -> ocean -> silent zeros + _, _, tile = srtm.get_tile_data(0, 0) + assert_allclose(tile, np.zeros((5, 5), dtype=np.float32)) + + # in the list but not on disk -> disk miss + with srtm.SrtmConf.set(on_missing='raise'): + with pytest.raises(srtm.TileNotAvailableOnDiskError): + srtm.get_tile_data(6, 45) + + +@pytest.mark.remote_data +@skip_rasterio +def test_copernicus_download(srtm_temp_dir): + + tdir = os.path.join(srtm_temp_dir, 'cop_dl') + os.makedirs(tdir, exist_ok=True) + + ilon, ilat = 6, 45 + tif_name = srtm._copernicus_tilename( + ilon, ilat, 'copernicus_glo90' + ) + '.tif' + + with srtm.SrtmConf.set( + srtm_dir=tdir, server='copernicus_glo90', download='missing', + ): + srtm.get_tile_interpolator.cache_clear() + + dl_path = srtm.get_copernicus_file(ilon, ilat) + assert dl_path is not None + assert dl_path.endswith(tif_name) + + +@pytest.mark.remote_data +@skip_rasterio +def test_copernicus_get_tile_data(srtm_temp_dir): + # value check against a real downloaded GLO-90 tile (Mont Blanc region), + # analogous to test_get_tile_data for SRTM. Note the pixel-centre (area) + # registration: the south-most row centre lies half a pixel plus the + # removed shared edge above the SW corner (45.00083), the north-most row + # centre exactly at 46.0. + + tdir = os.path.join(srtm_temp_dir, 'cop_tiledata') + os.makedirs(tdir, exist_ok=True) + + with srtm.SrtmConf.set( + srtm_dir=tdir, server='copernicus_glo90', download='missing', + ): + srtm.get_tile_interpolator.cache_clear() + + ilon, ilat = 6, 45 + lons, lats, tile = srtm.get_tile_data(ilon, ilat) + + assert tile.shape == (1200, 1200) + + assert_allclose(lons[::250, 0], np.array([ + 6., 6.20833333, 6.41666667, 6.625, 6.83333333 + ])) + assert_allclose(lats[0, ::250], np.array([ + 45.00083333, 45.20916667, 45.4175, 45.62583333, 45.83416667 + ])) + assert_allclose(tile[::250, ::250], np.array([ + [2443.5603, 3120.2942, 2775.819, 2283.7593, 1530.1218], + [2173.671, 2129.5164, 1736.1337, 1493.1702, 1649.1984], + [248.32872, 2015.355, 2559.966, 1861.7087, 2698.8232], + [1093.4408, 1177.1719, 1019.66125, 2303.7588, 1472.1299], + [370., 445.5, 1134.9672, 1501.6011, 3701.4497] + ]), rtol=1.e-6) + + +@pytest.mark.remote_data +@skip_rasterio +def test_copernicus_registration_real(srtm_temp_dir): + # Mont Blanc summit (45.8326 N, 6.8652 E); GLO-90 reads ~4790 m there + # (a wrong half-pixel registration would miss the narrow summit) + + tdir = os.path.join(srtm_temp_dir, 'cop_real') + os.makedirs(tdir, exist_ok=True) + + with srtm.SrtmConf.set( + srtm_dir=tdir, server='copernicus_glo90', + download='missing', interp='linear', + ): + srtm.get_tile_interpolator.cache_clear() + + h = srtm.srtm_height_data( + [6.8652] * apu.deg, [45.8326] * apu.deg + ).to_value(apu.m)[0] + + assert 4700. < h < 4810. + + +@pytest.mark.remote_data +@skip_rasterio +def test_copernicus_vs_srtm_overlap(srtm_temp_dir): + # cross-check GLO-90 against the viewpano SRTM tile on an overlapping + # land cell (French Alps). The two independent DEMs must be highly + # correlated; we use the correlation coefficient rather than an absolute + # tolerance, so the result is insensitive to the EGM96/EGM2008 datum + # offset (a few metres) between them. + + lons, lats = np.meshgrid( + np.arange(6.2, 6.8, 0.02), + np.arange(45.2, 45.8, 0.02), + ) + lons = lons.flatten() * apu.deg + lats = lats.flatten() * apu.deg + + tdir_s = os.path.join(srtm_temp_dir, 'xcheck_srtm') + tdir_c = os.path.join(srtm_temp_dir, 'xcheck_cop') + os.makedirs(tdir_s, exist_ok=True) + os.makedirs(tdir_c, exist_ok=True) + + with srtm.SrtmConf.set( + srtm_dir=tdir_s, server='viewpano', download='missing', + ): + srtm.get_tile_interpolator.cache_clear() + h_srtm = srtm.srtm_height_data(lons, lats).to_value(apu.m) + + with srtm.SrtmConf.set( + srtm_dir=tdir_c, server='copernicus_glo90', download='missing', + ): + srtm.get_tile_interpolator.cache_clear() + h_cop = srtm.srtm_height_data(lons, lats).to_value(apu.m) + + corr = np.corrcoef(h_srtm, h_cop)[0, 1] + assert corr > 0.99 + # and the median absolute difference should be small (metre-to-decametre + # level over this relief) + assert np.median(np.abs(h_srtm - h_cop)) < 40.