From 903453f2d85ecb7ea7315b71a49c71b1bdb14d53 Mon Sep 17 00:00:00 2001 From: Tariq Soliman Date: Mon, 2 Dec 2024 17:52:35 -0800 Subject: [PATCH] Add base_every_zoom to gdal2customtiles --- .../gdal2customtiles/gdal2customtiles.py | 110 +++++++++++++++--- .../gdal2customtiles/rasters2customtiles.py | 42 ++++++- 2 files changed, 136 insertions(+), 16 deletions(-) diff --git a/auxiliary/gdal2customtiles/gdal2customtiles.py b/auxiliary/gdal2customtiles/gdal2customtiles.py index 9adeb3de..8a96125c 100644 --- a/auxiliary/gdal2customtiles/gdal2customtiles.py +++ b/auxiliary/gdal2customtiles/gdal2customtiles.py @@ -56,9 +56,13 @@ from xml.etree import ElementTree from osgeo import gdal, osr +gdal.UseExceptions() Options = Any +# Unfortunately modifying this breaks other things +PLANET_RADIUS = 6378137 + # MMGIS def binary(num): return ''.join(bin(c).replace('0b', '').rjust(8, '0') for c in struct.pack('!f', num)) @@ -252,7 +256,6 @@ def ZoomForPixelSize(self, pixelSize, overriden_tile_size): def PixelsToMeters(self, px, py, zoom, overriden_tile_size): "Converts pixel coordinates in given zoom level of pyramid to EPSG:3857" - res = self.resolution * self.tile_size / \ overriden_tile_size / (2**zoom) mx = px * res + self.topleft_x @@ -511,9 +514,9 @@ class GlobalMercator(object): def __init__(self, tile_size: int = 256) -> None: "Initialize the TMS Global Mercator pyramid" self.tile_size = tile_size - self.initialResolution = 2 * math.pi * 6378137 / self.tile_size + self.initialResolution = 2 * math.pi * PLANET_RADIUS / self.tile_size # 156543.03392804062 for tile_size 256 pixels - self.originShift = 2 * math.pi * 6378137 / 2.0 + self.originShift = 2 * math.pi * PLANET_RADIUS / 2.0 # 20037508.342789244 def LatLonToMeters(self, lat, lon): @@ -541,7 +544,6 @@ def MetersToLatLon(self, mx, my): def PixelsToMeters(self, px, py, zoom): "Converts pixel coordinates in given zoom level of pyramid to EPSG:3857" - res = self.Resolution(zoom) mx = px * res - self.originShift my = py * res - self.originShift @@ -596,7 +598,7 @@ def TileLatLonBounds(self, tx, ty, zoom): def Resolution(self, zoom): "Resolution (meters/pixel) for given zoom level (measured at Equator)" - # return (2 * math.pi * 6378137) / (self.tile_size * 2**zoom) + # return (2 * math.pi * PLANET_RADIUS) / (self.tile_size * 2**zoom) return self.initialResolution / (2**zoom) def ZoomForPixelSize(self, pixelSize): @@ -1855,6 +1857,12 @@ def optparse_init() -> optparse.OptionParser: choices=resampling_list, help="Resampling method (%s) - default 'average'" % ",".join(resampling_list), ) + p.add_option( + "--base_every_zoom", + dest="base_every_zoom", + action="store_true", + help="Instead of making base tiles and then overview tiles from them, generate every zoom level as if it was the base zoom level. This may prevent inevitable rounding errors from deeper zoom levels from accumulating and propagating to lower zoom levels.", + ) p.add_option( "-s", "--s_srs", @@ -2366,7 +2374,6 @@ def open_input(self) -> None: ) # Open the input file - if self.input_file: input_dataset: gdal.Dataset = gdal.Open( self.input_file, gdal.GA_ReadOnly) @@ -3062,7 +3069,7 @@ def generate_base_tiles(self) -> Tuple[TileJobInfo, List[TileDetail]]: """ if not self.options.quiet: - print("Generating Base Tiles:") + print("Generating Base Tiles at zoom " + str(self.tmaxz) + ":") if self.options.verbose: print("") @@ -3144,6 +3151,7 @@ def generate_base_tiles(self) -> Tuple[TileJobInfo, List[TileDetail]]: rx, ry, rxsize, rysize = rb wx, wy, wxsize, wysize = wb + # print(tx, ty, tz, b, rb, wb, querysize) # MMGIS elif self.isRasterBounded: # 'raster' profile: @@ -3295,7 +3303,7 @@ def geo_query(self, ds, ulx, uly, lrx, lry, querysize=0): ry = int((uly - geotran[3]) / geotran[5] + 0.001) rxsize = max(1, int((lrx - ulx) / geotran[1] + 0.5)) rysize = max(1, int((lry - uly) / geotran[5] + 0.5)) - + if not querysize: wxsize, wysize = rxsize, rysize else: @@ -4743,14 +4751,79 @@ def getYTile(ty, tz, options): return ty +def ZoomForPixelSize(pixelSize): + MAXZOOMLEVEL = 32 + tile_size = 256 + initialResolution = 2 * math.pi * PLANET_RADIUS / tile_size + for i in range(MAXZOOMLEVEL): + if pixelSize > initialResolution / (2**i): + return max(0, i) # We don't want to scale up + return MAXZOOMLEVEL def worker_tile_details( - input_file: str, output_folder: str, options: Options + input_file: str, output_folder: str, options: Options, pool ) -> Tuple[TileJobInfo, List[TileDetail]]: - gdal2tiles = GDAL2Tiles(input_file, output_folder, options) - gdal2tiles.open_input() - gdal2tiles.generate_metadata() - tile_job_info, tile_details = gdal2tiles.generate_base_tiles() + + tile_job_info = [] + tile_details = [] + + if options.base_every_zoom: + + print("Tiling each zoom level as if it is the base zoom level...") + + minZoom = options.zoom[0] + maxZoom = options.zoom[1] + if minZoom is None: + if maxZoom is None: + minZoom = 0 + else: + minZoom = maxZoom + if maxZoom is None: + ds = gdal.Open(input_file) + gt = ds.GetGeoTransform() + maxZoom = ZoomForPixelSize(gt[1]) + + print("Tiling zoom " + str(minZoom) + " - " + str(maxZoom) + "...") + + firstPass = True + for base_tz in range(minZoom, maxZoom + 1, 1): + options.zoom[1] = base_tz + gdal2tiles = GDAL2Tiles(input_file, output_folder, options) + gdal2tiles.open_input() + if firstPass: + gdal2tiles.generate_metadata() + + tile_job_info_bez, tile_details_bez = gdal2tiles.generate_base_tiles() + + if not options.verbose and not options.quiet: + base_progress_bar = ProgressBar(len(tile_details_bez)) + base_progress_bar.start() + + if pool: + nb_processes = options.nb_processes or 1 + chunksize = max(1, min(128, len(tile_details) // nb_processes)) + for _ in pool.imap_unordered( + partial(create_base_tile, tile_job_info_bez), tile_details_bez, chunksize=chunksize + ): + if not options.verbose and not options.quiet: + base_progress_bar.log_progress() + else: + for tile_detail in tile_details_bez: + create_base_tile(tile_job_info_bez, tile_detail) + + if not options.verbose and not options.quiet: + base_progress_bar.log_progress() + + if getattr(threadLocal, "cached_ds", None): + del threadLocal.cached_ds + + firstPass = True + else: + gdal2tiles = GDAL2Tiles(input_file, output_folder, options) + gdal2tiles.open_input() + gdal2tiles.generate_metadata() + tile_job_info, tile_details = gdal2tiles.generate_base_tiles() + return tile_job_info, tile_details @@ -4852,6 +4925,10 @@ def single_threaded_tiling( conf, tile_details = worker_tile_details( input_file, output_folder, options) + if options.base_every_zoom: + sys.exit('Done!') + return + if options.verbose: print("Tiles details calc complete.") @@ -4900,8 +4977,13 @@ def multi_threaded_tiling( print("Begin tiles details calc") conf, tile_details = worker_tile_details( - input_file, output_folder, options) + input_file, output_folder, options, pool) + if options.base_every_zoom: + shutil.rmtree(os.path.dirname(conf.src_file)) + sys.exit('Done!') + return + if options.verbose: print("Tiles details calc complete.") diff --git a/auxiliary/gdal2customtiles/rasters2customtiles.py b/auxiliary/gdal2customtiles/rasters2customtiles.py index d89444d9..79661811 100644 --- a/auxiliary/gdal2customtiles/rasters2customtiles.py +++ b/auxiliary/gdal2customtiles/rasters2customtiles.py @@ -2,7 +2,9 @@ import subprocess import optparse from osgeo import gdal, osr +import math +PLANET_RADIUS = 3396190 # 6378137 def optparse_init() -> optparse.OptionParser: """Prepare the option parser for input (argv)""" @@ -48,6 +50,14 @@ def optparse_init() -> optparse.OptionParser: metavar="NODATA", help="Value in the input dataset considered as transparent", ) + + p.add_option( + "-d", + "--dbez", + dest="dbez", + action="store_true", + help="Don't base every zoom. If true, instead of tiling each zoom level as the base zoom level, builds overview tiles from the four deeper zoom tiles from which it is made.", + ) return p @@ -95,7 +105,6 @@ def ReprojectCoords(coords, src_srs, tgt_srs): trans_coords.append([x, y]) return trans_coords - def AutoGdalTranslate(geo_extent, cols, rows, raster): gdal_translate = "gdal_translate -of VRT -a_srs EPSG:4326 -gcp 0 0 " + str(geo_extent[0][0]) + " " + str(geo_extent[0][1]) + " -gcp " + str(cols) + " 0 " + str(geo_extent[3][0]) + " " + str( geo_extent[3][1]) + " -gcp " + str(cols) + " " + str(rows) + " " + str(geo_extent[2][0]) + " " + str(geo_extent[2][1]) + " " + raster + " " + raster[:-4] + ".vrt" @@ -122,13 +131,24 @@ def AutoGdal2Tiles(raster, options, outputdir): srcnodata = " --srcnodata=0,0,0" if options.srcnodata is not None: srcnodata = f" --srcnodata={options.srcnodata}" + base_every_zoom = " --base_every_zoom" + if options.dbez is True: + base_every_zoom = "" output = "" if outputdir is not None: output = f" {outputdir}" - gdal2tiles = f"python gdal2customtiles.py -n{dem}{processes}{tilesize}{zoom}{resume}{srcnodata} {raster[:-4]}.vrt{output}" + gdal2tiles = f"python gdal2customtiles.py -n{dem}{processes}{tilesize}{zoom}{resume}{srcnodata}{base_every_zoom} {raster[:-4]}.vrt{output}" print(f"Running: {gdal2tiles}\n") subprocess.Popen(gdal2tiles) +def ZoomForPixelSize(pixelSize): + MAXZOOMLEVEL = 32 + tile_size = 256 + initialResolution = 2 * math.pi * PLANET_RADIUS / tile_size + for i in range(MAXZOOMLEVEL): + if pixelSize > initialResolution / (2**i): + return max(0, i - 1) # We don't want to scale up + return MAXZOOMLEVEL - 1 parser = optparse_init() options, args = parser.parse_args(args=sys.argv) @@ -139,8 +159,26 @@ def AutoGdal2Tiles(raster, options, outputdir): gt = ds.GetGeoTransform() cols = ds.RasterXSize rows = ds.RasterYSize + +""" +x_res = gt[1] +y_res = -gt[5] + +# oversample if needed +oversample_zooms = ZoomForPixelSize(gt[1]) +print(oversample_zooms) + +if options.zoom: + oversample_factor = 1 << (19 - oversample_zooms) + x_res = x_res / oversample_factor + y_res = y_res / oversample_factor + cols = cols * oversample_factor + rows = rows * oversample_factor +print(cols, ds.RasterXSize, rows, ds.RasterYSize) +""" extent = GetExtent(gt, cols, rows) + src_srs = osr.SpatialReference() src_srs.ImportFromWkt(ds.GetProjection()) tgt_srs = src_srs.CloneGeogCS()