From f1abb634cbe48561855f9448577823ef57789fe5 Mon Sep 17 00:00:00 2001
From: Barry A Bragg Jr <hokieg3n1us@gmail.com>
Date: Fri, 15 Oct 2021 21:16:48 -0400
Subject: [PATCH 01/12] Implement MapboxTilesRenderer

Create MapboxTilesRenderer for outputting TMS tile sets as a mbtiles file (a SQLite database following the MBTiles 1.3 specification).
---
 datashader/tiles.py | 108 +++++++++++++++++++++++++++++++++++++---
 datashader/utils.py | 119 +++++++++++++++++++++++++++-----------------
 2 files changed, 174 insertions(+), 53 deletions(-)

diff --git a/datashader/tiles.py b/datashader/tiles.py
index 920cccd3a..ba813237d 100644
--- a/datashader/tiles.py
+++ b/datashader/tiles.py
@@ -1,16 +1,17 @@
 from __future__ import absolute_import, division, print_function
-from io import BytesIO
 
 import math
 import os
+import sqlite3
+from io import BytesIO
 
 import dask
 import dask.bag as db
-
 import numpy as np
-
 from PIL.Image import fromarray
 
+from .utils import meters_to_lnglat
+
 __all__ = ['render_tiles', 'MercatorTileDefinition']
 
 
@@ -60,8 +61,11 @@ def calculate_zoom_level_stats(super_tiles, load_data_func,
 
 def render_tiles(full_extent, levels, load_data_func,
                  rasterize_func, shader_func,
-                 post_render_func, output_path, color_ranging_strategy='fullscan'):
+                 post_render_func, output_path, color_ranging_strategy='fullscan', num_workers=4):
     results = dict()
+
+    validate_output_path(output_path, full_extent, levels)
+
     for level in levels:
         print('calculating statistics for level {}'.format(level))
         super_tiles, span = calculate_zoom_level_stats(list(gen_super_tiles(full_extent, level)),
@@ -69,7 +73,7 @@ def render_tiles(full_extent, levels, load_data_func,
                                                        color_ranging_strategy=color_ranging_strategy)
         print('rendering {} supertiles for zoom level {} with span={}'.format(len(super_tiles), level, span))
         b = db.from_sequence(super_tiles)
-        b.map(render_super_tile, span, output_path, shader_func, post_render_func).compute()
+        b.map(render_super_tile, span, output_path, shader_func, post_render_func).compute(num_workers=num_workers)
         results[level] = dict(success=True, stats=span, supertile_count=len(super_tiles))
 
     return results
@@ -98,17 +102,25 @@ def render_super_tile(tile_info, span, output_path, shader_func, post_render_fun
     return create_sub_tiles(ds_img, level, tile_info, output_path, post_render_func)
 
 
-def create_sub_tiles(data_array, level, tile_info, output_path, post_render_func=None):
+def validate_output_path(output_path, full_extent, levels):
     # validate / createoutput_dir
-    _create_dir(output_path)
+    if os.path.isdir(output_path):
+        _create_dir(output_path)
+    elif output_path.endswith("mbtiles"):
+        MapboxTileRenderer.setup(output_path, full_extent, levels[0], levels[len(levels) - 1])
 
+
+def create_sub_tiles(data_array, level, tile_info, output_path, post_render_func=None):
     # create tile source
     tile_def = MercatorTileDefinition(x_range=tile_info['x_range'],
                                       y_range=tile_info['y_range'],
                                       tile_size=256)
 
     # create Tile Renderer
-    if output_path.startswith('s3:'):
+    if output_path.endswith("mbtiles"):
+        renderer = MapboxTileRenderer(tile_def, output_location=output_path,
+                                      post_render_func=post_render_func)
+    elif output_path.startswith('s3:'):
         renderer = S3TileRenderer(tile_def, output_location=output_path,
                                   post_render_func=post_render_func)
     else:
@@ -407,3 +419,83 @@ def render(self, da, level):
             client.put_object(Body=output_buf, Bucket=bucket, Key=key, ACL='public-read')
 
         return 'https://{}.s3.amazonaws.com/{}'.format(bucket, s3_info.path)
+
+
+class MapboxTileRenderer(TileRenderer):
+
+    @staticmethod
+    def setup(output_location, full_extent, min_zoom, max_zoom, tile_format="PNG"):
+        con = sqlite3.connect(output_location)
+        cur = con.cursor()
+
+        # Create MBTiles tables.
+        cur.execute("""
+            create table tiles (
+                zoom_level integer,
+                tile_column integer,
+                tile_row integer,
+                tile_data blob);
+                """)
+        cur.execute("""create table metadata
+            (name text, value text);""")
+        cur.execute("""create unique index name on metadata (name);""")
+        cur.execute("""create unique index tile_index on tiles
+            (zoom_level, tile_column, tile_row);""")
+
+        # Compute Extents in WGS84
+        min_lon, min_lat = meters_to_lnglat(full_extent[0], full_extent[1])
+        max_lon, max_lat = meters_to_lnglat(full_extent[2], full_extent[3])
+
+        # Compute the Center & Zoom Level
+        lat_diff = max_lat - min_lat
+        lon_diff = max_lon - min_lon
+
+        lat_center = min_lat + (lat_diff / 2)
+        lon_center = min_lon + (lon_diff / 2)
+
+        max_diff = max(lon_diff, lat_diff)
+        zoom_level = None
+        if max_diff < (360.0 / math.pow(2, 20)):
+            zoom_level = 21
+        else:
+            zoom_level = int(-1 * ((math.log(max_diff) / math.log(2.0)) - (math.log(360.0) / math.log(2))))
+            if zoom_level < 1:
+                zoom_level = 1
+
+        # Setup MBTiles metadata table.
+        cur.execute("""insert into metadata (name, value) values (?, ?) """,
+                    ("name", "tiles"))
+        cur.execute("""insert into metadata (name, value) values (?, ?) """,
+                    ("format", tile_format.lower()))
+        cur.execute("""insert into metadata (name, value) values (?, ?) """,
+                    ("bounds ", "{},{},{},{}".format(min_lon, min_lat, max_lon, max_lat)))
+        cur.execute("""insert into metadata (name, value) values (?, ?) """,
+                    ("center ", "{},{},{}".format(lon_center, lat_center, zoom_level)))
+        cur.execute("""insert into metadata (name, value) values (?, ?) """,
+                    ("minzoom ", min_zoom))
+        cur.execute("""insert into metadata (name, value) values (?, ?) """,
+                    ("maxzoom ", max_zoom))
+
+        cur.close()
+        con.commit()
+        con.close()
+
+    def render(self, da, level):
+        con = sqlite3.connect(self.output_location, isolation_level=None)
+        cur = con.cursor()
+
+        for img, x, y, z in super(MapboxTileRenderer, self).render(da, level):
+            image_bytes = BytesIO()
+            img.save(image_bytes, self.tile_format)
+            image_bytes.seek(0)
+
+            tile_row = (2 ** z) - 1 - y
+
+            cur.execute("""insert into tiles (zoom_level,
+                tile_column, tile_row, tile_data) values
+                (?, ?, ?, ?);""",
+                        (z, x, tile_row, sqlite3.Binary(image_bytes.getvalue())))
+
+        cur.close()
+        con.commit()
+        con.close()
diff --git a/datashader/utils.py b/datashader/utils.py
index 2d90ffbf6..a91097ce2 100644
--- a/datashader/utils.py
+++ b/datashader/utils.py
@@ -45,8 +45,8 @@ class VisibleDeprecationWarning(UserWarning):
 
 # Get and save the Numba version, will be used to limit functionality
 numba_version = tuple([int(x) for x in re.match(
-                            r"([0-9]+)\.([0-9]+)\.([0-9]+)",
-                            nb.__version__).groups()])
+    r"([0-9]+)\.([0-9]+)\.([0-9]+)",
+    nb.__version__).groups()])
 
 
 class Expr(object):
