OSM tiling

examples/working_with_tiled_raster_data.ipynb

@@ -151,12 +151,10 @@
     "from matplotlib import pyplot as plt\n",
     "\n",
     "\n",
-    "name_png = f\"{source}_png_cmap\"\n",
+    "name_png = f\"{source}_png\"\n",
     "root_png = join(\"tmpdir\", name_png)\n",
     "da0.raster.to_slippy_tiles(\n",
     "    root=root_png,\n",
-    "    cmap=\"terrain\",\n",
-    "    norm=plt.Normalize(vmin=0, vmax=2000),\n",
     ")"
    ]
   },
@@ -199,7 +197,59 @@
     "\n",
     "If one were to see what this would look like in e.g. QGIS, a local sever is needed. \n",
     "With python this is easily done with the command `python -m http.server 8000` from the command line while within the folder where the tiles are located. In this case that would be 'root_png'.\n",
-    "In QGIS, make a new XYZ Tiles connection. For this new connection the URL becomes 'http://localhost:8000/{z}/{x}/{y}.png' and the interpolation is set to Terrarium Terrain RGB.\n"
+    "In QGIS, make a new XYZ Tiles connection. For this new connection the URL becomes 'http://localhost:8000/{z}/{x}/{y}.png' and the interpolation is set to Terrarium Terrain RGB.\n",
+    "\n",
+    "However, if the images are meant to be viewed as is, then a custom colormap can be defined to make them look nice!\n",
+    "\n",
+    "Let's make another dataset of png's!\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from matplotlib import pyplot as plt\n",
+    "\n",
+    "\n",
+    "name_png_cmap = f\"{source}_png_cmap\"\n",
+    "root_png_cmap = join(\"tmpdir\", name_png_cmap)\n",
+    "# let us define a nice color for a terrain image\n",
+    "da0.raster.to_slippy_tiles(\n",
+    "    root=root_png_cmap,\n",
+    "    cmap=\"terrain\",\n",
+    "    norm=plt.Normalize(vmin=0, vmax=2000),\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now lets take a look at the improved standard visuals!"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# View the data\n",
+    "from PIL import Image\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "# Create a figure to show the image\n",
+    "fig = plt.figure()\n",
+    "fig.set_size_inches(2.5, 2.5)\n",
+    "ax = fig.add_subplot(111)\n",
+    "ax.set_position([0, 0, 1, 1])\n",
+    "ax.axis(\"off\")\n",
+    "\n",
+    "# Show the image\n",
+    "im = Image.open(join(root_png_cmap, \"9\", \"273\", \"182.png\"))\n",
+    "ax.imshow(im)"
    ]
   },
   {
@@ -258,7 +308,7 @@
    "outputs": [],
    "source": [
     "# And for OSM\n",
-    "da1 = data_catalog.get_rasterdataset(name_osm)\n",
+    "da1 = data_catalog.get_rasterdataset(name_osm, zoom_level=11)\n",
     "da1.raster.shape"
    ]
   },
