94 lines
4.7 KiB
Python
94 lines
4.7 KiB
Python
"""Prepare the two downloaded Copernicus COGs; no browser/runtime GIS dependency.
|
|
|
|
python -m pip install --target .tools/terrain -r scripts/terrain-requirements.txt
|
|
python scripts/prepare-terrain.py
|
|
"""
|
|
import hashlib
|
|
import json
|
|
from pathlib import Path
|
|
import subprocess
|
|
import sys
|
|
|
|
ROOT = Path(__file__).resolve().parents[1]
|
|
sys.path.insert(0, str(ROOT / '.tools' / 'terrain'))
|
|
import numpy as np
|
|
import rasterio
|
|
from rasterio.merge import merge
|
|
from rasterio.warp import reproject, Resampling
|
|
from rasterio.transform import from_origin
|
|
from pyproj import CRS, Transformer
|
|
|
|
RADIUS = 12000
|
|
SPACING = 30
|
|
SIZE = 2 * RADIUS // SPACING + 1
|
|
ATTRIBUTION = ('produced using Copernicus WorldDEM-30 © 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')
|
|
|
|
|
|
def main():
|
|
# The app remains the single source of observer coordinates.
|
|
result = subprocess.run(['node', '--input-type=module', '-e',
|
|
"import {observer} from './src/config/observer.ts'; console.log(JSON.stringify(observer))"],
|
|
cwd=ROOT, capture_output=True, text=True, check=True)
|
|
observer = json.loads(result.stdout)
|
|
lat, lon = observer['latitude'], observer['longitude']
|
|
projection = f'+proj=aeqd +lat_0={lat} +lon_0={lon} +datum=WGS84 +units=m +no_defs'
|
|
local_crs = CRS.from_proj4(projection)
|
|
to_geo = Transformer.from_crs(local_crs, 'EPSG:4326', always_xy=True)
|
|
transform = from_origin(-RADIUS - SPACING / 2, RADIUS + SPACING / 2, SPACING, SPACING)
|
|
output = ROOT / 'public' / 'terrain'
|
|
output.mkdir(parents=True, exist_ok=True)
|
|
sources, manifest = [], []
|
|
try:
|
|
for east in (126, 127):
|
|
name = f'Copernicus_DSM_COG_10_N37_00_E{east}_00_DEM'
|
|
path = ROOT / 'data' / 'raw' / 'copernicus' / f'{name}.tif'
|
|
dataset = rasterio.open(path)
|
|
sources.append(dataset)
|
|
if dataset.crs.to_epsg() != 4326 or dataset.shape != (3600, 3600):
|
|
raise ValueError(f'Unexpected source grid: {path}')
|
|
manifest.append({
|
|
'file': path.name, 'url': f'https://copernicus-dem-30m.s3.amazonaws.com/{name}/{name}.tif',
|
|
'bytes': path.stat().st_size, 'sha256': hashlib.sha256(path.read_bytes()).hexdigest(),
|
|
'crs': dataset.crs.to_string(), 'bounds': list(dataset.bounds),
|
|
'shape': list(dataset.shape), 'pixelSizeDegrees': list(dataset.res),
|
|
'nodata': dataset.nodata,
|
|
})
|
|
# Merge before resampling, so bilinear interpolation spans the 127° seam.
|
|
mosaic, mosaic_transform = merge(sources, nodata=-9999)
|
|
heights = np.full((SIZE, SIZE), -9999, dtype=np.float32)
|
|
reproject(mosaic[0], heights, src_transform=mosaic_transform, src_crs=sources[0].crs,
|
|
src_nodata=-9999, dst_transform=transform, dst_crs=local_crs,
|
|
dst_nodata=-9999, resampling=Resampling.bilinear)
|
|
missing = int(np.count_nonzero(~np.isfinite(heights) | (heights == -9999)))
|
|
if missing:
|
|
raise ValueError(f'Refusing to publish a grid with {missing} missing values')
|
|
heights.astype('<f4').tofile(output / 'gwanaksan.f32')
|
|
bounds = [to_geo.transform(x, y) for x, y in
|
|
[(-RADIUS, -RADIUS), (RADIUS, -RADIUS), (RADIUS, RADIUS), (-RADIUS, RADIUS)]]
|
|
metadata = {
|
|
'version': 1, 'dataset': 'Copernicus DEM GLO-30 Public (2021 release)',
|
|
'surfaceType': 'DSM', 'sourceResolution': '1 arc-second (~30m)',
|
|
'verticalDatum': 'EGM2008 orthometric height (m)',
|
|
'origin': {'latitude': lat, 'longitude': lon},
|
|
'projection': projection, 'projectionWkt': local_crs.to_wkt(), 'width': SIZE, 'height': SIZE,
|
|
'spacingMeters': SPACING, 'radiusMeters': RADIUS,
|
|
'layout': 'row-major north-to-south, west-to-east, float32 little-endian',
|
|
'heightFile': 'gwanaksan.f32', 'geographicCorners': bounds,
|
|
'observerSurfaceElevation': float(heights[SIZE // 2, SIZE // 2]),
|
|
'minElevation': float(heights.min()), 'maxElevation': float(heights.max()),
|
|
'missingSamples': missing, 'sources': manifest, 'attribution': ATTRIBUTION,
|
|
}
|
|
(output / 'gwanaksan.json').write_text(json.dumps(metadata, ensure_ascii=False, indent=2) + '\n', encoding='utf-8')
|
|
print(json.dumps({k: metadata[k] for k in ['width', 'height', 'observerSurfaceElevation', 'minElevation', 'maxElevation', 'missingSamples']}, indent=2))
|
|
for item in manifest:
|
|
print(item['file'], item['sha256'])
|
|
finally:
|
|
for dataset in sources:
|
|
dataset.close()
|
|
|
|
|
|
if __name__ == '__main__':
|
|
main()
|