@@ -56,6 +56,7 @@ class Expr(object):
     ``inputs`` attribute/property, containing a tuple of everything that fully
     defines that expression.
     """
+
     def __hash__(self):
         return hash((type(self), self._hashable_inputs()))
 
@@ -85,6 +86,7 @@ def _hashable_inputs(self):
 
 class Dispatcher(object):
     """Simple single dispatch."""
+
     def __init__(self):
         self._lookup = {}
 
@@ -212,17 +214,17 @@ def calc_bbox(xs, ys, res):
 
     xmin = ymin = np.inf
     xmax = ymax = -np.inf
-    Ab = np.array([[res[0],  0.,      xbound],
-                   [0.,      -res[1], ybound],
-                   [0.,      0.,      1.]])
+    Ab = np.array([[res[0], 0., xbound],
+                   [0., -res[1], ybound],
+                   [0., 0., 1.]])
     for x_, y_ in [(0, 0), (0, len(ys)), (len(xs), 0), (len(xs), len(ys))]:
         x, y, _ = np.dot(Ab, np.array([x_, y_, 1.]))
         if x < xmin: xmin = x
         if x > xmax: xmax = x
         if y < ymin: ymin = y
         if y > ymax: ymax = y
-    xpad, ypad = res[0]/2., res[1]/2.
-    return xmin-xpad, ymin+ypad, xmax-xpad, ymax+ypad
+    xpad, ypad = res[0] / 2., res[1] / 2.
+    return xmin - xpad, ymin + ypad, xmax - xpad, ymax + ypad
 
 
 def get_indices(start, end, coords, res):
@@ -241,11 +243,11 @@ def get_indices(start, end, coords, res):
         Resolution along an axis (aka "grid/cell sizes")
     """
     size = len(coords)
-    half = abs(res)/2.
+    half = abs(res) / 2.
     vmin, vmax = coords.min(), coords.max()
-    span = vmax-vmin
-    start, end = start+half-vmin, end-half-vmin
-    sidx, eidx = int((start/span)*size), int((end/span)*size)
+    span = vmax - vmin
+    start, end = start + half - vmin, end - half - vmin
+    sidx, eidx = int((start / span) * size), int((end / span) * size)
     if eidx < sidx:
         return sidx, sidx
     return sidx, eidx
@@ -274,7 +276,7 @@ def orient_array(raster, res=None, layer=None):
     if res is None:
         res = calc_res(raster)
     array = raster.data
-    if layer is not None: array = array[layer-1]
+    if layer is not None: array = array[layer - 1]
     r0zero = np.timedelta64(0, 'ns') if isinstance(res[0], np.timedelta64) else 0
     r1zero = np.timedelta64(0, 'ns') if isinstance(res[1], np.timedelta64) else 0
     if array.ndim == 2:
@@ -312,11 +314,11 @@ def compute_coords(width, height, x_range, y_range, res):
         1D array of y-coordinates
     """
     (x0, x1), (y0, y1) = x_range, y_range
-    xd = (x1-x0)/float(width)
-    yd = (y1-y0)/float(height)
-    xpad, ypad = abs(xd/2.), abs(yd/2.)
-    x0, x1 = x0+xpad, x1-xpad
-    y0, y1 = y0+ypad, y1-ypad
+    xd = (x1 - x0) / float(width)
+    yd = (y1 - y0) / float(height)
+    xpad, ypad = abs(xd / 2.), abs(yd / 2.)
+    x0, x1 = x0 + xpad, x1 - xpad
+    y0, y1 = y0 + ypad, y1 - ypad
     xs = np.linspace(x0, x1, width)
     ys = np.linspace(y0, y1, height)
     if res[0] < 0: xs = xs[::-1]
@@ -327,10 +329,10 @@ def compute_coords(width, height, x_range, y_range, res):
 def downsample_aggregate(aggregate, factor, how='mean'):
     """Create downsampled aggregate factor in pixels units"""
     ys, xs = aggregate.shape[:2]
-    crarr = aggregate[:ys-(ys % int(factor)), :xs-(xs % int(factor))]
+    crarr = aggregate[:ys - (ys % int(factor)), :xs - (xs % int(factor))]
     concat = np.concatenate([[crarr[i::factor, j::factor]
-                            for i in range(factor)]
-                            for j in range(factor)])
+                              for i in range(factor)]
+                             for j in range(factor)])
 
     if how == 'mean':
         return np.nanmean(concat, axis=0)
@@ -381,6 +383,7 @@ def _(*args):
         if not last or last[0] != args:
             last[:] = args, f(*args)
         return last[1]
+
     return _
 
 
@@ -427,6 +430,32 @@ def lnglat_to_meters(longitude, latitude):
     return (easting, northing)
 
 
+def meters_to_lnglat(easting, northing):
+    """
+    Projects the given easting, northing values into
+    longitude, latitude coordinates.
+    easting and northing values are assumed to be in Web Mercator
+    (aka Pseudo-Mercator or EPSG:3857) coordinates.
+    Args:
+        easting
+        northing
+    Returns:
+        (longitude, latitude)
+    """
+    if isinstance(easting, (list, tuple)):
+        easting = np.array(easting)
+    if isinstance(northing, (list, tuple)):
+        northing = np.array(northing)
+
+    origin_shift = np.pi * 6378137
+    longitude = easting * 180.0 / origin_shift
+    with np.errstate(divide='ignore'):
+        latitude = np.arctan(
+            np.exp(northing * np.pi / origin_shift)
+        ) * 360.0 / np.pi - 90
+    return longitude, latitude
+
+
 # Heavily inspired by odo
 def dshape_from_pandas_helper(col):
     """Return an object from datashape.coretypes given a column from a pandas
@@ -500,33 +529,33 @@ def dshape_from_xarray_dataset(xr_ds):
 
 
 def dataframe_from_multiple_sequences(x_values, y_values):
-   """
-   Converts a set of multiple sequences (eg: time series), stored as a 2 dimensional
-   numpy array into a pandas dataframe that can be plotted by datashader.
-   The pandas dataframe eventually contains two columns ('x' and 'y') with the data.
-   Each time series is separated by a row of NaNs.
-   Discussion at: https://github.com/bokeh/datashader/issues/286#issuecomment-334619499
+    """
+    Converts a set of multiple sequences (eg: time series), stored as a 2 dimensional
+    numpy array into a pandas dataframe that can be plotted by datashader.
+    The pandas dataframe eventually contains two columns ('x' and 'y') with the data.
+    Each time series is separated by a row of NaNs.
+    Discussion at: https://github.com/bokeh/datashader/issues/286#issuecomment-334619499
 
-   x_values: 1D numpy array with the values to be plotted on the x axis (eg: time)
-   y_values: 2D numpy array with the sequences to be plotted of shape (num sequences X length of each sequence)
+    x_values: 1D numpy array with the values to be plotted on the x axis (eg: time)
+    y_values: 2D numpy array with the sequences to be plotted of shape (num sequences X length of each sequence)
 
-   """
+    """
 
-   # Add a NaN at the end of the array of x values
-   x = np.zeros(x_values.shape[0] + 1)
-   x[-1] = np.nan
-   x[:-1] = x_values
+    # Add a NaN at the end of the array of x values
+    x = np.zeros(x_values.shape[0] + 1)
+    x[-1] = np.nan
+    x[:-1] = x_values
 