@@ -355,7 +405,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.4"
+   "version": "3.11.5"
   },
   "orig_nbformat": 4,
   "vscode": {

hydromt/raster.py

@@ -2271,8 +2271,8 @@ def to_slippy_tiles(
         min_lvl: int = None,
         max_lvl: int = None,
         driver="png",
-        cmap=None,
-        norm=None,
+        cmap: str | object = None,
+        norm: object = None,
         **kwargs,
     ):
         """Produce tiles in /zoom/x/y.<ext> structure (EPSG:3857).
@@ -2292,6 +2292,13 @@ def to_slippy_tiles(
             If None, the zoomlevels will be determined based on the data resolution
         driver : str, optional
             file output driver, one of 'png', 'netcdf4' or 'GTiff'
+        cmap : str | object, optional
+            A colormap, either defined by a string and imported from matplotlib
+            via that string or as a ListedColormap object from matplotlib itself.
+        norm : object, optional
+            A matplotlib Normalize object that defines a range between a maximum
+            and minimum value
+
         **kwargs
             Key-word arguments to write file
             for netcdf4, these are passed to ~:py:meth:xarray.DataArray.to_netcdf:
@@ -2309,6 +2316,10 @@ def to_slippy_tiles(
 
         # Fixed pixel size and CRS for XYZ tiles
         pxs = 256
+        # Extent in y-direction for pseudo mercator (EPSG:3857)
+        y_ext = math.atan(math.sinh(math.pi)) * (180 / math.pi)
+        y_ext_pm = mct.xy(0, y_ext)[1]
+
         crs = CRS.from_epsg(3857)
         ext = {"png": "png", "netcdf4": "nc", "gtiff": "tif"}.get(driver.lower(), None)
         if ext is None:
@@ -2322,6 +2333,8 @@ def to_slippy_tiles(
         # make sure dataarray has a name
         name = obj.name or "data"
         obj.name = name
+        obj_res = obj.raster.res[0]
+        obj_bounds = list(obj.raster.bounds)
 
         # colormap output
         if cmap is not None and driver != "png":
@@ -2340,28 +2353,61 @@ def to_slippy_tiles(
         }.get(driver.lower(), {})
         kwargs = {**kwargs0, **kwargs}
 
+        # Setting up information for zoomlevel calculation and
+        # determination of tile windows
+        # This section is purely for the resolution
+        bounds_4326_clip = list(obj.raster.transform_bounds("EPSG:4326"))
+        bounds_4326_clip[:2] = map(
+            max,
+            zip(bounds_4326_clip[:2], (-180, -y_ext)),
+        )
+        bounds_4326_clip[2:] = map(
+            min,
+            zip(bounds_4326_clip[2:], (180, y_ext)),
+        )
+        obj_clipped_to_pseudo = obj.raster.clip_bbox(
+            bounds_4326_clip,
+            crs="EPSG:4326",
+        )
+        tr_3857 = rasterio.warp.calculate_default_transform(
+            obj_clipped_to_pseudo.raster.crs,
+            "EPSG:3857",
+            *obj_clipped_to_pseudo.shape,
+            *obj_clipped_to_pseudo.raster.bounds,
+        )[0]
+
+        del obj_clipped_to_pseudo
+
+        # This section is for dealing with rounding errors
+        obj_bounds_cor = [
+            obj_bounds[0] + 0.5 * obj_res,
+            obj_bounds[1] + 0.5 * obj_res,
+            obj_bounds[2] - 0.5 * obj_res,
+            obj_bounds[3] - 0.5 * obj_res,
+        ]
+        bounds_4326 = rasterio.warp.transform_bounds(
+            obj.raster.crs,
+            "EPSG:4326",
+            *obj_bounds_cor,
+        )
+
         # Calculate min/max zoomlevel based
-        bounds_wgsg84 = obj.raster.transform_bounds("EPSG:4326")
         if max_lvl is None:  # calculate max zoomlevel close to native resolution
-            max_lat = max(np.abs(bounds_wgsg84[1]), np.abs(bounds_wgsg84[3]))
-            if self.crs.is_projected and self.crs.axis_info[0].unit_name == "metre":
-                dx_m = obj.raster.res[0]
-            else:
-                dx_m = gis_utils.cellres(max_lat, *obj.raster.res)[0]
-            C = 2 * np.pi * 6378137  # circumference of the earth
+            # Determine the max number of zoom levels with the resolution
+            dres = tr_3857[0]
             max_lvl = int(
-                np.ceil(np.log2(C * np.cos(np.deg2rad(max_lat)) / dx_m / pxs))
+                math.ceil((math.log10((y_ext_pm * 2) / (dres * pxs)) / math.log10(2)))
             )
         if min_lvl is None:  # calculate min zoomlevel based on the data extent
-            min_lvl = mct.bounding_tile(*bounds_wgsg84).z
+            min_lvl = mct.bounding_tile(*bounds_4326).z
 
         # Loop through the zoom levels
         zoom_levels = {}
         logger.info(f"Producing tiles from zoomlevel {min_lvl} to {max_lvl}")
         for zl in range(max_lvl, min_lvl - 1, -1):
             fns = []
             # Go through the zoomlevels
-            for i, tile in enumerate(mct.tiles(*bounds_wgsg84, zl, truncate=True)):
+            for i, tile in enumerate(mct.tiles(*bounds_4326, zl, truncate=True)):
                 ssd = Path(root, str(zl), f"{tile.x}")
                 os.makedirs(ssd, exist_ok=True)
                 tile_bounds = mct.xy_bounds(tile)
@@ -2460,7 +2506,7 @@ def to_slippy_tiles(
                 "zoom_levels": zoom_levels,
             }
             name = os.path.basename(root)
-            with open(join(root, "data.yml"), "w") as f:
+            with open(join(root, f"{name}.yml"), "w") as f:
                 yaml.dump({name: yml}, f, default_flow_style=False, sort_keys=False)
 
     def to_raster(