From 960c2899e1001c25b7a3fd97839d02ff29929850 Mon Sep 17 00:00:00 2001 From: Andrew Kirmse Date: Mon, 11 Nov 2024 16:18:28 -0800 Subject: [PATCH] Use full projection string for gdalbuildvrt. Sometimes two projections have the same name, but different internal parameters, which would cause gdalbuiltvrt to fail because it can handle only one projection in a VRT file. So now we use the full WKT string to describe the projection. --- scripts/run_prominence.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/scripts/run_prominence.py b/scripts/run_prominence.py index f5c3958..0f75d75 100644 --- a/scripts/run_prominence.py +++ b/scripts/run_prominence.py @@ -44,9 +44,6 @@ def round_up(coord, interval): """Return coord rounded up to the nearest interval""" return math.ceil(coord / interval) * interval -def sign(a): - return bool(a > 0) - bool(a < 0) - def filename_for_coordinates(x, y, degrees_per_tile): """Return output filename for the given coordinates""" y += degrees_per_tile # Name uses upper left corner @@ -122,7 +119,11 @@ def create_vrts(tile_dir, input_files, skip_boundary): ds = gdal.Open(input_file) prj = ds.GetProjection() srs = osr.SpatialReference(wkt=prj) - projection_name = srs.GetAttrValue('projcs') + # gdalbuildvrt needs the projections to match pretty much exactly within a VRT file. + # It's not enough to use the projection name; sometimes two projections have the + # same name but slightly different datums, and gdalbuildvrt will fail. + # So instead we use the whole WKT string. + projection_name = srs.ExportToWkt() projection_map.setdefault(projection_name, []).append(input_file) # For each projection, build raw and warped VRTs