-   # Tile this array of x values: number of repeats = number of sequences/time series in the data
-   x = np.tile(x, y_values.shape[0])
+    # Tile this array of x values: number of repeats = number of sequences/time series in the data
+    x = np.tile(x, y_values.shape[0])
 
-   # Add a NaN at the end of every sequence in y_values
-   y = np.zeros((y_values.shape[0], y_values.shape[1] + 1))
-   y[:, -1] = np.nan
-   y[:, :-1] = y_values
+    # Add a NaN at the end of every sequence in y_values
+    y = np.zeros((y_values.shape[0], y_values.shape[1] + 1))
+    y[:, -1] = np.nan
+    y[:, :-1] = y_values
 
-   # Return a dataframe with this new set of x and y values
-   return pd.DataFrame({'x': x, 'y': y.flatten()})
+    # Return a dataframe with this new set of x and y values
+    return pd.DataFrame({'x': x, 'y': y.flatten()})
 
 
 def _pd_mesh(vertices, simplices):
@@ -537,7 +566,7 @@ def _pd_mesh(vertices, simplices):
     winding = [0, 1, 2]
     first_tri = vertices.values[simplices.values[0, winding].astype(np.int64), :2]
     a, b, c = first_tri
-    if np.cross(b-a, c-a).item() >= 0:
+    if np.cross(b - a, c - a).item() >= 0:
         winding = [0, 2, 1]
 
     # Construct mesh by indexing into vertices with simplex indices
@@ -568,7 +597,7 @@ def _dd_mesh(vertices, simplices):
     # Compute a chunksize that will not split the vertices of a single
     # triangle across partitions
     approx_npartitions = max(vertices.npartitions, simplices.npartitions)
-    chunksize = int(np.ceil(len(res) / (3*approx_npartitions)) * 3)
+    chunksize = int(np.ceil(len(res) / (3 * approx_npartitions)) * 3)
 
     # Create dask dataframe
     res = dd.from_pandas(res, chunksize=chunksize)
@@ -592,8 +621,8 @@ def mesh(vertices, simplices):
                                 'consider casting simplices to integers '
                                 'with ".astype(int)"')
 
-    assert len(vertices.columns) > 2 or simplices.values.shape[1] > 3, 'If no vertex weight column is provided, a triangle weight column is required.'
-
+    assert len(vertices.columns) > 2 or simplices.values.shape[
+        1] > 3, 'If no vertex weight column is provided, a triangle weight column is required.'
 
     if isinstance(vertices, dd.DataFrame) and isinstance(simplices, dd.DataFrame):
         return _dd_mesh(vertices, simplices)

From e508a9a0b5ecdad9e5033a16a136396ee2db6107 Mon Sep 17 00:00:00 2001
From: Barry A Bragg Jr <hokieg3n1us@gmail.com>
Date: Fri, 15 Oct 2021 21:19:49 -0400
Subject: [PATCH 02/12] Undo Formatting

---
 datashader/utils.py | 95 ++++++++++++++++++++++-----------------------
 1 file changed, 46 insertions(+), 49 deletions(-)

diff --git a/datashader/utils.py b/datashader/utils.py
index a91097ce2..dd34f6de4 100644
--- a/datashader/utils.py
+++ b/datashader/utils.py
@@ -45,8 +45,8 @@ class VisibleDeprecationWarning(UserWarning):
 
 # Get and save the Numba version, will be used to limit functionality
 numba_version = tuple([int(x) for x in re.match(
-    r"([0-9]+)\.([0-9]+)\.([0-9]+)",
-    nb.__version__).groups()])
+                            r"([0-9]+)\.([0-9]+)\.([0-9]+)",
+                            nb.__version__).groups()])
 
 
 class Expr(object):
@@ -56,7 +56,6 @@ class Expr(object):
     ``inputs`` attribute/property, containing a tuple of everything that fully
     defines that expression.
     """
-
     def __hash__(self):
         return hash((type(self), self._hashable_inputs()))
 
@@ -86,7 +85,6 @@ def _hashable_inputs(self):
 
 class Dispatcher(object):
     """Simple single dispatch."""
-
     def __init__(self):
         self._lookup = {}
 
@@ -214,17 +212,17 @@ def calc_bbox(xs, ys, res):
 
     xmin = ymin = np.inf
     xmax = ymax = -np.inf
-    Ab = np.array([[res[0], 0., xbound],
-                   [0., -res[1], ybound],
-                   [0., 0., 1.]])
+    Ab = np.array([[res[0],  0.,      xbound],
+                   [0.,      -res[1], ybound],
+                   [0.,      0.,      1.]])
     for x_, y_ in [(0, 0), (0, len(ys)), (len(xs), 0), (len(xs), len(ys))]:
         x, y, _ = np.dot(Ab, np.array([x_, y_, 1.]))
         if x < xmin: xmin = x
         if x > xmax: xmax = x
         if y < ymin: ymin = y
         if y > ymax: ymax = y
-    xpad, ypad = res[0] / 2., res[1] / 2.
-    return xmin - xpad, ymin + ypad, xmax - xpad, ymax + ypad
+    xpad, ypad = res[0]/2., res[1]/2.
+    return xmin-xpad, ymin+ypad, xmax-xpad, ymax+ypad
 
 
 def get_indices(start, end, coords, res):
@@ -243,11 +241,11 @@ def get_indices(start, end, coords, res):
         Resolution along an axis (aka "grid/cell sizes")
     """
     size = len(coords)
-    half = abs(res) / 2.
+    half = abs(res)/2.
     vmin, vmax = coords.min(), coords.max()
-    span = vmax - vmin
-    start, end = start + half - vmin, end - half - vmin
-    sidx, eidx = int((start / span) * size), int((end / span) * size)
+    span = vmax-vmin
+    start, end = start+half-vmin, end-half-vmin
+    sidx, eidx = int((start/span)*size), int((end/span)*size)
     if eidx < sidx:
         return sidx, sidx
     return sidx, eidx
@@ -276,7 +274,7 @@ def orient_array(raster, res=None, layer=None):
     if res is None:
         res = calc_res(raster)
     array = raster.data
-    if layer is not None: array = array[layer - 1]
+    if layer is not None: array = array[layer-1]
     r0zero = np.timedelta64(0, 'ns') if isinstance(res[0], np.timedelta64) else 0
     r1zero = np.timedelta64(0, 'ns') if isinstance(res[1], np.timedelta64) else 0
     if array.ndim == 2:
