diff --git a/opendm/tiles/gdal2tiles.py b/opendm/tiles/gdal2tiles.py index 081c335a5..f34377034 100644 --- a/opendm/tiles/gdal2tiles.py +++ b/opendm/tiles/gdal2tiles.py @@ -51,6 +51,8 @@ from osgeo import gdal from osgeo import osr +_EPSG_SRS_CACHE = {} + try: from PIL import Image import numpy @@ -209,7 +211,6 @@ def __init__(self, tileSize=256): self.initialResolution = 2 * math.pi * 6378137 / self.tileSize # 156543.03392804062 for tileSize 256 pixels self.originShift = 2 * math.pi * 6378137 / 2.0 - # 20037508.342789244 def LatLonToMeters(self, lat, lon): "Converts given lat/lon in WGS84 Datum to XY in Spherical Mercator EPSG:3857" @@ -223,10 +224,11 @@ def LatLonToMeters(self, lat, lon): def MetersToLatLon(self, mx, my): "Converts XY point from Spherical Mercator EPSG:3857 to lat/lon in WGS84 Datum" - lon = (mx / self.originShift) * 180.0 - lat = (my / self.originShift) * 180.0 - - lat = 180 / math.pi * (2 * math.atan(math.exp(lat * math.pi / 180.0)) - math.pi / 2.0) + inv_originShift = 180.0 / self.originShift + lon = mx * inv_originShift + lat = my * inv_originShift + pi = math.pi + lat = 180.0 / pi * (2.0 * math.atan(math.exp(lat * pi / 180.0)) - pi / 2.0) return lat, lon def PixelsToMeters(self, px, py, zoom): @@ -413,7 +415,7 @@ def TileBounds(self, tx, ty, zoom): ) def TileLatLonBounds(self, tx, ty, zoom): - "Returns bounds of the given tile in the SWNE form" + """Returns bounds of the given tile in the SWNE form""" b = self.TileBounds(tx, ty, zoom) return (b[1], b[0], b[3], b[2]) @@ -712,12 +714,14 @@ def setup_output_srs(input_srs, options): """ Setup the desired SRS (based on options) """ - output_srs = osr.SpatialReference() - - if options.profile == 'mercator': - output_srs.ImportFromEPSG(3857) - elif options.profile == 'geodetic': - output_srs.ImportFromEPSG(4326) + profile = options.profile + + # Use cached instances for standard profiles; else copy input_srs + if profile == 'mercator': + # Use a clone to avoid potential mutation of cached SRS + output_srs = _get_epsg_srs(3857).Clone() + elif profile == 'geodetic': + output_srs = _get_epsg_srs(4326).Clone() else: output_srs = input_srs @@ -2825,27 +2829,61 @@ def log_progress(self, nb_items=1): def get_tile_swne(tile_job_info, options): + # Optimize by avoiding redundant class constructions and SRS parsing if called repeatedly with same profile if options.profile == 'mercator': - mercator = GlobalMercator() + # GlobalMercator is stateless except for __init__ parameters (tileSize) + # If pattern of repeated tile calls is observed with default tileSize, we can use a static instance + # Create an attribute on the function for persistent reuse if possible (thread-safety not considered) + if not hasattr(get_tile_swne, "_global_mercator") or getattr(get_tile_swne, "_gm_tileSize", None) != 256: + get_tile_swne._global_mercator = GlobalMercator() + get_tile_swne._gm_tileSize = 256 + mercator = get_tile_swne._global_mercator tile_swne = mercator.TileLatLonBounds + elif options.profile == 'geodetic': - geodetic = GlobalGeodetic(options.tmscompatible) + tmscompatible = options.tmscompatible + if not hasattr(get_tile_swne, "_global_geodetic") or getattr(get_tile_swne, "_gg_tmscompatible", object()) != tmscompatible: + get_tile_swne._global_geodetic = GlobalGeodetic(tmscompatible) + get_tile_swne._gg_tmscompatible = tmscompatible + geodetic = get_tile_swne._global_geodetic tile_swne = geodetic.TileLatLonBounds + elif options.profile == 'raster': - srs4326 = osr.SpatialReference() - srs4326.ImportFromEPSG(4326) + # These OSR operations are very costly. Cache SRS/CT if possible, based on in_srs_wkt + srs4326 = getattr(get_tile_swne, "_cached_srs4326", None) + if srs4326 is None: + srs4326 = osr.SpatialReference() + srs4326.ImportFromEPSG(4326) + get_tile_swne._cached_srs4326 = srs4326 + if tile_job_info.kml and tile_job_info.in_srs_wkt: - in_srs = osr.SpatialReference() - in_srs.ImportFromWkt(tile_job_info.in_srs_wkt) - ct = osr.CoordinateTransformation(in_srs, srs4326) + cached_key = (tile_job_info.in_srs_wkt, id(srs4326)) + cache = getattr(get_tile_swne, "_ct_cache", None) + if cache is None: + cache = {} + get_tile_swne._ct_cache = cache + if cached_key in cache: + ct = cache[cached_key][0] + in_srs = cache[cached_key][1] + else: + in_srs = osr.SpatialReference() + in_srs.ImportFromWkt(tile_job_info.in_srs_wkt) + ct = osr.CoordinateTransformation(in_srs, srs4326) + cache[cached_key] = (ct, in_srs) # Avoid further lookups + + tilesize = tile_job_info.tilesize + tmaxz = tile_job_info.tmaxz + out_geo_trans = tile_job_info.out_geo_trans + ominy = tile_job_info.ominy + is_epsg_4326 = tile_job_info.is_epsg_4326 def rastertileswne(x, y, z): - pixelsizex = (2 ** (tile_job_info.tmaxz - z) * tile_job_info.out_geo_trans[1]) - west = tile_job_info.out_geo_trans[0] + x * tile_job_info.tilesize * pixelsizex - east = west + tile_job_info.tilesize * pixelsizex - south = tile_job_info.ominy + y * tile_job_info.tilesize * pixelsizex - north = south + tile_job_info.tilesize * pixelsizex - if not tile_job_info.is_epsg_4326: + pixelsizex = (2 ** (tmaxz - z) * out_geo_trans[1]) + west = out_geo_trans[0] + x * tilesize * pixelsizex + east = west + tilesize * pixelsizex + south = ominy + y * tilesize * pixelsizex + north = south + tilesize * pixelsizex + if not is_epsg_4326: # Transformation to EPSG:4326 (WGS84 datum) west, south = ct.TransformPoint(west, south)[:2] east, north = ct.TransformPoint(east, north)[:2] @@ -2942,6 +2980,14 @@ def main(): else: multi_threaded_tiling(input_file, output_folder, options) +def _get_epsg_srs(epsg_code): + srs = _EPSG_SRS_CACHE.get(epsg_code) + if srs is None: + srs = osr.SpatialReference() + srs.ImportFromEPSG(epsg_code) + _EPSG_SRS_CACHE[epsg_code] = srs + return srs + if __name__ == '__main__': main()