@@ -314,11 +312,11 @@ def compute_coords(width, height, x_range, y_range, res):
         1D array of y-coordinates
     """
     (x0, x1), (y0, y1) = x_range, y_range
-    xd = (x1 - x0) / float(width)
-    yd = (y1 - y0) / float(height)
-    xpad, ypad = abs(xd / 2.), abs(yd / 2.)
-    x0, x1 = x0 + xpad, x1 - xpad
-    y0, y1 = y0 + ypad, y1 - ypad
+    xd = (x1-x0)/float(width)
+    yd = (y1-y0)/float(height)
+    xpad, ypad = abs(xd/2.), abs(yd/2.)
+    x0, x1 = x0+xpad, x1-xpad
+    y0, y1 = y0+ypad, y1-ypad
     xs = np.linspace(x0, x1, width)
     ys = np.linspace(y0, y1, height)
     if res[0] < 0: xs = xs[::-1]
@@ -329,10 +327,10 @@ def compute_coords(width, height, x_range, y_range, res):
 def downsample_aggregate(aggregate, factor, how='mean'):
     """Create downsampled aggregate factor in pixels units"""
     ys, xs = aggregate.shape[:2]
-    crarr = aggregate[:ys - (ys % int(factor)), :xs - (xs % int(factor))]
+    crarr = aggregate[:ys-(ys % int(factor)), :xs-(xs % int(factor))]
     concat = np.concatenate([[crarr[i::factor, j::factor]
-                              for i in range(factor)]
-                             for j in range(factor)])
+                            for i in range(factor)]
+                            for j in range(factor)])
 
     if how == 'mean':
         return np.nanmean(concat, axis=0)
@@ -383,7 +381,6 @@ def _(*args):
         if not last or last[0] != args:
             last[:] = args, f(*args)
         return last[1]
-
     return _
 
 
@@ -529,33 +526,33 @@ def dshape_from_xarray_dataset(xr_ds):
 
 
 def dataframe_from_multiple_sequences(x_values, y_values):
-    """
-    Converts a set of multiple sequences (eg: time series), stored as a 2 dimensional
-    numpy array into a pandas dataframe that can be plotted by datashader.
-    The pandas dataframe eventually contains two columns ('x' and 'y') with the data.
-    Each time series is separated by a row of NaNs.
-    Discussion at: https://github.com/bokeh/datashader/issues/286#issuecomment-334619499
+   """
+   Converts a set of multiple sequences (eg: time series), stored as a 2 dimensional
+   numpy array into a pandas dataframe that can be plotted by datashader.
+   The pandas dataframe eventually contains two columns ('x' and 'y') with the data.
+   Each time series is separated by a row of NaNs.
+   Discussion at: https://github.com/bokeh/datashader/issues/286#issuecomment-334619499
 
-    x_values: 1D numpy array with the values to be plotted on the x axis (eg: time)
-    y_values: 2D numpy array with the sequences to be plotted of shape (num sequences X length of each sequence)
+   x_values: 1D numpy array with the values to be plotted on the x axis (eg: time)
+   y_values: 2D numpy array with the sequences to be plotted of shape (num sequences X length of each sequence)
 
-    """
+   """
 
-    # Add a NaN at the end of the array of x values
-    x = np.zeros(x_values.shape[0] + 1)
-    x[-1] = np.nan
-    x[:-1] = x_values
+   # Add a NaN at the end of the array of x values
+   x = np.zeros(x_values.shape[0] + 1)
+   x[-1] = np.nan
+   x[:-1] = x_values
 
-    # Tile this array of x values: number of repeats = number of sequences/time series in the data
-    x = np.tile(x, y_values.shape[0])
+   # Tile this array of x values: number of repeats = number of sequences/time series in the data
+   x = np.tile(x, y_values.shape[0])
 
-    # Add a NaN at the end of every sequence in y_values
-    y = np.zeros((y_values.shape[0], y_values.shape[1] + 1))
-    y[:, -1] = np.nan
-    y[:, :-1] = y_values
+   # Add a NaN at the end of every sequence in y_values
+   y = np.zeros((y_values.shape[0], y_values.shape[1] + 1))
+   y[:, -1] = np.nan
+   y[:, :-1] = y_values
 
-    # Return a dataframe with this new set of x and y values
-    return pd.DataFrame({'x': x, 'y': y.flatten()})
+   # Return a dataframe with this new set of x and y values
+   return pd.DataFrame({'x': x, 'y': y.flatten()})
 
 
 def _pd_mesh(vertices, simplices):
@@ -566,7 +563,7 @@ def _pd_mesh(vertices, simplices):
     winding = [0, 1, 2]
     first_tri = vertices.values[simplices.values[0, winding].astype(np.int64), :2]
     a, b, c = first_tri
-    if np.cross(b - a, c - a).item() >= 0:
+    if np.cross(b-a, c-a).item() >= 0:
         winding = [0, 2, 1]
 
     # Construct mesh by indexing into vertices with simplex indices
@@ -597,7 +594,7 @@ def _dd_mesh(vertices, simplices):
     # Compute a chunksize that will not split the vertices of a single
     # triangle across partitions
     approx_npartitions = max(vertices.npartitions, simplices.npartitions)
-    chunksize = int(np.ceil(len(res) / (3 * approx_npartitions)) * 3)
+    chunksize = int(np.ceil(len(res) / (3*approx_npartitions)) * 3)
 
     # Create dask dataframe
     res = dd.from_pandas(res, chunksize=chunksize)
@@ -621,10 +618,10 @@ def mesh(vertices, simplices):
                                 'consider casting simplices to integers '
                                 'with ".astype(int)"')
 
-    assert len(vertices.columns) > 2 or simplices.values.shape[
-        1] > 3, 'If no vertex weight column is provided, a triangle weight column is required.'
+    assert len(vertices.columns) > 2 or simplices.values.shape[1] > 3, 'If no vertex weight column is provided, a triangle weight column is required.'
+
 
     if isinstance(vertices, dd.DataFrame) and isinstance(simplices, dd.DataFrame):
         return _dd_mesh(vertices, simplices)
 
-    return _pd_mesh(vertices, simplices)
+    return _pd_mesh(vertices, simplices)
\ No newline at end of file

From 85fb11ac2bfc36693b09a197f4a9a9bc4eb00207 Mon Sep 17 00:00:00 2001
From: Barry A Bragg Jr <hokieg3n1us@gmail.com>
Date: Fri, 15 Oct 2021 21:34:07 -0400
Subject: [PATCH 03/12] Create directories for mbtiles file.

---
 datashader/tiles.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/datashader/tiles.py b/datashader/tiles.py
index ba813237d..18666ed11 100644
--- a/datashader/tiles.py
+++ b/datashader/tiles.py
@@ -107,6 +107,8 @@ def validate_output_path(output_path, full_extent, levels):
     if os.path.isdir(output_path):
         _create_dir(output_path)
     elif output_path.endswith("mbtiles"):
+        _create_dir(os.path.dirname(output_path))
+        # Create mbtiles file and setup sqlite tables.
         MapboxTileRenderer.setup(output_path, full_extent, levels[0], levels[len(levels) - 1])
 
 

From c5d6d09fce2eadecc7afc770455c1c0977463395 Mon Sep 17 00:00:00 2001
From: Barry A Bragg Jr <hokieg3n1us@gmail.com>
Date: Fri, 15 Oct 2021 21:54:28 -0400
Subject: [PATCH 04/12] Use filename as name in mbtiles metadata

---
 datashader/tiles.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/datashader/tiles.py b/datashader/tiles.py
index 18666ed11..6d96361d9 100644
--- a/datashader/tiles.py
+++ b/datashader/tiles.py
@@ -466,7 +466,7 @@ def setup(output_location, full_extent, min_zoom, max_zoom, tile_format="PNG"):
 
         # Setup MBTiles metadata table.
         cur.execute("""insert into metadata (name, value) values (?, ?) """,
-                    ("name", "tiles"))
+                    ("name", os.path.splitext(os.path.basename(output_location))[0]))
         cur.execute("""insert into metadata (name, value) values (?, ?) """,
                     ("format", tile_format.lower()))
         cur.execute("""insert into metadata (name, value) values (?, ?) """,

From 64d28dca168dfc70511219481b2c3e2f4571dca0 Mon Sep 17 00:00:00 2001
From: Barry A Bragg Jr <hokieg3n1us@gmail.com>
Date: Mon, 18 Oct 2021 19:06:02 -0400
Subject: [PATCH 05/12] Minor updates

Compute min/max of span in-place instead of building a list of all values, then computing using dask. Use dimensions of data array for label based indexing, instead of hard-coded x, y labels (which forced input coordinate columns to be named x, y).
---
 datashader/tiles.py | 20 ++++++++++++++------
 1 file changed, 14 insertions(+), 6 deletions(-)

diff --git a/datashader/tiles.py b/datashader/tiles.py
index 6d96361d9..fb1027a8d 100644
--- a/datashader/tiles.py
+++ b/datashader/tiles.py
@@ -39,7 +39,8 @@ def calculate_zoom_level_stats(super_tiles, load_data_func,
                                rasterize_func,
                                color_ranging_strategy='fullscan'):
     if color_ranging_strategy == 'fullscan':
-        stats = []
+        span_min = None
+        span_max = None
         is_bool = False
         for super_tile in super_tiles:
             agg = _get_super_tile_min_max(super_tile, load_data_func, rasterize_func)
@@ -47,13 +48,19 @@ def calculate_zoom_level_stats(super_tiles, load_data_func,
             if agg.dtype.kind == 'b':
                 is_bool = True
             else:
-                stats.append(np.nanmin(agg.data))
-                stats.append(np.nanmax(agg.data))
+                if span_min is None:
+                    span_min = np.nanmin(agg.data)
+                else:
+                    span_min = min(span_min, np.nanmin(agg.data))
+
+                if span_max is None:
+                    span_max = np.nanmax(agg.data)
+                else:
+                    span_max = max(span_max, np.nanmax(agg.data))
         if is_bool:
             span = (0, 1)
         else:
-            b = db.from_sequence(stats)
-            span = dask.compute(b.min(), b.max())
+            span = (span_min, span_max)
         return super_tiles, span
     else:
         raise ValueError('Invalid color_ranging_strategy option')
@@ -314,7 +321,8 @@ def render(self, da, level):
         for t in tiles:
             x, y, z, data_extent = t
             dxmin, dymin, dxmax, dymax = data_extent
-            arr = da.loc[{'x': slice(dxmin, dxmax), 'y': slice(dymin, dymax)}]
+
+            arr = da.loc[{da.dims[1]: slice(dxmin, dxmax), da.dims[0]: slice(dymin, dymax)}]
 
             if 0 in arr.shape:
                 continue

From 3b1c9d4afdb7aa30e2d1cfa1177a44a1ffec5097 Mon Sep 17 00:00:00 2001
From: Barry A Bragg Jr <hokieg3n1us@gmail.com>
Date: Wed, 20 Oct 2021 07:43:00 -0400
Subject: [PATCH 06/12] Out-of-Core Processing

Include capability to provide a local_cache_path. If provided, the aggregation for the super tiles will be stored locally as a NetCDF file (instead of keeping all these in memory, which eventually overflows memory at high zoom levels). This allows most of the processing to now be done completely out of core.
---
 datashader/tiles.py | 59 ++++++++++++++++++++++++++++++---------------
 1 file changed, 40 insertions(+), 19 deletions(-)

diff --git a/datashader/tiles.py b/datashader/tiles.py
index fb1027a8d..afef7ef9a 100644
--- a/datashader/tiles.py
+++ b/datashader/tiles.py
@@ -5,9 +5,9 @@
 import sqlite3
 from io import BytesIO
 
-import dask
 import dask.bag as db
 import numpy as np
+import xarray
 from PIL.Image import fromarray
 
 from .utils import meters_to_lnglat
@@ -37,14 +37,22 @@ def _get_super_tile_min_max(tile_info, load_data_func, rasterize_func):
 
 def calculate_zoom_level_stats(super_tiles, load_data_func,
                                rasterize_func,
-                               color_ranging_strategy='fullscan'):
+                               color_ranging_strategy='fullscan', local_cache_path=None):
     if color_ranging_strategy == 'fullscan':
         span_min = None
         span_max = None
         is_bool = False
+        index = 0
         for super_tile in super_tiles:
             agg = _get_super_tile_min_max(super_tile, load_data_func, rasterize_func)
-            super_tile['agg'] = agg
+
+            if local_cache_path:
+                cache_file_path = os.path.join(local_cache_path, 'super_tile_' + str(index) + '.nc')
+                agg.to_netcdf(cache_file_path)
+                super_tile['agg'] = cache_file_path
+            else:
+                super_tile['agg'] = agg
+
             if agg.dtype.kind == 'b':
                 is_bool = True
             else:
@@ -57,6 +65,8 @@ def calculate_zoom_level_stats(super_tiles, load_data_func,
                     span_max = np.nanmax(agg.data)
                 else:
                     span_max = max(span_max, np.nanmax(agg.data))
+            index = index + 1
+
         if is_bool:
             span = (0, 1)
         else:
@@ -68,19 +78,22 @@ def calculate_zoom_level_stats(super_tiles, load_data_func,
 
 def render_tiles(full_extent, levels, load_data_func,
                  rasterize_func, shader_func,
-                 post_render_func, output_path, color_ranging_strategy='fullscan', num_workers=4):
+                 post_render_func, output_path, color_ranging_strategy='fullscan', num_workers=4,
+                 local_cache_path=None):
     results = dict()
 
-    validate_output_path(output_path, full_extent, levels)
+    validate_output_path(output_path, full_extent, levels, local_cache_path)
 
     for level in levels:
         print('calculating statistics for level {}'.format(level))
         super_tiles, span = calculate_zoom_level_stats(list(gen_super_tiles(full_extent, level)),
                                                        load_data_func, rasterize_func,
-                                                       color_ranging_strategy=color_ranging_strategy)
+                                                       color_ranging_strategy=color_ranging_strategy,
+                                                       local_cache_path=local_cache_path)
         print('rendering {} supertiles for zoom level {} with span={}'.format(len(super_tiles), level, span))
         b = db.from_sequence(super_tiles)
-        b.map(render_super_tile, span, output_path, shader_func, post_render_func).compute(num_workers=num_workers)
+        b.map(render_super_tile, span, output_path, shader_func, post_render_func, local_cache_path).compute(
+            num_workers=num_workers)
         results[level] = dict(success=True, stats=span, supertile_count=len(super_tiles))
 
     return results
@@ -91,8 +104,8 @@ def gen_super_tiles(extent, zoom_level, span=None):
     super_tile_size = min(2 ** 4 * 256,
                           (2 ** zoom_level) * 256)
     super_tile_def = MercatorTileDefinition(x_range=(xmin, xmax), y_range=(ymin, ymax), tile_size=super_tile_size)
-    super_tiles = super_tile_def.get_tiles_by_extent(extent, zoom_level)
-    for s in super_tiles:
+
+    for s in super_tile_def.get_tiles_by_extent(extent, zoom_level):
         st_extent = s[3]
         x_range = (st_extent[0], st_extent[2])
         y_range = (st_extent[1], st_extent[3])
@@ -103,13 +116,23 @@ def gen_super_tiles(extent, zoom_level, span=None):
                'span': span}
 
 
-def render_super_tile(tile_info, span, output_path, shader_func, post_render_func):
+def render_super_tile(tile_info, span, output_path, shader_func, post_render_func, local_cache_path):
     level = tile_info['level']
-    ds_img = shader_func(tile_info['agg'], span=span)
+
+    agg = None
+
+    if local_cache_path is not None:
+        agg = xarray.open_dataarray(tile_info['agg']).astype('uint32')
+        os.remove(tile_info['agg'])
+    else:
+        agg = tile_info['agg']
+
+    ds_img = shader_func(agg, span=span)
+
     return create_sub_tiles(ds_img, level, tile_info, output_path, post_render_func)
 
 
-def validate_output_path(output_path, full_extent, levels):
+def validate_output_path(output_path, full_extent, levels, local_cache_path):
     # validate / createoutput_dir
     if os.path.isdir(output_path):
         _create_dir(output_path)
@@ -118,6 +141,9 @@ def validate_output_path(output_path, full_extent, levels):
         # Create mbtiles file and setup sqlite tables.
         MapboxTileRenderer.setup(output_path, full_extent, levels[0], levels[len(levels) - 1])
 
+    if local_cache_path:
+        _create_dir(local_cache_path)
+
 
 def create_sub_tiles(data_array, level, tile_info, output_path, post_render_func=None):
     # create tile source
@@ -281,15 +307,11 @@ def get_tiles_by_extent(self, extent, level):
         txmin, tymax = self.meters_to_tile(xmin, ymin, level)
         txmax, tymin = self.meters_to_tile(xmax, ymax, level)
 
-        # TODO: vectorize?
-        tiles = []
         for ty in range(tymin, tymax + 1):
             for tx in range(txmin, txmax + 1):
                 if self.is_valid_tile(tx, ty, level):
                     t = (tx, ty, level, self.get_tile_meters(tx, ty, level))
-                    tiles.append(t)
-
-        return tiles
+                    yield t
 
     def get_tile_meters(self, tx, ty, level):
         ty = invert_y_tile(ty, level)  # convert to TMS for conversion to meters
@@ -317,8 +339,7 @@ def render(self, da, level):
         ymin, ymax = self.tile_def.y_range
         extent = xmin, ymin, xmax, ymax
 
-        tiles = self.tile_def.get_tiles_by_extent(extent, level)
-        for t in tiles:
+        for t in self.tile_def.get_tiles_by_extent(extent, level):
             x, y, z, data_extent = t
             dxmin, dymin, dxmax, dymax = data_extent
 

From 9ff9ca647c46224d48aa955c46074550f2d09602 Mon Sep 17 00:00:00 2001
From: Barry A Bragg Jr <hokieg3n1us@gmail.com>
Date: Wed, 20 Oct 2021 20:41:18 -0400
Subject: [PATCH 07/12] Use netCDF4 for caching

Use the netCDF4 library for writing/loading XArray DataArrays from cache, since it properly handles unsigned data types. Optional dependency that will raise an error if not installed when the caching feature is used.
---
 datashader/tiles.py | 10 ++++++++--
 1 file changed, 8 insertions(+), 2 deletions(-)

diff --git a/datashader/tiles.py b/datashader/tiles.py
index afef7ef9a..3f7b0a8d1 100644
--- a/datashader/tiles.py
+++ b/datashader/tiles.py
@@ -12,6 +12,11 @@
 
 from .utils import meters_to_lnglat
 
+try:
+    import netCDF4
+except Exception:
+    netCDF4 = None
+
 __all__ = ['render_tiles', 'MercatorTileDefinition']
 
 
@@ -48,7 +53,7 @@ def calculate_zoom_level_stats(super_tiles, load_data_func,
 
             if local_cache_path:
                 cache_file_path = os.path.join(local_cache_path, 'super_tile_' + str(index) + '.nc')
-                agg.to_netcdf(cache_file_path)
+                agg.to_netcdf(cache_file_path, engine='netcdf4', format='NETCDF4')
                 super_tile['agg'] = cache_file_path
             else:
                 super_tile['agg'] = agg
@@ -122,7 +127,7 @@ def render_super_tile(tile_info, span, output_path, shader_func, post_render_fun
     agg = None
 
     if local_cache_path is not None:
-        agg = xarray.open_dataarray(tile_info['agg']).astype('uint32')
+        agg = xarray.load_dataarray(tile_info['agg'], engine='netcdf4')
         os.remove(tile_info['agg'])
     else:
         agg = tile_info['agg']
@@ -142,6 +147,7 @@ def validate_output_path(output_path, full_extent, levels, local_cache_path):
         MapboxTileRenderer.setup(output_path, full_extent, levels[0], levels[len(levels) - 1])
 
     if local_cache_path:
+        assert netCDF4, 'netcdf4 library must be installed for use with local_cache.'
         _create_dir(local_cache_path)
 
 

From 6d0173b49cdfde3eec7ecbe38ccfdfc2c640f45b Mon Sep 17 00:00:00 2001
From: Barry A Bragg Jr <hokieg3n1us@gmail.com>
Date: Mon, 25 Oct 2021 18:55:25 -0400
Subject: [PATCH 08/12] Check if path is a directory, before calling makedirs

---
 datashader/tiles.py | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/datashader/tiles.py b/datashader/tiles.py
index 3f7b0a8d1..655035747 100644
--- a/datashader/tiles.py
+++ b/datashader/tiles.py
@@ -25,7 +25,8 @@ def _create_dir(path):
     import os, errno
 
     try:
-        os.makedirs(path)
+        if os.path.isdir(path):
+            os.makedirs(path)
     except OSError as e:
         if e.errno != errno.EEXIST:
             raise
@@ -139,9 +140,9 @@ def render_super_tile(tile_info, span, output_path, shader_func, post_render_fun
 
 def validate_output_path(output_path, full_extent, levels, local_cache_path):
     # validate / createoutput_dir
-    if os.path.isdir(output_path):
-        _create_dir(output_path)
-    elif output_path.endswith("mbtiles"):
+    _create_dir(output_path)
+
+    if output_path.endswith("mbtiles"):
         _create_dir(os.path.dirname(output_path))
         # Create mbtiles file and setup sqlite tables.
         MapboxTileRenderer.setup(output_path, full_extent, levels[0], levels[len(levels) - 1])

From b6f3e80ebe9d3b05ad3b57e4be908854ac220ec6 Mon Sep 17 00:00:00 2001
From: Barry A Bragg Jr <hokieg3n1us@gmail.com>
Date: Tue, 26 Oct 2021 16:24:24 -0400
Subject: [PATCH 09/12] SQL consistency.

Allow updates to an existing mbtiles file, inserting or replacing rows in the SQLite database and make SQL statements consistent.
---
 datashader/tiles.py | 34 +++++++++++++++++-----------------
 1 file changed, 17 insertions(+), 17 deletions(-)

diff --git a/datashader/tiles.py b/datashader/tiles.py
index 655035747..9c6d1743e 100644
--- a/datashader/tiles.py
+++ b/datashader/tiles.py
@@ -88,7 +88,7 @@ def render_tiles(full_extent, levels, load_data_func,
                  local_cache_path=None):
     results = dict()
 
-    validate_output_path(output_path, full_extent, levels, local_cache_path)
+    _setup(full_extent, levels, output_path, local_cache_path)
 
     for level in levels:
         print('calculating statistics for level {}'.format(level))
@@ -138,7 +138,7 @@ def render_super_tile(tile_info, span, output_path, shader_func, post_render_fun
     return create_sub_tiles(ds_img, level, tile_info, output_path, post_render_func)
 
 
-def validate_output_path(output_path, full_extent, levels, local_cache_path):
+def _setup(full_extent, levels, output_path, local_cache_path):
     # validate / createoutput_dir
     _create_dir(output_path)
 
@@ -468,16 +468,16 @@ def setup(output_location, full_extent, min_zoom, max_zoom, tile_format="PNG"):
 
         # Create MBTiles tables.
         cur.execute("""
-            create table tiles (
+            create table if not exists tiles (
                 zoom_level integer,
                 tile_column integer,
                 tile_row integer,
                 tile_data blob);
                 """)
-        cur.execute("""create table metadata
+        cur.execute("""create table if not exists metadata
             (name text, value text);""")
-        cur.execute("""create unique index name on metadata (name);""")
-        cur.execute("""create unique index tile_index on tiles
+        cur.execute("""create unique index if not exists name on metadata (name);""")
+        cur.execute("""create unique index if not exists tile_index on tiles
             (zoom_level, tile_column, tile_row);""")
 
         # Compute Extents in WGS84
@@ -501,18 +501,18 @@ def setup(output_location, full_extent, min_zoom, max_zoom, tile_format="PNG"):
                 zoom_level = 1
 
         # Setup MBTiles metadata table.
-        cur.execute("""insert into metadata (name, value) values (?, ?) """,
+        cur.execute("""insert or replace into metadata (name, value) values (?, ?);""",
                     ("name", os.path.splitext(os.path.basename(output_location))[0]))
-        cur.execute("""insert into metadata (name, value) values (?, ?) """,
+        cur.execute("""insert or replace into metadata (name, value) values (?, ?);""",
                     ("format", tile_format.lower()))
-        cur.execute("""insert into metadata (name, value) values (?, ?) """,
-                    ("bounds ", "{},{},{},{}".format(min_lon, min_lat, max_lon, max_lat)))
-        cur.execute("""insert into metadata (name, value) values (?, ?) """,
-                    ("center ", "{},{},{}".format(lon_center, lat_center, zoom_level)))
-        cur.execute("""insert into metadata (name, value) values (?, ?) """,
-                    ("minzoom ", min_zoom))
-        cur.execute("""insert into metadata (name, value) values (?, ?) """,
-                    ("maxzoom ", max_zoom))
+        cur.execute("""insert or replace into metadata (name, value) values (?, ?);""",
+                    ("bounds", "{},{},{},{}".format(min_lon, min_lat, max_lon, max_lat)))
+        cur.execute("""insert or replace into metadata (name, value) values (?, ?);""",
+                    ("center", "{},{},{}".format(lon_center, lat_center, zoom_level)))
+        cur.execute("""insert or replace into metadata (name, value) values (?, ?);""",
+                    ("minzoom", min_zoom))
+        cur.execute("""insert or replace into metadata (name, value) values (?, ?);""",
+                    ("maxzoom", max_zoom))
 
         cur.close()
         con.commit()
@@ -529,7 +529,7 @@ def render(self, da, level):
 
             tile_row = (2 ** z) - 1 - y
 
-            cur.execute("""insert into tiles (zoom_level,
+            cur.execute("""insert or replace into tiles (zoom_level,
                 tile_column, tile_row, tile_data) values
                 (?, ?, ?, ?);""",
                         (z, x, tile_row, sqlite3.Binary(image_bytes.getvalue())))

From 9435272c09eb9781180a7daf00719d08710e0e57 Mon Sep 17 00:00:00 2001
From: Barry A Bragg Jr <hokieg3n1us@gmail.com>
Date: Fri, 21 Jan 2022 00:03:17 -0500
Subject: [PATCH 10/12] Memory Overhead Tradeoff

calculate_zoom_level_stats moved to a parallel implementation with dask bags, that does not cache the super tiles in memory. If a local_cache_path is provided, the super tiles will be cached to disk as NetCDF files. Super tiles will then be either recomputed during the render_super_tiles function, or loaded from the local cache.

Note: If using the Datashader render_tiles, the scheduler for Dask should be configured to 'threads' (when running on local) to avoid Dask bag computations from copying the input DataFrame to multiple processes, which will increase memory overhead. If outputting tiles to the MBTiles format, the num_workers should be tuned to prevent the SQLite database from locking during transactions. Both of these can be configured using dask.config.set(scheduler='threads', num_workers=4).
---
 datashader/tests/test_tiles.py | 26 +++++++------
 datashader/tiles.py            | 69 ++++++++++++++++++----------------
 2 files changed, 51 insertions(+), 44 deletions(-)

diff --git a/datashader/tests/test_tiles.py b/datashader/tests/test_tiles.py
index a1dd2af71..2e240ce78 100644
--- a/datashader/tests/test_tiles.py
+++ b/datashader/tests/test_tiles.py
@@ -18,6 +18,8 @@
 MERCATOR_CONST = 20037508.34
 
 df = None
+
+
 def mock_load_data_func(x_range, y_range):
     global df
     if df is None:
@@ -89,39 +91,39 @@ def assert_is_numeric(value):
     assert any([is_int_or_float, is_numpy_int_or_float])
 
 
-
 def test_get_super_tile_min_max():
-
     tile_info = {'level': 0,
-                'x_range': (-MERCATOR_CONST, MERCATOR_CONST),
-                'y_range': (-MERCATOR_CONST, MERCATOR_CONST),
-                'tile_size': 256,
-                'span': (0, 1000)}
+                 'x_range': (-MERCATOR_CONST, MERCATOR_CONST),
+                 'y_range': (-MERCATOR_CONST, MERCATOR_CONST),
+                 'tile_size': 256,
+                 'span': (0, 1000)}
 
-    agg = _get_super_tile_min_max(tile_info, mock_load_data_func, mock_rasterize_func)
+    tile_info = _get_super_tile_min_max(tile_info, mock_load_data_func, mock_rasterize_func)
 
-    result = [np.nanmin(agg.data), np.nanmax(agg.data)]
+    result = [tile_info['span_min'], tile_info['span_max']]
 
     assert isinstance(result, list)
     assert len(result) == 2
     assert_is_numeric(result[0])
     assert_is_numeric(result[1])
 
+
 def test_calculate_zoom_level_stats_with_fullscan_ranging_strategy():
     full_extent = (-MERCATOR_CONST, -MERCATOR_CONST,
                    MERCATOR_CONST, MERCATOR_CONST)
     level = 0
     color_ranging_strategy = 'fullscan'
     super_tiles, span = calculate_zoom_level_stats(list(gen_super_tiles(full_extent, level)),
-                                        mock_load_data_func,
-                                        mock_rasterize_func,
-                                        color_ranging_strategy=color_ranging_strategy)
+                                                   mock_load_data_func,
+                                                   mock_rasterize_func,
+                                                   color_ranging_strategy=color_ranging_strategy)
 
     assert isinstance(span, (list, tuple))
     assert len(span) == 2
     assert_is_numeric(span[0])
     assert_is_numeric(span[1])
 
+
 def test_meters_to_tile():
     # Part of NYC (used in taxi demo)
     full_extent_of_data = (-8243206.93436, 4968192.04221, -8226510.539480001, 4982886.20438)
@@ -129,4 +131,4 @@ def test_meters_to_tile():
     zoom = 12
     tile_def = MercatorTileDefinition((xmin, xmax), (ymin, ymax), tile_size=256)
     tile = tile_def.meters_to_tile(xmin, ymin, zoom)
-    assert tile == (1205, 1540) # using Google tile coordinates, not TMS
+    assert tile == (1205, 1540)  # using Google tile coordinates, not TMS
diff --git a/datashader/tiles.py b/datashader/tiles.py
index 9c6d1743e..8940407c4 100644
--- a/datashader/tiles.py
+++ b/datashader/tiles.py
@@ -3,6 +3,7 @@
 import math
 import os
 import sqlite3
+import uuid
 from io import BytesIO
 
 import dask.bag as db
@@ -32,46 +33,48 @@ def _create_dir(path):
             raise
 
 
-def _get_super_tile_min_max(tile_info, load_data_func, rasterize_func):
-    tile_size = tile_info['tile_size']
+def _get_super_tile_min_max(tile_info, load_data_func, rasterize_func, local_cache_path=None):
     df = load_data_func(tile_info['x_range'], tile_info['y_range'])
     agg = rasterize_func(df, x_range=tile_info['x_range'],
                          y_range=tile_info['y_range'],
-                         height=tile_size, width=tile_size)
-    return agg
+                         height=tile_info['tile_size'], width=tile_info['tile_size'])
+    tile_info['agg_type'] = agg.dtype.kind
+    tile_info['span_min'] = np.nanmin(agg.data)
+    tile_info['span_max'] = np.nanmax(agg.data)
+
+    if local_cache_path:
+        cache_file_path = os.path.join(local_cache_path, 'super_tile_' + str(uuid.uuid4()) + '.nc')
+        agg.to_netcdf(cache_file_path, engine='netcdf4', format='NETCDF4')
+        tile_info['cache_file'] = cache_file_path
+
+    return tile_info
 
 
 def calculate_zoom_level_stats(super_tiles, load_data_func,
                                rasterize_func,
                                color_ranging_strategy='fullscan', local_cache_path=None):
     if color_ranging_strategy == 'fullscan':
+
+        super_tiles = db.from_sequence(super_tiles).map(_get_super_tile_min_max, load_data_func,
+                                                        rasterize_func, local_cache_path).compute()
+
         span_min = None
         span_max = None
         is_bool = False
-        index = 0
-        for super_tile in super_tiles:
-            agg = _get_super_tile_min_max(super_tile, load_data_func, rasterize_func)
-
-            if local_cache_path:
-                cache_file_path = os.path.join(local_cache_path, 'super_tile_' + str(index) + '.nc')
-                agg.to_netcdf(cache_file_path, engine='netcdf4', format='NETCDF4')
-                super_tile['agg'] = cache_file_path
-            else:
-                super_tile['agg'] = agg
 
-            if agg.dtype.kind == 'b':
+        for super_tile in super_tiles:
+            if super_tile['agg_type'] == 'b':
                 is_bool = True
             else:
                 if span_min is None:
-                    span_min = np.nanmin(agg.data)
+                    span_min = super_tile['span_min']
                 else:
-                    span_min = min(span_min, np.nanmin(agg.data))
+                    span_min = min(span_min, super_tile['span_min'])
 
                 if span_max is None:
-                    span_max = np.nanmax(agg.data)
+                    span_max = super_tile['span_max']
                 else:
-                    span_max = max(span_max, np.nanmax(agg.data))
-            index = index + 1
+                    span_max = max(span_max, super_tile['span_max'])
 
         if is_bool:
             span = (0, 1)
@@ -84,7 +87,7 @@ def calculate_zoom_level_stats(super_tiles, load_data_func,
 
 def render_tiles(full_extent, levels, load_data_func,
                  rasterize_func, shader_func,
-                 post_render_func, output_path, color_ranging_strategy='fullscan', num_workers=4,
+                 post_render_func, output_path, color_ranging_strategy='fullscan',
                  local_cache_path=None):
     results = dict()
 
@@ -98,8 +101,8 @@ def render_tiles(full_extent, levels, load_data_func,
                                                        local_cache_path=local_cache_path)
         print('rendering {} supertiles for zoom level {} with span={}'.format(len(super_tiles), level, span))
         b = db.from_sequence(super_tiles)
-        b.map(render_super_tile, span, output_path, shader_func, post_render_func, local_cache_path).compute(
-            num_workers=num_workers)
+        b.map(render_super_tile, span, output_path, load_data_func, rasterize_func, shader_func, post_render_func,
+              local_cache_path).compute()
         results[level] = dict(success=True, stats=span, supertile_count=len(super_tiles))
 
     return results
@@ -122,20 +125,22 @@ def gen_super_tiles(extent, zoom_level, span=None):
                'span': span}
 
 
-def render_super_tile(tile_info, span, output_path, shader_func, post_render_func, local_cache_path):
-    level = tile_info['level']
-
+def render_super_tile(tile_info, span, output_path, load_data_func, rasterize_func, shader_func, post_render_func,
+                      local_cache_path):
     agg = None
 
     if local_cache_path is not None:
-        agg = xarray.load_dataarray(tile_info['agg'], engine='netcdf4')
-        os.remove(tile_info['agg'])
+        agg = xarray.load_dataarray(tile_info['cache_file'], engine='netcdf4')
+        os.remove(tile_info['cache_file'])
     else:
-        agg = tile_info['agg']
+        df = load_data_func(tile_info['x_range'], tile_info['y_range'])
+        agg = rasterize_func(df, x_range=tile_info['x_range'],
+                             y_range=tile_info['y_range'],
+                             height=tile_info['tile_size'], width=tile_info['tile_size'])
 
     ds_img = shader_func(agg, span=span)
 
-    return create_sub_tiles(ds_img, level, tile_info, output_path, post_render_func)
+    return create_sub_tiles(ds_img, tile_info, output_path, post_render_func)
 
 
 def _setup(full_extent, levels, output_path, local_cache_path):
@@ -152,7 +157,7 @@ def _setup(full_extent, levels, output_path, local_cache_path):
         _create_dir(local_cache_path)
 
 
-def create_sub_tiles(data_array, level, tile_info, output_path, post_render_func=None):
+def create_sub_tiles(data_array, tile_info, output_path, post_render_func=None):
     # create tile source
     tile_def = MercatorTileDefinition(x_range=tile_info['x_range'],
                                       y_range=tile_info['y_range'],
@@ -169,7 +174,7 @@ def create_sub_tiles(data_array, level, tile_info, output_path, post_render_func
         renderer = FileSystemTileRenderer(tile_def, output_location=output_path,
                                           post_render_func=post_render_func)
 
-    return renderer.render(data_array, level=level)
+    return renderer.render(data_array, level=tile_info['level'])
 
 
 def invert_y_tile(y, z):

From 11bc3d6a550133846ed71a6d7689608232a6dad1 Mon Sep 17 00:00:00 2001
From: Barry A Bragg Jr <hokieg3n1us@gmail.com>
Date: Fri, 21 Jan 2022 00:38:45 -0500
Subject: [PATCH 11/12] Path handling for MBTiles output.

Handle setup of output paths.
---
 datashader/tiles.py | 15 +++++++++------
 1 file changed, 9 insertions(+), 6 deletions(-)

diff --git a/datashader/tiles.py b/datashader/tiles.py
index 8940407c4..038c83c4b 100644
--- a/datashader/tiles.py
+++ b/datashader/tiles.py
@@ -26,8 +26,7 @@ def _create_dir(path):
     import os, errno
 
     try:
-        if os.path.isdir(path):
-            os.makedirs(path)
+        os.makedirs(path)
     except OSError as e:
         if e.errno != errno.EEXIST:
             raise
@@ -144,13 +143,17 @@ def render_super_tile(tile_info, span, output_path, load_data_func, rasterize_fu
 
 
 def _setup(full_extent, levels, output_path, local_cache_path):
-    # validate / createoutput_dir
-    _create_dir(output_path)
+    if os.path.splitext(output_path)[1] == '.mbtiles':
+
+        output_dir = os.path.dirname(output_path)
+
+        if output_dir:
+            _create_dir(output_dir)
 
-    if output_path.endswith("mbtiles"):
-        _create_dir(os.path.dirname(output_path))
         # Create mbtiles file and setup sqlite tables.
         MapboxTileRenderer.setup(output_path, full_extent, levels[0], levels[len(levels) - 1])
+    else:
+        _create_dir(output_path)
 
     if local_cache_path:
         assert netCDF4, 'netcdf4 library must be installed for use with local_cache.'

From a0e6f1cbb0fdb66f99f23fa5b7f5b324fd0df7f8 Mon Sep 17 00:00:00 2001
From: Barry A Bragg Jr <36244490+hokieg3n1us@users.noreply.github.com>
Date: Thu, 18 May 2023 21:00:23 -0400
Subject: [PATCH 12/12] Update tiles.py

---
 datashader/tiles.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/datashader/tiles.py b/datashader/tiles.py
index 0eecd777c..8e122707a 100644
--- a/datashader/tiles.py
+++ b/datashader/tiles.py
@@ -1,12 +1,10 @@
 from __future__ import annotations
 from io import BytesIO
 
-
 import math
 import os
 import sqlite3
 import uuid
-from io import BytesIO
 
 import dask.bag as db
 import numpy as np