From 9eaa6994158d393b246818baefd6e1b06e44b921 Mon Sep 17 00:00:00 2001 From: Brendan Date: Wed, 15 Feb 2023 10:50:17 +0100 Subject: [PATCH 01/35] Working on osm... --- hydromt/raster.py | 216 +++++++++++++++++++++++----------------------- 1 file changed, 109 insertions(+), 107 deletions(-) diff --git a/hydromt/raster.py b/hydromt/raster.py index 65fa70d75..31e0cec86 100644 --- a/hydromt/raster.py +++ b/hydromt/raster.py @@ -1964,113 +1964,115 @@ def tile_window(shape, px): with open(join(root, "..", f"{mName}.yml"), "w") as f: yaml.dump({mName: yml}, f, default_flow_style=False, sort_keys=False) - # def to_osm( - # self, - # root: str, - # zl: int, - # bbox: tuple = (), - # ): - # """Generate tiles from raster according to the osm scheme - - # Parameters - # ---------- - # root : str - # Path to folder where the database will be created - # zl : int - # Maximum zoom level of the database - # Everyting is generated incrementally up until this level - # E.g. zl = 8, levels generated: 0 to 7 - # bbox : tuple, optional - # Bounding Box in the objects crs - - # """ - - # assert self._obj.ndim == 2, "Only 2d datasets are accepted..." - # obj = self._obj.transpose(self.y_dim, self.x_dim) - - # mName = os.path.normpath(os.path.basename(root)) - - # def create_folder(path): - # if not os.path.exists(path): - # os.makedirs(path) - - # def transform_res(dres, transformer): - # return transformer.transform(0, dres)[0] - - # create_folder(root) - - # dres = abs(self._obj.raster.res[0]) - # if bbox: - # minx, miny, maxx, maxy = bbox - # else: - # minx, miny, maxx, maxy = self._obj.raster.transform_bounds( - # dst_crs=self._obj.raster.crs - # ) - - # transformer = pyproj.Transformer.from_crs(self._obj.raster.crs.to_epsg(), 3857) - # minx, miny = map( - # max, zip(transformer.transform(miny, minx), [-20037508.34] * 2) - # ) - # maxx, maxy = map(min, zip(transformer.transform(maxy, maxx), [20037508.34] * 2)) - - # dres = transform_res(dres, transformer) - # nzl = int(np.ceil((np.log10((20037508.34 * 2) / (dres * 256)) / np.log10(2)))) - - # if zl > nzl: - # zl = nzl - - # def tile_window(zl, minx, miny, maxx, maxy): - # # Basic stuff - # dx = (20037508.34 * 2) / (2**zl) - # # Origin displacement - # odx = np.floor(abs(-20037508.34 - minx) / dx) - # ody = np.floor(abs(20037508.34 - maxy) / dx) - - # # Set the new origin - # minx = -20037508.34 + odx * dx - # maxy = 20037508.34 - ody * dx - - # # Create window generator - # lu = product(np.arange(minx, maxx, dx), np.arange(maxy, miny, -dx)) - # for l, u in lu: - # col = int(odx + (l - minx) / dx) - # row = int(ody + (maxy - u) / dx) - # yield Affine(dx / 256, 0, l, 0, -dx / 256, u), col, row - - # for zlvl in range(zl): - # sd = f"{root}\\{zlvl}" - # create_folder(sd) - # file = open(f"{sd}\\filelist.txt", "w") - - # for transform, col, row in tile_window(zlvl, minx, miny, maxx, maxy): - # ssd = f"{sd}\\{col}" - # create_folder(ssd) - - # temp = obj.load() - # temp = temp.raster.reproject( - # dst_transform=transform, - # dst_crs=3857, - # dst_width=256, - # dst_height=256, - # ) - - # temp.raster.to_raster(f"{ssd}\\{row}.tif", driver="GTiff") - - # file.write(f"{ssd}\\{row}.tif\n") - - # del temp - - # file.close() - - # gis_utils.create_vrt(sd, mName) - # # Write a quick yaml for the database - # with open(f"{root}\\..\\{mName}.yml", "w") as w: - # w.write(f"{mName}:\n") - # crs = 3857 - # w.write(f" crs: {crs}\n") - # w.write(" data_type: RasterDataset\n") - # w.write(" driver: raster\n") - # w.write(f" path: {mName}/{{zoom_level}}/{mName}.vrt\n") + def to_osm( + self, + root: str, + zl: int, + bbox: tuple = (), + ): + """Generate tiles from raster according to the osm scheme + + Parameters + ---------- + root : str + Path to folder where the database will be created + zl : int + Maximum zoom level of the database + Everyting is generated incrementally up until this level + E.g. zl = 8, levels generated: 0 to 7 + bbox : tuple, optional + Bounding Box in the objects crs + + """ + + assert self._obj.ndim == 2, "Only 2d datasets are accepted..." + obj = self._obj.copy() + obj = obj.transpose(self.y_dim, self.x_dim) + obj = obj.raster.reproject(dst_crs=3857) + + mName = os.path.normpath(os.path.basename(root)) + + def create_folder(path): + if not os.path.exists(path): + os.makedirs(path) + + def transform_res(dres, transformer): + return transformer.transform(0, dres)[0] + + create_folder(root) + + dres = abs(self._obj.raster.res[0]) + if bbox: + minx, miny, maxx, maxy = bbox + else: + minx, miny, maxx, maxy = self._obj.raster.transform_bounds( + dst_crs=self._obj.raster.crs + ) + + transformer = pyproj.Transformer.from_crs(self._obj.raster.crs.to_epsg(), 3857) + minx, miny = map( + max, zip(transformer.transform(miny, minx), [-20037508.34] * 2) + ) + maxx, maxy = map(min, zip(transformer.transform(maxy, maxx), [20037508.34] * 2)) + + dres = transform_res(dres, transformer) + nzl = int(np.ceil((np.log10((20037508.34 * 2) / (dres * 256)) / np.log10(2)))) + + if zl > nzl: + zl = nzl + + def tile_window(zl, minx, miny, maxx, maxy): + # Basic stuff + dx = (20037508.34 * 2) / (2**zl) + # Origin displacement + odx = np.floor(abs(-20037508.34 - minx) / dx) + ody = np.floor(abs(20037508.34 - maxy) / dx) + + # Set the new origin + minx = -20037508.34 + odx * dx + maxy = 20037508.34 - ody * dx + + # Create window generator + lu = product(np.arange(minx, maxx, dx), np.arange(maxy, miny, -dx)) + for l, u in lu: + col = int(odx + (l - minx) / dx) + row = int(ody + (maxy - u) / dx) + yield Affine(dx / 256, 0, l, 0, -dx / 256, u), col, row + + for zlvl in range(zl): + sd = f"{root}\\{zlvl}" + create_folder(sd) + file = open(f"{sd}\\filelist.txt", "w") + + for transform, col, row in tile_window(zlvl, minx, miny, maxx, maxy): + ssd = f"{sd}\\{col}" + create_folder(ssd) + + temp = obj.load() + temp = temp.raster.reproject( + dst_transform=transform, + dst_crs=3857, + dst_width=256, + dst_height=256, + ) + + temp.raster.to_raster(f"{ssd}\\{row}.tif", driver="GTiff") + + file.write(f"{ssd}\\{row}.tif\n") + + del temp + + file.close() + + gis_utils.create_vrt(sd, mName) + # Write a quick yaml for the database + with open(f"{root}\\..\\{mName}.yml", "w") as w: + w.write(f"{mName}:\n") + crs = 3857 + w.write(f" crs: {crs}\n") + w.write(" data_type: RasterDataset\n") + w.write(" driver: raster\n") + w.write(f" path: {mName}/{{zoom_level}}/{mName}.vrt\n") def to_raster( self, From 955e859b6f11161eff03d86b81641df086605ea8 Mon Sep 17 00:00:00 2001 From: Brendan Date: Wed, 12 Jul 2023 15:01:53 +0200 Subject: [PATCH 02/35] Added 'to_osm_tiles' --- docs/api.rst | 1 + docs/changelog.rst | 1 + examples/tmpdir/.gitkeep | 0 examples/working_with_tiled_raster_data.ipynb | 68 +++++- hydromt/raster.py | 226 ++++++++++++++---- 5 files changed, 244 insertions(+), 52 deletions(-) delete mode 100644 examples/tmpdir/.gitkeep diff --git a/docs/api.rst b/docs/api.rst index 97f60dc04..f86e92d7e 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -459,6 +459,7 @@ Raster writing methods :toctree: _generated DataArray.raster.to_xyz_tiles + DataArray.raster.to_osm_tiles DataArray.raster.to_raster Dataset.raster.to_mapstack diff --git a/docs/changelog.rst b/docs/changelog.rst index ebac50351..70e446487 100644 --- a/docs/changelog.rst +++ b/docs/changelog.rst @@ -14,6 +14,7 @@ Added - New generic methods for ``GridModel``: ``setup_grid``, ``setup_grid_from_constant``, ``setup_grid_from_rasterdataset``, ``setup_grid_from_raster_reclass``, ``setup_grid_from_geodataframe``. PR #333 - New ``grid`` workflow methods to support the setup methods in ``GridModel``: ``grid_from_constant``, ``grid_from_rasterdataset``, ``grid_from_raster_reclass``, ``grid_from_geodataframe``. PR #333 - New raster method ``rasterize_geometry``. +- New raster method ``to_osm_tiles``: tiling of a raster dataset according to the osm structure. Changed ------- diff --git a/examples/tmpdir/.gitkeep b/examples/tmpdir/.gitkeep deleted file mode 100644 index e69de29bb..000000000 diff --git a/examples/working_with_tiled_raster_data.ipynb b/examples/working_with_tiled_raster_data.ipynb index 2394e468f..7684192de 100644 --- a/examples/working_with_tiled_raster_data.ipynb +++ b/examples/working_with_tiled_raster_data.ipynb @@ -52,7 +52,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### tiling raster with XYZ stucture\n", + "### Tiling raster with XYZ stucture\n", "\n", "First let's have a look at the XYZ structure.\n", "an xarray.DataArray is simple written to a tile database in XYZ structure via .raster.to_xyz_tiles\n", @@ -92,12 +92,55 @@ "At last, a .yml file is produced which can be read by the [DataCatalog](https://deltares.github.io/hydromt/latest/_generated/hydromt.data_catalog.DataCatalog.html) of HydroMT." ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Tiling raster with OSM stucture\n", + "\n", + "Now let's have a look at tiling according to the OSM structure\n", + "an xarray.DataArray is simple written to a tile database in OSM structure via .raster.to_osm_tiles" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Write the database in XYZ stucture\n", + "name_osm = f\"{source}_osm\"\n", + "root_osm = join(\"tmpdir\", name_osm)\n", + "zoom_levels = [0, 1, 2, 5] # note that you have to present all levels\n", + "da0.raster.to_osm_tiles(\n", + " root=root_osm,\n", + " zoom_levels=zoom_levels,\n", + " min_lvl=9, # minimal level from where to upscale from (prevent memory issues)\n", + " driver=\"GTiff\", # try also 'netcdf4'\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The tiles in the 'merit_hydro_OSM' folder now contains all the zoom_levels from the minimal starting level (9) until the highest level (0).\n", + "\n", + "Every tile, regardless of the zoomlevel, has a resolution of 256 by 256 pixels.\n", + "\n", + "Zoomlevel 0 is at the scale of the entire world (one tile) and is the most downscaled. Zoomlevel 9 contains the highest resolution (most tiles) in regards to this tile database.\n", + "\n", + "A mosaic is created per zoomlevel of these tiles in a .vrt file.\n", + "\n", + "At last, a .yml file is produced which can be read by the [DataCatalog](https://deltares.github.io/hydromt/latest/_generated/hydromt.data_catalog.DataCatalog.html) of HydroMT." + ] + }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ - "### reading tiled raster data with zoom levels\n", + "### Reading tiled raster data with zoom levels\n", "\n", "With DataCatalog.[get_rasterdataset](https://deltares.github.io/hydromt/latest/_generated/hydromt.data_catalog.DataCatalog.get_rasterdataset.html) a raster (.vrt) can be retrieved. In case of a tile database it can be done for a certain zoomlevel. E.g." ] @@ -112,10 +155,10 @@ "from hydromt import DataCatalog\n", "\n", "# Load the yml into a DataCatalog\n", - "data_catalog = DataCatalog(join(root, f\"{name}.yml\"), logger=logger)\n", + "data_catalog = DataCatalog([join(root, f\"{name}.yml\"), join(root_osm, f\"{name_osm}.yml\")], logger=logger)\n", "\n", "# View the structure of the DataCatalog\n", - "data_catalog[name]" + "data_catalog" ] }, { @@ -129,6 +172,17 @@ "da0.raster.shape" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# And for OSM\n", + "da1 = data_catalog.get_rasterdataset(name_osm)\n", + "da1.raster.shape" + ] + }, { "cell_type": "code", "execution_count": null, @@ -137,8 +191,10 @@ "source": [ "# Request a raster from the Datacatolog based on zoom resolution & unit\n", "da = data_catalog.get_rasterdataset(name, zoom_level=(1/600, 'degree'))\n", + "da = data_catalog.get_rasterdataset(name_osm, zoom_level=(1/600, 'degree'))\n", "\n", - "da = data_catalog.get_rasterdataset(name, zoom_level=(1e3, 'meter'))\n" + "da = data_catalog.get_rasterdataset(name, zoom_level=(1e3, 'degree'))\n", + "da = data_catalog.get_rasterdataset(name_osm, zoom_level=(1e4, 'meter'))\n" ] }, { @@ -220,7 +276,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.9" + "version": "3.11.4" }, "orig_nbformat": 4, "vscode": { diff --git a/hydromt/raster.py b/hydromt/raster.py index f4b59aa73..124cb7205 100644 --- a/hydromt/raster.py +++ b/hydromt/raster.py @@ -9,6 +9,7 @@ import logging import math +import gc import os import tempfile from itertools import product @@ -2233,31 +2234,48 @@ def tile_window(shape, px): with open(join(root, f"{mName}.yml"), "w") as f: yaml.dump({mName: yml}, f, default_flow_style=False, sort_keys=False) - def to_osm( + def to_osm_tiles( self, root: str, - zl: int, + zoom_levels: list, bbox: tuple = (), + method: str = "nearest", + min_lvl: int = 5, + driver="GTiff", + **kwargs, ): - """Generate tiles from raster according to the osm scheme + """Export rasterdataset to tiles in an osm structure. Parameters ---------- root : str - Path to folder where the database will be created - zl : int - Maximum zoom level of the database - Everyting is generated incrementally up until this level - E.g. zl = 8, levels generated: 0 to 7 + Path where the database will be saved + Database yml will be put one directory above + zoom_levels : list + Zoom levels to be put in the database bbox : tuple, optional - Bounding Box in the objects crs - + A region where the tiles should be generated for + (if absent, the entire input raster dataset will be tiled) + method : str, optional + How to resample the data when downscaling. + E.g. 'nearest' for resampling with the nearest value + min_lvl : int, optional + Minimal starting level from to where to downscale from. + Increase if memory issues are encountered. + driver : str, optional + GDAL driver (e.g., 'GTiff' for geotif files), or 'netcdf4' for netcdf files. + **kwargs + Key-word arguments to write raster files """ + resampling = getattr(Resampling, method, None) + if resampling is None: + raise ValueError(f"Resampling method unknown: {method}.") + px_size = 256 # Fixed resolution for osm + assert self._obj.ndim == 2, "Only 2d datasets are accepted..." obj = self._obj.copy() obj = obj.transpose(self.y_dim, self.x_dim) - obj = obj.raster.reproject(dst_crs=3857) mName = os.path.normpath(os.path.basename(root)) @@ -2268,31 +2286,30 @@ def create_folder(path): def transform_res(dres, transformer): return transformer.transform(0, dres)[0] - create_folder(root) + transformer = pyproj.Transformer.from_crs(self._obj.raster.crs.to_epsg(), 3857) - dres = abs(self._obj.raster.res[0]) if bbox: minx, miny, maxx, maxy = bbox else: - minx, miny, maxx, maxy = self._obj.raster.transform_bounds( - dst_crs=self._obj.raster.crs - ) + minx, miny, maxx, maxy = obj.raster.transform_bounds("EPSG:3857") - transformer = pyproj.Transformer.from_crs(self._obj.raster.crs.to_epsg(), 3857) minx, miny = map( - max, zip(transformer.transform(miny, minx), [-20037508.34] * 2) + max, + zip((minx, miny), [-20037508.34] * 2), + ) + maxx, maxy = map( + min, + zip((maxx, maxy), [20037508.34] * 2), ) - maxx, maxy = map(min, zip(transformer.transform(maxy, maxx), [20037508.34] * 2)) + dres = abs(self._obj.raster.res[0]) dres = transform_res(dres, transformer) nzl = int(np.ceil((np.log10((20037508.34 * 2) / (dres * 256)) / np.log10(2)))) - if zl > nzl: - zl = nzl + if max(zoom_levels) > nzl: + raise OSError("") def tile_window(zl, minx, miny, maxx, maxy): - # Basic stuff - dx = (20037508.34 * 2) / (2**zl) # Origin displacement odx = np.floor(abs(-20037508.34 - minx) / dx) ody = np.floor(abs(20037508.34 - maxy) / dx) @@ -2302,46 +2319,163 @@ def tile_window(zl, minx, miny, maxx, maxy): maxy = 20037508.34 - ody * dx # Create window generator - lu = product(np.arange(minx, maxx, dx), np.arange(maxy, miny, -dx)) + lu = product( + np.arange(round(minx, 2), round(maxx, 2), dx), + np.arange(round(maxy, 2), round(miny, 2), -dx), + ) for l, u in lu: - col = int(odx + (l - minx) / dx) - row = int(ody + (maxy - u) / dx) - yield Affine(dx / 256, 0, l, 0, -dx / 256, u), col, row + col = int(round(odx + (l - minx) / dx, 2)) + row = int(round(ody + (maxy - u) / dx, 2)) + bounds = [ + l, + u - dx, + l + dx, + u, + ] + yield Affine(px, 0, l, 0, -px, u), bounds, col, row + + zoom_levels.sort() + if zoom_levels[-1] > min_lvl: + min_lvl = zoom_levels[-1] + zoom_levels = list(range(zoom_levels[0], min_lvl + 1, 1)) + zoom_levels.reverse() + + create_folder(root) + + vrt_fn = None + nodata = self.nodata + prev = 0 + zls = {} + for zl in zoom_levels: + diff = abs(zl - prev) + dx = (20037508.34 * 2) / (2**zl) + px = dx / px_size + + # read data from previous zoomlevel + if vrt_fn is not None: + obj = xr.open_dataarray(vrt_fn, engine="rasterio").squeeze( + "band", drop=True + ) + x_dim, y_dim = obj.raster.x_dim, obj.raster.y_dim + obj = obj.chunk({x_dim: 1024, y_dim: 1024}) + + obj_bounds = obj.raster.bounds + obj_res = obj.raster.res[0] + obj_shape = obj.raster.shape - for zlvl in range(zl): - sd = f"{root}\\{zlvl}" + t_r = pyproj.Transformer.from_crs( + 3857, obj.raster.crs.to_epsg(), always_xy=True + ) + + sd = f"{root}\\{zl}" create_folder(sd) - file = open(f"{sd}\\filelist.txt", "w") + txt_path = join(sd, "filelist.txt") + file = open(txt_path, "w") - for transform, col, row in tile_window(zlvl, minx, miny, maxx, maxy): + for transform, bounds, col, row in tile_window(zl, minx, miny, maxx, maxy): ssd = f"{sd}\\{col}" create_folder(ssd) - temp = obj.load() - temp = temp.raster.reproject( + ul = t_r.transform(bounds[0], bounds[3]) + lr = t_r.transform(bounds[2], bounds[1]) + out_bounds = [ul[0], lr[1], lr[0], ul[1]] + + temp_bounds = [ + max(out_bounds[0], obj_bounds[0]), + max(out_bounds[1], obj_bounds[1]), + min(out_bounds[2], obj_bounds[2]), + min(out_bounds[3], obj_bounds[3]), + ] + + _window = [ + abs(x - y) / obj_res for x, y in zip(temp_bounds, obj_bounds) + ] + + u = math.floor(round(_window[3], 2)) + l = math.floor(round(_window[0], 2)) + h = obj_shape[0] - u - math.floor(round(_window[1], 2)) + w = obj_shape[1] - l - math.floor(round(_window[2], 2)) + + temp = obj[u : u + h, l : l + w] + + tile = np.empty((px_size, px_size), np.float64) + tile[:] = nodata + + temp_transform = Affine( + obj_res, + 0, + obj_bounds[0] + l * obj_res, + 0, + -obj_res, + obj_bounds[3] - u * obj_res, + ) + + rasterio.warp.reproject( + temp.values, + tile, + src_transform=temp_transform, + src_crs=obj.raster.crs, dst_transform=transform, - dst_crs=3857, - dst_width=256, - dst_height=256, + dst_crs="EPSG:3857", + dst_nodata=nodata, + resampling=resampling, ) - temp.raster.to_raster(f"{ssd}\\{row}.tif", driver="GTiff") + temp = xr.DataArray( + name=temp.name, + dims=("y", "x"), + coords={ + "y": np.arange( + bounds[3] + 0.5 * transform[4], bounds[1], transform[4] + )[:px_size], + "x": np.arange( + bounds[0] + 0.5 * transform[0], + bounds[2] + 0.5 * transform[0], + transform[0], + )[:px_size], + }, + data=tile, + ) + + del tile + + temp.raster.set_crs(3857) - file.write(f"{ssd}\\{row}.tif\n") + temp.raster.set_nodata(nodata) + + if driver == "netcdf4": + path = join(ssd, f"{row}.nc") + temp = temp.raster.gdal_compliant() + temp.to_netcdf(path, engine="netcdf4", **kwargs) + elif driver in gis_utils.GDAL_EXT_CODE_MAP: + ext = gis_utils.GDAL_EXT_CODE_MAP.get(driver) + path = join(ssd, f"{row}.{ext}") + temp.raster.to_raster(path, driver=driver, **kwargs) + else: + raise ValueError(f"Unkown file driver {driver}") + file.write(f"{path}\n") del temp + gc.collect() file.close() + # Create a vrt using GDAL + vrt_fn = join(root, f"{mName}_zl{zl}.vrt") + gis_utils.create_vrt(vrt_fn, file_list_path=txt_path) + prev = zl + zls.update({zl: round(float(px), 2)}) + del obj - gis_utils.create_vrt(sd, mName) # Write a quick yaml for the database - with open(f"{root}\\..\\{mName}.yml", "w") as w: - w.write(f"{mName}:\n") - crs = 3857 - w.write(f" crs: {crs}\n") - w.write(" data_type: RasterDataset\n") - w.write(" driver: raster\n") - w.write(f" path: {mName}/{{zoom_level}}/{mName}.vrt\n") + yml = { + "crs": 3857, + "data_type": "RasterDataset", + "driver": "raster", + "path": f"{mName}_zl{{zoom_level}}.vrt", + "zoom_levels": zls, + } + with open(join(root, f"{mName}.yml"), "w") as f: + yaml.dump({mName: yml}, f, default_flow_style=False, sort_keys=False) def to_raster( self, From 9bc42d0022032e4d152ac1266cf824779cc66f31 Mon Sep 17 00:00:00 2001 From: Brendan Date: Wed, 12 Jul 2023 15:19:41 +0200 Subject: [PATCH 03/35] Cleared output, deleted unused variable --- examples/prep_data_catalog.ipynb | 2 +- hydromt/raster.py | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/examples/prep_data_catalog.ipynb b/examples/prep_data_catalog.ipynb index ad465b72f..c74066a24 100644 --- a/examples/prep_data_catalog.ipynb +++ b/examples/prep_data_catalog.ipynb @@ -773,7 +773,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.6" + "version": "3.11.4" }, "vscode": { "interpreter": { diff --git a/hydromt/raster.py b/hydromt/raster.py index 124cb7205..bcba499fd 100644 --- a/hydromt/raster.py +++ b/hydromt/raster.py @@ -2347,7 +2347,6 @@ def tile_window(zl, minx, miny, maxx, maxy): prev = 0 zls = {} for zl in zoom_levels: - diff = abs(zl - prev) dx = (20037508.34 * 2) / (2**zl) px = dx / px_size From 951d3e1f92af20c0ff32dead45a7510da08a476b Mon Sep 17 00:00:00 2001 From: Brendan Date: Wed, 12 Jul 2023 15:22:09 +0200 Subject: [PATCH 04/35] Restored showing of DataCatalog --- examples/working_with_tiled_raster_data.ipynb | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/examples/working_with_tiled_raster_data.ipynb b/examples/working_with_tiled_raster_data.ipynb index 7684192de..b5e6b7ada 100644 --- a/examples/working_with_tiled_raster_data.ipynb +++ b/examples/working_with_tiled_raster_data.ipynb @@ -158,7 +158,17 @@ "data_catalog = DataCatalog([join(root, f\"{name}.yml\"), join(root_osm, f\"{name_osm}.yml\")], logger=logger)\n", "\n", "# View the structure of the DataCatalog\n", - "data_catalog" + "data_catalog[name]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# For the osm dataset\n", + "data_catalog[name_osm]" ] }, { From c1ab0e3167005d62519ab6ee3c85261bb794e293 Mon Sep 17 00:00:00 2001 From: Brendan Date: Wed, 12 Jul 2023 15:24:20 +0200 Subject: [PATCH 05/35] Linting... --- hydromt/raster.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/hydromt/raster.py b/hydromt/raster.py index bcba499fd..878328ae3 100644 --- a/hydromt/raster.py +++ b/hydromt/raster.py @@ -7,9 +7,9 @@ """Extension for xarray to provide rasterio capabilities to xarray datasets/arrays.""" from __future__ import annotations +import gc import logging import math -import gc import os import tempfile from itertools import product @@ -2344,7 +2344,6 @@ def tile_window(zl, minx, miny, maxx, maxy): vrt_fn = None nodata = self.nodata - prev = 0 zls = {} for zl in zoom_levels: dx = (20037508.34 * 2) / (2**zl) @@ -2461,7 +2460,6 @@ def tile_window(zl, minx, miny, maxx, maxy): # Create a vrt using GDAL vrt_fn = join(root, f"{mName}_zl{zl}.vrt") gis_utils.create_vrt(vrt_fn, file_list_path=txt_path) - prev = zl zls.update({zl: round(float(px), 2)}) del obj From ed12d3df95cc33337b98bc732730daa6c096a8cb Mon Sep 17 00:00:00 2001 From: Brendan Date: Wed, 12 Jul 2023 15:27:05 +0200 Subject: [PATCH 06/35] Removed whitespace --- hydromt/raster.py | 1 - 1 file changed, 1 deletion(-) diff --git a/hydromt/raster.py b/hydromt/raster.py index 878328ae3..87dde28e8 100644 --- a/hydromt/raster.py +++ b/hydromt/raster.py @@ -2267,7 +2267,6 @@ def to_osm_tiles( **kwargs Key-word arguments to write raster files """ - resampling = getattr(Resampling, method, None) if resampling is None: raise ValueError(f"Resampling method unknown: {method}.") From 149d4203c2155da817c66a2434b0de861b4efb75 Mon Sep 17 00:00:00 2001 From: Brendan Date: Wed, 12 Jul 2023 15:33:11 +0200 Subject: [PATCH 07/35] Readded tmpdir --- examples/tmpdir/.gitkeep | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 examples/tmpdir/.gitkeep diff --git a/examples/tmpdir/.gitkeep b/examples/tmpdir/.gitkeep new file mode 100644 index 000000000..e69de29bb From da11a4eeb29b0d663da1022cd54e4540b3ff292c Mon Sep 17 00:00:00 2001 From: Brendan Date: Wed, 12 Jul 2023 15:37:31 +0200 Subject: [PATCH 08/35] Added OSM tiling test --- tests/test_raster.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/tests/test_raster.py b/tests/test_raster.py index 942d425c8..05e557ba2 100644 --- a/tests/test_raster.py +++ b/tests/test_raster.py @@ -500,5 +500,14 @@ def test_to_xyz_tiles(tmpdir, rioda_large): assert len(f.readlines()) == 1 -# def test_to_osm(tmpdir, dummy): -# dummy.raster.to_osm( +def test_to_osm_tiles(tmpdir, rioda_large): + path = str(tmpdir) + rioda_large.raster.to_osm_tiles( + os.path.join(path, "dummy_osm"), + zoom_levels=[5], + min_lvl=8, + ) + with open(os.path.join(path, "dummy_osm", "5", "filelist.txt"), "r") as f: + assert len(f.readlines()) == 1 + with open(os.path.join(path, "dummy_osm", "8", "filelist.txt"), "r") as f: + assert len(f.readlines()) == 12 From bb0827313322ff66c99702d66e0adea1f91e4932 Mon Sep 17 00:00:00 2001 From: Brendan Date: Wed, 12 Jul 2023 15:45:58 +0200 Subject: [PATCH 09/35] No concats in paths.. --- hydromt/raster.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/hydromt/raster.py b/hydromt/raster.py index 87dde28e8..a7d8a7da4 100644 --- a/hydromt/raster.py +++ b/hydromt/raster.py @@ -2364,13 +2364,13 @@ def tile_window(zl, minx, miny, maxx, maxy): 3857, obj.raster.crs.to_epsg(), always_xy=True ) - sd = f"{root}\\{zl}" + sd = join(root, f"{zl}") create_folder(sd) txt_path = join(sd, "filelist.txt") file = open(txt_path, "w") for transform, bounds, col, row in tile_window(zl, minx, miny, maxx, maxy): - ssd = f"{sd}\\{col}" + ssd = join(sd, f"{col}") create_folder(ssd) ul = t_r.transform(bounds[0], bounds[3]) From f884bbf4570609172606adb36c018bc43e7ac7bb Mon Sep 17 00:00:00 2001 From: Brendan Date: Thu, 13 Jul 2023 09:57:22 +0200 Subject: [PATCH 10/35] Chunking like this blows up memory, removed it --- hydromt/raster.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/hydromt/raster.py b/hydromt/raster.py index a7d8a7da4..469c50a5c 100644 --- a/hydromt/raster.py +++ b/hydromt/raster.py @@ -2353,8 +2353,6 @@ def tile_window(zl, minx, miny, maxx, maxy): obj = xr.open_dataarray(vrt_fn, engine="rasterio").squeeze( "band", drop=True ) - x_dim, y_dim = obj.raster.x_dim, obj.raster.y_dim - obj = obj.chunk({x_dim: 1024, y_dim: 1024}) obj_bounds = obj.raster.bounds obj_res = obj.raster.res[0] From 7836a401b3590189f76b3aa7bcfebdb01252367f Mon Sep 17 00:00:00 2001 From: Brendan Date: Fri, 14 Jul 2023 19:37:10 +0200 Subject: [PATCH 11/35] Properly induce the reprojected pixel size of the source data --- hydromt/raster.py | 22 ++++++++++++---------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/hydromt/raster.py b/hydromt/raster.py index 469c50a5c..68b243cc8 100644 --- a/hydromt/raster.py +++ b/hydromt/raster.py @@ -2282,11 +2282,6 @@ def create_folder(path): if not os.path.exists(path): os.makedirs(path) - def transform_res(dres, transformer): - return transformer.transform(0, dres)[0] - - transformer = pyproj.Transformer.from_crs(self._obj.raster.crs.to_epsg(), 3857) - if bbox: minx, miny, maxx, maxy = bbox else: @@ -2301,13 +2296,15 @@ def transform_res(dres, transformer): zip((maxx, maxy), [20037508.34] * 2), ) - dres = abs(self._obj.raster.res[0]) - dres = transform_res(dres, transformer) + _obj_tr = rasterio.warp.calculate_default_transform( + obj.raster.crs, + "EPSG:3857", + *obj.shape, + *obj.raster.bounds, + )[0] + dres = _obj_tr[0] nzl = int(np.ceil((np.log10((20037508.34 * 2) / (dres * 256)) / np.log10(2)))) - if max(zoom_levels) > nzl: - raise OSError("") - def tile_window(zl, minx, miny, maxx, maxy): # Origin displacement odx = np.floor(abs(-20037508.34 - minx) / dx) @@ -2348,6 +2345,11 @@ def tile_window(zl, minx, miny, maxx, maxy): dx = (20037508.34 * 2) / (2**zl) px = dx / px_size + if max(zoom_levels) > nzl: + logger.warning( + f"Pixels at zoomlevel {zl} are smaller than the original data" + ) + # read data from previous zoomlevel if vrt_fn is not None: obj = xr.open_dataarray(vrt_fn, engine="rasterio").squeeze( From 1e8fa3f03d8818b8506884eac21ad5a1d0d27fde Mon Sep 17 00:00:00 2001 From: Brendan Date: Wed, 26 Jul 2023 10:30:46 +0200 Subject: [PATCH 12/35] Updated changelog post v0.8.0 release --- docs/changelog.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/changelog.rst b/docs/changelog.rst index d8a07cbfb..7b9e578e0 100644 --- a/docs/changelog.rst +++ b/docs/changelog.rst @@ -11,7 +11,7 @@ Unreleased Added ----- -- +- New raster method ``to_osm_tiles``: tiling of a raster dataset according to the osm structure. Changed ------- From a85d4588f312081b71137291d5ca4d10548535e9 Mon Sep 17 00:00:00 2001 From: Brendan Date: Fri, 11 Aug 2023 13:26:04 +0200 Subject: [PATCH 13/35] snake_case to_osm --- hydromt/raster.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/hydromt/raster.py b/hydromt/raster.py index b505ced3b..635a4c520 100644 --- a/hydromt/raster.py +++ b/hydromt/raster.py @@ -2290,7 +2290,7 @@ def to_osm_tiles( obj = self._obj.copy() obj = obj.transpose(self.y_dim, self.x_dim) - mName = os.path.normpath(os.path.basename(root)) + m_name = os.path.normpath(os.path.basename(root)) def create_folder(path): if not os.path.exists(path): @@ -2471,7 +2471,7 @@ def tile_window(zl, minx, miny, maxx, maxy): file.close() # Create a vrt using GDAL - vrt_fn = join(root, f"{mName}_zl{zl}.vrt") + vrt_fn = join(root, f"{m_name}_zl{zl}.vrt") gis_utils.create_vrt(vrt_fn, file_list_path=txt_path) zls.update({zl: round(float(px), 2)}) del obj @@ -2481,11 +2481,11 @@ def tile_window(zl, minx, miny, maxx, maxy): "crs": 3857, "data_type": "RasterDataset", "driver": "raster", - "path": f"{mName}_zl{{zoom_level}}.vrt", + "path": f"{m_name}_zl{{zoom_level}}.vrt", "zoom_levels": zls, } - with open(join(root, f"{mName}.yml"), "w") as f: - yaml.dump({mName: yml}, f, default_flow_style=False, sort_keys=False) + with open(join(root, f"{m_name}.yml"), "w") as f: + yaml.dump({m_name: yml}, f, default_flow_style=False, sort_keys=False) def to_raster( self, From dc01c39fe5668b84858b841cec9073d3160c4e27 Mon Sep 17 00:00:00 2001 From: Brendan Date: Fri, 11 Aug 2023 13:34:06 +0200 Subject: [PATCH 14/35] moved 'create_folder' --- hydromt/raster.py | 9 +-------- hydromt/utils.py | 6 ++++++ 2 files changed, 7 insertions(+), 8 deletions(-) diff --git a/hydromt/raster.py b/hydromt/raster.py index ce892d119..60d769dd1 100644 --- a/hydromt/raster.py +++ b/hydromt/raster.py @@ -37,6 +37,7 @@ from shapely.geometry import LineString, Polygon, box from . import _compat, gis_utils +from .utils import create_folder logger = logging.getLogger(__name__) XDIMS = ("x", "longitude", "lon", "long") @@ -2155,10 +2156,6 @@ def to_xyz_tiles( """ mName = os.path.normpath(os.path.basename(root)) - def create_folder(path): - if not os.path.exists(path): - os.makedirs(path) - def tile_window(shape, px): """Yield (left, upper, width, height).""" nr, nc = shape @@ -2292,10 +2289,6 @@ def to_osm_tiles( m_name = os.path.normpath(os.path.basename(root)) - def create_folder(path): - if not os.path.exists(path): - os.makedirs(path) - if bbox: minx, miny, maxx, maxy = bbox else: diff --git a/hydromt/utils.py b/hydromt/utils.py index bb4d6d0db..451a34893 100644 --- a/hydromt/utils.py +++ b/hydromt/utils.py @@ -1,4 +1,5 @@ """Utility functions for hydromt that have no other home.""" +import os def partition_dictionaries(left, right): @@ -36,6 +37,11 @@ def partition_dictionaries(left, right): return common, left_less_right, right_less_left +def create_folder(path): + if not os.path.exists(path): + os.makedirs(path) + + def _dict_pprint(d): import json From c15ee9b5e06d18f04a2b1d26aa34323dab98301d Mon Sep 17 00:00:00 2001 From: Brendan Date: Fri, 11 Aug 2023 14:09:36 +0200 Subject: [PATCH 15/35] small osm fixes --- hydromt/raster.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/hydromt/raster.py b/hydromt/raster.py index 60d769dd1..90a4bbe9e 100644 --- a/hydromt/raster.py +++ b/hydromt/raster.py @@ -2402,8 +2402,7 @@ def tile_window(zl, minx, miny, maxx, maxy): temp = obj[u : u + h, l : l + w] - tile = np.empty((px_size, px_size), np.float64) - tile[:] = nodata + tile = np.full((px_size, px_size), nodata, dtype=np.float64) temp_transform = Affine( obj_res, @@ -2438,7 +2437,7 @@ def tile_window(zl, minx, miny, maxx, maxy): transform[0], )[:px_size], }, - data=tile, + data=tile.copy(), ) del tile From d21e47a481eec07abc4809d9a78b72a87d36fcaa Mon Sep 17 00:00:00 2001 From: Brendan Date: Fri, 11 Aug 2023 14:56:51 +0200 Subject: [PATCH 16/35] Added small docstring to 'create_folder' --- hydromt/utils.py | 1 + 1 file changed, 1 insertion(+) diff --git a/hydromt/utils.py b/hydromt/utils.py index 451a34893..2ae41a01c 100644 --- a/hydromt/utils.py +++ b/hydromt/utils.py @@ -38,6 +38,7 @@ def partition_dictionaries(left, right): def create_folder(path): + """Create a folder if there is none""" if not os.path.exists(path): os.makedirs(path) From 3ed491a9662f9038766c9da16054e13d3d68ec1b Mon Sep 17 00:00:00 2001 From: Brendan Date: Fri, 11 Aug 2023 15:08:29 +0200 Subject: [PATCH 17/35] Added a period (apparently...) --- hydromt/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hydromt/utils.py b/hydromt/utils.py index 2ae41a01c..b133a29ba 100644 --- a/hydromt/utils.py +++ b/hydromt/utils.py @@ -38,7 +38,7 @@ def partition_dictionaries(left, right): def create_folder(path): - """Create a folder if there is none""" + """Create a folder if there is none.""" if not os.path.exists(path): os.makedirs(path) From 670710716097d29d176a932a97af8c770a54df53 Mon Sep 17 00:00:00 2001 From: Brendan Date: Thu, 14 Sep 2023 17:11:19 +0200 Subject: [PATCH 18/35] Improved readability; added a tiling function that caters to webviewers --- docs/api.rst | 1 + docs/changelog.rst | 1 + examples/working_with_tiled_raster_data.ipynb | 68 ++- hydromt/raster.py | 437 +++++++++++++++--- hydromt/utils.py | 17 + tests/test_raster.py | 23 + 6 files changed, 487 insertions(+), 60 deletions(-) diff --git a/docs/api.rst b/docs/api.rst index 3be0d5051..7f57bb7b5 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -466,6 +466,7 @@ Raster writing methods DataArray.raster.to_xyz_tiles DataArray.raster.to_osm_tiles + DataArray.raster.to_webviewer_tiles DataArray.raster.to_raster Dataset.raster.to_mapstack diff --git a/docs/changelog.rst b/docs/changelog.rst index 666bf0945..d015d72f7 100644 --- a/docs/changelog.rst +++ b/docs/changelog.rst @@ -12,6 +12,7 @@ Unreleased Added ----- - New raster method ``to_osm_tiles``: tiling of a raster dataset according to the osm structure. +- New raster method ``to_webviewer_tiles``: tiling of a raster dataset for webviewers (png files) - docs now include a dropdown for selecting older versions of the docs. (#457) - Support for loading the same data source but from different places (e.g. local & aws) - Add support for reading and writing tabular data in ``parquet`` format. (PR #445) diff --git a/examples/working_with_tiled_raster_data.ipynb b/examples/working_with_tiled_raster_data.ipynb index b5e6b7ada..be41b6a52 100644 --- a/examples/working_with_tiled_raster_data.ipynb +++ b/examples/working_with_tiled_raster_data.ipynb @@ -135,6 +135,72 @@ "At last, a .yml file is produced which can be read by the [DataCatalog](https://deltares.github.io/hydromt/latest/_generated/hydromt.data_catalog.DataCatalog.html) of HydroMT." ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Tiling raster for a webviewer\n", + "\n", + "Finally, let's take a look at tiling of a raster dataset with its use being to view the data in a webviewer.\n", + "This is easily done with the .raster.to_webviewer_tiles method. \n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Write the data in OSM structure, but to images!\n", + "name_png = f\"{source}_png\"\n", + "root_png = join(\"tmpdir\", name_png)\n", + "da0.raster.to_webviewer_tiles(\n", + " root = root_png\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now that the images are created, let's take a look at an individual tile." + ] + }, + { + "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, \"9\", \"273\", \"182.png\"))\n", + "ax.imshow(im)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The image itself tends to look like an oil painting, but this is correct.\n", + "\n", + "The colors for this image are determined so that they are correctly visually represented according to Terrarium Terrain RGB.\n", + "\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" + ] + }, { "attachments": {}, "cell_type": "markdown", @@ -286,7 +352,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.4" + "version": "3.11.5" }, "orig_nbformat": 4, "vscode": { diff --git a/hydromt/raster.py b/hydromt/raster.py index 90a4bbe9e..e16b55bfa 100644 --- a/hydromt/raster.py +++ b/hydromt/raster.py @@ -18,6 +18,7 @@ import dask import geopandas as gpd +import mercantile as mct import numpy as np import pandas as pd import pyproj @@ -28,6 +29,8 @@ import xarray as xr import yaml from affine import Affine +from pathlib import Path +from PIL import Image from pyproj import CRS from rasterio import features from rasterio.enums import MergeAlg, Resampling @@ -37,7 +40,7 @@ from shapely.geometry import LineString, Polygon, box from . import _compat, gis_utils -from .utils import create_folder +from .utils import create_folder, elevation2rgb, rgb2elevation logger = logging.getLogger(__name__) XDIMS = ("x", "longitude", "lon", "long") @@ -200,6 +203,75 @@ def full_from_transform( return da +def tile_window_xyz(shape, px): + """Yield (left, upper, width, height).""" + nr, nc = shape + lu = product(range(0, nc, px), range(0, nr, px)) + + ## create the window + for l, u in lu: + h = min(px, nr - u) + w = min(px, nc - l) + yield (l, u, w, h) + + +def tile_window_osm( + zl, + y_ext, + px_size, + minx, + miny, + maxx, + maxy, +): + """Generate tiles in the xyz structure (according to OSM) + Returns: transform, bounds, x and y of the tile + """ + + # Calculate the dx, px parameters + dx = (y_ext * 2) / (2**zl) + px = dx / px_size + + # Origin displacement + odx = np.floor(abs(-y_ext - minx) / dx) + ody = np.floor(abs(y_ext - maxy) / dx) + + # Set the new origin + minx = -y_ext + odx * dx + maxy = y_ext - ody * dx + + # Create window generator + lu = product( + np.arange(round(minx, 2), round(maxx, 2), dx), + np.arange(round(maxy, 2), round(miny, 2), -dx), + ) + # Loop through the windows and return the transform, bounds, x and y number + for l, u in lu: + col = int(round(odx + (l - minx) / dx, 2)) + row = int(round(ody + (maxy - u) / dx, 2)) + bounds = [ + l, + u - dx, + l + dx, + u, + ] + yield Affine(px, 0, l, 0, -px, u), bounds, col, row + + +def tile_window_png(zl, minx, maxx, miny, maxy): + """Create tile windows according to OSM using mercantile + Returns: x, y and bounds of the tile + """ + ul = mct.tile(minx, maxy, zl) + lr = mct.tile(maxx, miny, zl) + + x_range = range(ul.x, lr.x + 1) + y_range = range(ul.y, lr.y + 1) + + for _x, _y in product(x_range, y_range): + yield _x, _y, mct.xy_bounds(_x, _y, zl) + + class XGeoBase(object): """Base class for the GIS extensions for xarray.""" @@ -2154,18 +2226,8 @@ def to_xyz_tiles( **kwargs Key-word arguments to write raster files """ - mName = os.path.normpath(os.path.basename(root)) - def tile_window(shape, px): - """Yield (left, upper, width, height).""" - nr, nc = shape - lu = product(range(0, nc, px), range(0, nr, px)) - - ## create the window - for l, u in lu: - h = min(px, nr - u) - w = min(px, nc - l) - yield (l, u, w, h) + mName = os.path.normpath(os.path.basename(root)) vrt_fn = None prev = 0 @@ -2196,7 +2258,7 @@ def tile_window(shape, px): txt_path = join(sd, "filelist.txt") file = open(txt_path, "w") - for l, u, w, h in tile_window(obj.shape, pxzl): + for l, u, w, h in tile_window_xyz(obj.shape, pxzl): col = int(np.ceil(l / pxzl)) row = int(np.ceil(u / pxzl)) ssd = join(sd, f"{col}") @@ -2278,112 +2340,125 @@ def to_osm_tiles( **kwargs Key-word arguments to write raster files """ + + # Ensure root directory exists and set name for datacatalog and stuff + create_folder(root) + m_name = os.path.normpath(os.path.basename(root)) + + # Set variables and flags + vrt_fn = None + + # Dimension count check before continuing + if self._obj.ndim != 2: + raise ValueError( + f"Only 2d datasets are accepted: {self._obj.ndim} dimensions are present" + ) + + # 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] + + # Fixed pixel size for XYZ tiles + px_size = 256 + + # Check whether the given resampling method is a valid one resampling = getattr(Resampling, method, None) if resampling is None: raise ValueError(f"Resampling method unknown: {method}.") - px_size = 256 # Fixed resolution for osm - assert self._obj.ndim == 2, "Only 2d datasets are accepted..." + # Object to local variable, also transpose it and extract some meta obj = self._obj.copy() obj = obj.transpose(self.y_dim, self.x_dim) + nodata = obj.raster.nodata - m_name = os.path.normpath(os.path.basename(root)) - + # Set bounds based on given bbox or extent of object if bbox: minx, miny, maxx, maxy = bbox else: minx, miny, maxx, maxy = obj.raster.transform_bounds("EPSG:3857") + # Make sure it lies within the extent of Pseudo Mercator minx, miny = map( max, - zip((minx, miny), [-20037508.34] * 2), + zip((minx, miny), [-y_ext_pm] * 2), ) maxx, maxy = map( min, - zip((maxx, maxy), [20037508.34] * 2), + zip((maxx, maxy), [y_ext_pm] * 2), ) + # Calculate the transform in 3857 for the given object + # (Mostly just for resolution) _obj_tr = rasterio.warp.calculate_default_transform( obj.raster.crs, "EPSG:3857", *obj.shape, *obj.raster.bounds, )[0] - dres = _obj_tr[0] - nzl = int(np.ceil((np.log10((20037508.34 * 2) / (dres * 256)) / np.log10(2)))) - - def tile_window(zl, minx, miny, maxx, maxy): - # Origin displacement - odx = np.floor(abs(-20037508.34 - minx) / dx) - ody = np.floor(abs(20037508.34 - maxy) / dx) - # Set the new origin - minx = -20037508.34 + odx * dx - maxy = 20037508.34 - ody * dx - - # Create window generator - lu = product( - np.arange(round(minx, 2), round(maxx, 2), dx), - np.arange(round(maxy, 2), round(miny, 2), -dx), - ) - for l, u in lu: - col = int(round(odx + (l - minx) / dx, 2)) - row = int(round(ody + (maxy - u) / dx, 2)) - bounds = [ - l, - u - dx, - l + dx, - u, - ] - yield Affine(px, 0, l, 0, -px, u), bounds, col, row + # Determine the max number of zoom levels with the resolution + dres = _obj_tr[0] + nzl = int(np.ceil((np.log10((y_ext_pm * 2) / (dres * px_size)) / np.log10(2)))) + # Set the zoom levels to be created zoom_levels.sort() if zoom_levels[-1] > min_lvl: min_lvl = zoom_levels[-1] zoom_levels = list(range(zoom_levels[0], min_lvl + 1, 1)) zoom_levels.reverse() - create_folder(root) + # Dict to save pixels size per zl, for the DataCatalog + zls_px_size = {} - vrt_fn = None - nodata = self.nodata - zls = {} + # Loop through the zoomlevels for zl in zoom_levels: - dx = (20037508.34 * 2) / (2**zl) + dx = (y_ext_pm * 2) / (2**zl) px = dx / px_size + # Simple warning for max zl exceedence if max(zoom_levels) > nzl: logger.warning( f"Pixels at zoomlevel {zl} are smaller than the original data" ) - # read data from previous zoomlevel + # read data from previous zoomlevel is there is one if vrt_fn is not None: obj = xr.open_dataarray(vrt_fn, engine="rasterio").squeeze( "band", drop=True ) + # Acquire some meta from the object + # Whether its the previous zl's vrt, or the original object obj_bounds = obj.raster.bounds obj_res = obj.raster.res[0] obj_shape = obj.raster.shape + # Get the transformer to warp coords to Pseudo Mercator t_r = pyproj.Transformer.from_crs( 3857, obj.raster.crs.to_epsg(), always_xy=True ) + # Ensure the directory is there per zoomlevel + # Also create the textfile needed for vrt creation sd = join(root, f"{zl}") create_folder(sd) txt_path = join(sd, "filelist.txt") file = open(txt_path, "w") - for transform, bounds, col, row in tile_window(zl, minx, miny, maxx, maxy): + # Loop through the windows + for transform, bounds, col, row in tile_window_osm( + zl, y_ext_pm, px_size, minx, miny, maxx, maxy + ): + # Ensure directory per x is there ssd = join(sd, f"{col}") create_folder(ssd) + # Transform bounding box of tile to 3857 ul = t_r.transform(bounds[0], bounds[3]) lr = t_r.transform(bounds[2], bounds[1]) out_bounds = [ul[0], lr[1], lr[0], ul[1]] + # Clip tile bounds to the bounds of the total object bounds temp_bounds = [ max(out_bounds[0], obj_bounds[0]), max(out_bounds[1], obj_bounds[1]), @@ -2391,17 +2466,18 @@ def tile_window(zl, minx, miny, maxx, maxy): min(out_bounds[3], obj_bounds[3]), ] + # This honestly is some weird stuff I wrote, but it works... + # Its needed because I index by number of rows and cols _window = [ abs(x - y) / obj_res for x, y in zip(temp_bounds, obj_bounds) ] - u = math.floor(round(_window[3], 2)) l = math.floor(round(_window[0], 2)) h = obj_shape[0] - u - math.floor(round(_window[1], 2)) w = obj_shape[1] - l - math.floor(round(_window[2], 2)) - temp = obj[u : u + h, l : l + w] + # Create empty tile tile = np.full((px_size, px_size), nodata, dtype=np.float64) temp_transform = Affine( @@ -2412,7 +2488,7 @@ def tile_window(zl, minx, miny, maxx, maxy): -obj_res, obj_bounds[3] - u * obj_res, ) - + # Warp the data and place it into the empty tile rasterio.warp.reproject( temp.values, tile, @@ -2423,7 +2499,7 @@ def tile_window(zl, minx, miny, maxx, maxy): dst_nodata=nodata, resampling=resampling, ) - + # Convert the tile (array) to a xr.DataArray in order to save properly temp = xr.DataArray( name=temp.name, dims=("y", "x"), @@ -2442,10 +2518,11 @@ def tile_window(zl, minx, miny, maxx, maxy): del tile + # Set some metadata temp.raster.set_crs(3857) - temp.raster.set_nodata(nodata) + # Save the data to the drive according to the driver if driver == "netcdf4": path = join(ssd, f"{row}.nc") temp = temp.raster.gdal_compliant() @@ -2465,7 +2542,7 @@ def tile_window(zl, minx, miny, maxx, maxy): # Create a vrt using GDAL vrt_fn = join(root, f"{m_name}_zl{zl}.vrt") gis_utils.create_vrt(vrt_fn, file_list_path=txt_path) - zls.update({zl: round(float(px), 2)}) + zls_px_size.update({zl: round(float(px), 2)}) del obj # Write a quick yaml for the database @@ -2474,11 +2551,253 @@ def tile_window(zl, minx, miny, maxx, maxy): "data_type": "RasterDataset", "driver": "raster", "path": f"{m_name}_zl{{zoom_level}}.vrt", - "zoom_levels": zls, + "zoom_levels": zls_px_size, } with open(join(root, f"{m_name}.yml"), "w") as f: yaml.dump({m_name: yml}, f, default_flow_style=False, sort_keys=False) + def to_webviewer_tiles( + self, + root: Path | str, + bbox: tuple | list = None, + method: str = "nearest", + ): + """Produce png's in XYZ structure (EPSG:3857) + Generally meant for webviewers + + Parameters + ---------- + root : Path | str + Path where the database will be saved + bbox : tuple | list, optional + A region where the tiles should be generated for + (if absent, the entire input raster dataset will be tiled) + method : str, optional + How to resample the data when downscaling. + E.g. 'nearest' for resampling with the nearest value + (This is only used for the first/ highest zoomlevel) + + Raises + ------ + ValueError + 2d arrays only + """ + + # Ensure the root directory is there + create_folder(root) + + # Set some flags and internal variables + _from_png = False + _count = 0 + _prev_zl = -1 + + # 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] + + # Fixed pixel size for XYZ tiles + px_size = 256 + + # Object to local variable, also transpose it and extract some meta + assert self._obj.ndim == 2, "Only 2d datasets are accepted..." + obj = self._obj.copy() + obj = obj.transpose(self.y_dim, self.x_dim) + obj_res = obj.raster.res[0] + obj_bounds = obj.raster.bounds + nodata = obj.raster.nodata + + # Calculate the transform and bounds in 3857 for the given object + # (Mostly just for the resolution) + # Bound is used for determining the lowest zoom level + tr_3857 = rasterio.warp.calculate_default_transform( + obj.raster.crs, + "EPSG:3857", + *obj.shape, + *obj.raster.bounds, + )[0] + bounds_3857 = rasterio.warp.transform_bounds( + obj.raster.crs, + "EPSG:3857", + *obj.raster.bounds, + ) + # bounds_3857 = [ + # min(max(_n, -y_ext_pm), y_ext_pm) for _n in bounds_3857 + # ] + + # Determine the max number of zoom levels with the resolution + dres = tr_3857[0] + nzl = int( + math.ceil((math.log10((y_ext_pm * 2) / (dres * px_size)) / math.log10(2))) + ) + + # Check whether the given resampling method is a valid one + resampling = getattr(Resampling, method, None) + if resampling is None: + raise ValueError(f"Resampling method unknown: {method}.") + + # Set bounds based on given bbox or extent of centroids from object + if bbox: + minx, miny, maxx, maxy = rasterio.warp.transform_bounds( + obj.crs, "EPSG:4326", *bbox + ) + else: + minx = obj_bounds[0] + 0.5 * obj_res + miny = obj_bounds[1] + 0.5 * obj_res + maxx = obj_bounds[2] - 0.5 * obj_res + maxy = obj_bounds[3] - 0.5 * obj_res + # For work to follow + minx, miny, maxx, maxy = rasterio.warp.transform_bounds( + obj.raster.crs, "EPSG:4326", minx, miny, maxx, maxy + ) + + # Make sure it lies within the extent of Pseudo Mercator + minx, miny = map( + max, + zip((minx, miny), (-180, -y_ext)), + ) + maxx, maxy = map( + min, + zip((maxx, maxy), (180, y_ext)), + ) + + # For determining the most course zoom level + bounds_3857 = rasterio.warp.transform_bounds( + "EPSG:4326", "EPSG:3857", minx, miny, maxx, maxy + ) + + # Determine the lowest zoom level (meaning only one tile, remains) + xsize = bounds_3857[2] - bounds_3857[0] + ysize = bounds_3857[3] - bounds_3857[1] + xlvl = math.floor(math.log((y_ext_pm * 2) / xsize) / math.log(2)) + ylvl = math.floor(math.log((y_ext_pm * 2) / ysize) / math.log(2)) + minl = int(min(xlvl, ylvl)) + + # Shift for all the tiles from the previous level + shift = ( + (0, 0), + (0, 1), + (1, 0), + (1, 1), + ) + + # Loop through the zoom levels + logger.info(f"Producing tiles from zoomlevel {nzl} to {minl}") + for zl in range(nzl, minl - 1, -1): + dst_res = (y_ext_pm * 2) / (px_size * 2**zl) # resolution per zoom level + # Create the sub directory for the zoom levels + sd = Path(root, f"{zl}") + create_folder(sd) + + # The highest zoomlevel has to be done from the original data + if not _from_png: + # Go through the zoomlevels + for x, y, bbox in tile_window_png(zl, minx, maxx, miny, maxy): + # Ensure the directory for every x location is there + ssd = Path(sd, f"{x}") + create_folder(ssd) + + # Get tile bounds an dst transfrom in 3857 + tile_bounds = rasterio.warp.transform_bounds( + "EPSG:3857", + obj.raster.crs, + *bbox, + ) + dst_transform = Affine(dst_res, 0, bbox[0], 0, -dst_res, bbox[3]) + + # Clip tile bounds with the boundary of the object + temp_bounds = [ + max(tile_bounds[0], obj_bounds[0]), + max(tile_bounds[1], obj_bounds[1]), + min(tile_bounds[2], obj_bounds[2]), + min(tile_bounds[3], obj_bounds[3]), + ] + # Select that data + z = obj.sel( + x=slice(temp_bounds[0] - obj_res, temp_bounds[2] + obj_res), + y=slice(temp_bounds[3] + obj_res, temp_bounds[1] - obj_res), + ) + z_bounds = z.raster.bounds + tile = np.full((px_size, px_size), np.nan, dtype=np.float64) + z_transform = Affine( + obj_res, + 0, + z_bounds[0], + 0, + -obj_res, + z_bounds[3], + ) + # Warp that data to Pseudo Mercator + rasterio.warp.reproject( + z.values, + tile, + src_transform=z_transform, + src_crs=obj.raster.crs, + dst_transform=dst_transform, + dst_crs="EPSG:3857", + dst_nodata=np.nan, + resampling=resampling, + ) + # Create RGB bands from the warped data and add transparency to nan values + rgb = elevation2rgb(tile.copy()) + _alpha = np.full((px_size, px_size), 255, dtype=np.uint8) + _alpha[np.where(np.isnan(tile))] = 0 + rgb = np.stack(rgb + (_alpha,), axis=2) + rgb_array = Image.fromarray(rgb) + rgb_array.save(Path(ssd, f"{y}.png")) + # Set the flag to True, as now png's can be used from here on out + _from_png = True + _prev_zl = zl + + # Use png's from the harddrive from the previous zoom level + else: + for x, y, _ in tile_window_png(zl, minx, maxx, miny, maxy): + # Ensure directory + ssd = Path(sd, f"{x}") + create_folder(ssd) + # Create a temporary array, 4 times the size of a tile + temp = np.full((px_size * 2, px_size * 2), np.nan, dtype=np.float64) + # Every tile from this level has 4 on the previous + # Go though the shift to acquire all of them + for sh in shift: + loc = (x * 2 + sh[0], y * 2 + sh[1]) + _im_file = Path( + root, f"{_prev_zl}", f"{loc[0]}", f"{loc[1]}.png" + ) + # Check if the file is really there, if not: it was not written + if not _im_file.exists(): + continue + # Open the file and load the data as rgba bands + _im = Image.open(_im_file) + data_raw = np.array(_im.resize((px_size, px_size))) + data = data_raw[:, :, :3] + _elev = rgb2elevation( + *np.split(data, data.shape[2], axis=2) + ).reshape(px_size, px_size) + # Use the alpha band to set data to nodata where the alpha is 0 + # (This means that these cells were also nodata the previous time around) + _elev[np.where(data_raw[:, :, 3] == 0)] = np.nan + del data_raw + # Add it to the temporary array + temp[ + px_size * sh[1] : px_size * (1 + sh[1]), + px_size * sh[0] : px_size * (1 + sh[0]), + ] = _elev + # Downscale the array to the size of a tile + temp = np.nanmean( + np.nanmean(temp.reshape((256, 2, 256, 2)), axis=3), axis=1 + ) + # Do the rgba funzies again + rgb = elevation2rgb(temp) + _alpha = np.full((px_size, px_size), 255, dtype=np.uint8) + _alpha[np.where(np.isnan(temp))] = 0 + del temp + rgb = np.stack(rgb + (_alpha,), axis=2) + rgb_array = Image.fromarray(rgb) + # Save it to the drive + rgb_array.save(Path(ssd, f"{y}.png")) + # On to the next zoomlevel + _prev_zl = zl + def to_raster( self, raster_path, diff --git a/hydromt/utils.py b/hydromt/utils.py index b133a29ba..11c618815 100644 --- a/hydromt/utils.py +++ b/hydromt/utils.py @@ -1,5 +1,6 @@ """Utility functions for hydromt that have no other home.""" import os +import numpy as np def partition_dictionaries(left, right): @@ -43,6 +44,22 @@ def create_folder(path): os.makedirs(path) +def elevation2rgb(val): + """Convert elevation to rgb tuple""" + val += 32768 + r = np.floor(val / 256).astype(np.uint8) + g = np.floor(val % 256).astype(np.uint8) + b = np.floor((val - np.floor(val)) * 256).astype(np.uint8) + + return (r, g, b) + + +def rgb2elevation(r, g, b): + """Convert rgb tuple to elevation""" + val = (r * 256 + g + b / 256) - 32768 + return val + + def _dict_pprint(d): import json diff --git a/tests/test_raster.py b/tests/test_raster.py index b421eba90..a82cb172f 100644 --- a/tests/test_raster.py +++ b/tests/test_raster.py @@ -526,6 +526,10 @@ def test_to_xyz_tiles(tmpdir, rioda_large): with open(os.path.join(path, "dummy_xyz", "2", "filelist.txt"), "r") as f: assert len(f.readlines()) == 1 + test_bounds = [2.13, -2.13, 3.2, -1.07] + _test_r = open_raster(os.path.join(path, "dummy_xyz", "0", "2", "1.tif")) + assert [round(_n, 2) for _n in _test_r.raster.bounds] == test_bounds + def test_to_osm_tiles(tmpdir, rioda_large): path = str(tmpdir) @@ -538,3 +542,22 @@ def test_to_osm_tiles(tmpdir, rioda_large): assert len(f.readlines()) == 1 with open(os.path.join(path, "dummy_osm", "8", "filelist.txt"), "r") as f: assert len(f.readlines()) == 12 + + test_bounds = [0.0, -313086.07, 313086.07, -0.0] + _test_r = open_raster(os.path.join(path, "dummy_osm", "7", "64", "64.tif")) + assert [round(_n, 2) for _n in _test_r.raster.bounds] == test_bounds + + +def test_to_png_tiles(tmpdir, rioda_large): + path = str(tmpdir) + rioda_large.raster.to_webviewer_tiles( + os.path.join(path, "dummy_png"), + ) + + _zl = os.listdir( + os.path.join(path, "dummy_png"), + ) + _zl = [int(_n) for _n in _zl] + assert len(_zl) == 4 + assert min(_zl) == 6 + assert max(_zl) == 9 From 47ff64d66d6f280bf1ffb068f048d82281e0bf6e Mon Sep 17 00:00:00 2001 From: Brendan Date: Thu, 14 Sep 2023 17:30:32 +0200 Subject: [PATCH 19/35] Fixed dependency; linting stuff --- hydromt/raster.py | 13 ++++--------- hydromt/utils.py | 5 +++-- pyproject.toml | 1 + 3 files changed, 8 insertions(+), 11 deletions(-) diff --git a/hydromt/raster.py b/hydromt/raster.py index e16b55bfa..d66b0651d 100644 --- a/hydromt/raster.py +++ b/hydromt/raster.py @@ -14,6 +14,7 @@ import tempfile from itertools import product from os.path import isdir, join +from pathlib import Path from typing import Any, Optional, Union import dask @@ -29,7 +30,6 @@ import xarray as xr import yaml from affine import Affine -from pathlib import Path from PIL import Image from pyproj import CRS from rasterio import features @@ -225,9 +225,8 @@ def tile_window_osm( maxy, ): """Generate tiles in the xyz structure (according to OSM) - Returns: transform, bounds, x and y of the tile + Returns: transform, bounds, x and y of the tile. """ - # Calculate the dx, px parameters dx = (y_ext * 2) / (2**zl) px = dx / px_size @@ -260,7 +259,7 @@ def tile_window_osm( def tile_window_png(zl, minx, maxx, miny, maxy): """Create tile windows according to OSM using mercantile - Returns: x, y and bounds of the tile + Returns: x, y and bounds of the tile. """ ul = mct.tile(minx, maxy, zl) lr = mct.tile(maxx, miny, zl) @@ -2226,7 +2225,6 @@ def to_xyz_tiles( **kwargs Key-word arguments to write raster files """ - mName = os.path.normpath(os.path.basename(root)) vrt_fn = None @@ -2340,7 +2338,6 @@ def to_osm_tiles( **kwargs Key-word arguments to write raster files """ - # Ensure root directory exists and set name for datacatalog and stuff create_folder(root) m_name = os.path.normpath(os.path.basename(root)) @@ -2563,7 +2560,7 @@ def to_webviewer_tiles( method: str = "nearest", ): """Produce png's in XYZ structure (EPSG:3857) - Generally meant for webviewers + Generally meant for webviewers. Parameters ---------- @@ -2582,7 +2579,6 @@ def to_webviewer_tiles( ValueError 2d arrays only """ - # Ensure the root directory is there create_folder(root) @@ -2604,7 +2600,6 @@ def to_webviewer_tiles( obj = obj.transpose(self.y_dim, self.x_dim) obj_res = obj.raster.res[0] obj_bounds = obj.raster.bounds - nodata = obj.raster.nodata # Calculate the transform and bounds in 3857 for the given object # (Mostly just for the resolution) diff --git a/hydromt/utils.py b/hydromt/utils.py index 15eb1cf6c..3537811ec 100644 --- a/hydromt/utils.py +++ b/hydromt/utils.py @@ -1,5 +1,6 @@ """Utility functions for hydromt that have no other home.""" import os + import numpy as np @@ -50,7 +51,7 @@ def create_folder(path): def elevation2rgb(val): - """Convert elevation to rgb tuple""" + """Convert elevation to rgb tuple.""" val += 32768 r = np.floor(val / 256).astype(np.uint8) g = np.floor(val % 256).astype(np.uint8) @@ -60,7 +61,7 @@ def elevation2rgb(val): def rgb2elevation(r, g, b): - """Convert rgb tuple to elevation""" + """Convert rgb tuple to elevation.""" val = (r * 256 + g + b / 256) - 32768 return val diff --git a/pyproject.toml b/pyproject.toml index f3dba74e3..5f87c8e24 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -17,6 +17,7 @@ dependencies = [ "entrypoints", # Provide access for Plugins "geopandas>=0.10", # pandas but geo, wraps fiona and shapely "gdal>=3.1", # enable geospacial data manipulation, both raster and victor + "mercantile", # for easilty tiling in XYZ structure (OSM) "numba", # speed up computations (used in e.g. stats) "numpy>=1.20", # pin necessary to ensure compatability with C headers "netcdf4", # netcfd IO From cbd1ff0e08ac704e968d4cf2a2a094d29b62d4b9 Mon Sep 17 00:00:00 2001 From: Brendan Date: Thu, 14 Sep 2023 17:38:19 +0200 Subject: [PATCH 20/35] linting..... --- hydromt/raster.py | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/hydromt/raster.py b/hydromt/raster.py index d66b0651d..9a57a86ab 100644 --- a/hydromt/raster.py +++ b/hydromt/raster.py @@ -224,7 +224,8 @@ def tile_window_osm( maxx, maxy, ): - """Generate tiles in the xyz structure (according to OSM) + """Generate tiles in the xyz structure (according to OSM). + Returns: transform, bounds, x and y of the tile. """ # Calculate the dx, px parameters @@ -258,7 +259,8 @@ def tile_window_osm( def tile_window_png(zl, minx, maxx, miny, maxy): - """Create tile windows according to OSM using mercantile + """Create tile windows according to OSM using mercantile. + Returns: x, y and bounds of the tile. """ ul = mct.tile(minx, maxy, zl) @@ -2348,7 +2350,8 @@ def to_osm_tiles( # Dimension count check before continuing if self._obj.ndim != 2: raise ValueError( - f"Only 2d datasets are accepted: {self._obj.ndim} dimensions are present" + f"Only 2d datasets are accepted: {self._obj.ndim} \ + dimensions are present" ) # Extent in y-direction for pseudo mercator (EPSG:3857) @@ -2559,7 +2562,8 @@ def to_webviewer_tiles( bbox: tuple | list = None, method: str = "nearest", ): - """Produce png's in XYZ structure (EPSG:3857) + """Produce png's in XYZ structure (EPSG:3857). + Generally meant for webviewers. Parameters @@ -2732,7 +2736,8 @@ def to_webviewer_tiles( dst_nodata=np.nan, resampling=resampling, ) - # Create RGB bands from the warped data and add transparency to nan values + # Create RGB bands from the warped data + # and add transparency to nan values rgb = elevation2rgb(tile.copy()) _alpha = np.full((px_size, px_size), 255, dtype=np.uint8) _alpha[np.where(np.isnan(tile))] = 0 @@ -2769,7 +2774,8 @@ def to_webviewer_tiles( *np.split(data, data.shape[2], axis=2) ).reshape(px_size, px_size) # Use the alpha band to set data to nodata where the alpha is 0 - # (This means that these cells were also nodata the previous time around) + # (This means that these cells were + # also nodata the previous time around) _elev[np.where(data_raw[:, :, 3] == 0)] = np.nan del data_raw # Add it to the temporary array From 5f9c17bedb5f1f8e0b08131152341036a318930b Mon Sep 17 00:00:00 2001 From: Brendan Date: Thu, 14 Sep 2023 17:39:19 +0200 Subject: [PATCH 21/35] remove whitespace --- hydromt/raster.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hydromt/raster.py b/hydromt/raster.py index 9a57a86ab..e7b1cbf8c 100644 --- a/hydromt/raster.py +++ b/hydromt/raster.py @@ -2774,7 +2774,7 @@ def to_webviewer_tiles( *np.split(data, data.shape[2], axis=2) ).reshape(px_size, px_size) # Use the alpha band to set data to nodata where the alpha is 0 - # (This means that these cells were + # (This means that these cells were # also nodata the previous time around) _elev[np.where(data_raw[:, :, 3] == 0)] = np.nan del data_raw From 45bf02efdd4e09d0b9d62b3ea189a02fecc5acfb Mon Sep 17 00:00:00 2001 From: Dirk Eilander Date: Wed, 20 Sep 2023 22:00:19 +0200 Subject: [PATCH 22/35] remove unnecessary create_folder method --- hydromt/raster.py | 29 +++++++++++++---------------- hydromt/utils.py | 7 ------- 2 files changed, 13 insertions(+), 23 deletions(-) diff --git a/hydromt/raster.py b/hydromt/raster.py index e7b1cbf8c..542176baf 100644 --- a/hydromt/raster.py +++ b/hydromt/raster.py @@ -13,7 +13,7 @@ import os import tempfile from itertools import product -from os.path import isdir, join +from os.path import join from pathlib import Path from typing import Any, Optional, Union @@ -40,7 +40,7 @@ from shapely.geometry import LineString, Polygon, box from . import _compat, gis_utils -from .utils import create_folder, elevation2rgb, rgb2elevation +from .utils import elevation2rgb, rgb2elevation logger = logging.getLogger(__name__) XDIMS = ("x", "longitude", "lon", "long") @@ -2254,7 +2254,7 @@ def to_xyz_tiles( # Write the raster paths to a text file sd = join(root, f"{zl}") - create_folder(sd) + os.makedirs(sd, exist_ok=True) txt_path = join(sd, "filelist.txt") file = open(txt_path, "w") @@ -2263,9 +2263,8 @@ def to_xyz_tiles( row = int(np.ceil(u / pxzl)) ssd = join(sd, f"{col}") - create_folder(ssd) - # create temp tile + os.makedirs(ssd, exist_ok=True) temp = obj[u : u + h, l : l + w] if zl != 0: temp = temp.coarsen( @@ -2341,7 +2340,7 @@ def to_osm_tiles( Key-word arguments to write raster files """ # Ensure root directory exists and set name for datacatalog and stuff - create_folder(root) + os.makedirs(root, exist_ok=True) m_name = os.path.normpath(os.path.basename(root)) # Set variables and flags @@ -2441,7 +2440,7 @@ def to_osm_tiles( # Ensure the directory is there per zoomlevel # Also create the textfile needed for vrt creation sd = join(root, f"{zl}") - create_folder(sd) + os.makedirs(sd, exist_ok=True) txt_path = join(sd, "filelist.txt") file = open(txt_path, "w") @@ -2451,7 +2450,7 @@ def to_osm_tiles( ): # Ensure directory per x is there ssd = join(sd, f"{col}") - create_folder(ssd) + os.makedirs(ssd, exist_ok=True) # Transform bounding box of tile to 3857 ul = t_r.transform(bounds[0], bounds[3]) @@ -2584,7 +2583,7 @@ def to_webviewer_tiles( 2d arrays only """ # Ensure the root directory is there - create_folder(root) + os.makedirs(root, exist_ok=True) # Set some flags and internal variables _from_png = False @@ -2685,7 +2684,7 @@ def to_webviewer_tiles( dst_res = (y_ext_pm * 2) / (px_size * 2**zl) # resolution per zoom level # Create the sub directory for the zoom levels sd = Path(root, f"{zl}") - create_folder(sd) + os.makedirs(sd, exist_ok=True) # The highest zoomlevel has to be done from the original data if not _from_png: @@ -2693,7 +2692,7 @@ def to_webviewer_tiles( for x, y, bbox in tile_window_png(zl, minx, maxx, miny, maxy): # Ensure the directory for every x location is there ssd = Path(sd, f"{x}") - create_folder(ssd) + os.makedirs(ssd, exist_ok=True) # Get tile bounds an dst transfrom in 3857 tile_bounds = rasterio.warp.transform_bounds( @@ -2753,7 +2752,7 @@ def to_webviewer_tiles( for x, y, _ in tile_window_png(zl, minx, maxx, miny, maxy): # Ensure directory ssd = Path(sd, f"{x}") - create_folder(ssd) + os.makedirs(ssd, exist_ok=True) # Create a temporary array, 4 times the size of a tile temp = np.full((px_size * 2, px_size * 2), np.nan, dtype=np.float64) # Every tile from this level has 4 on the previous @@ -3220,8 +3219,7 @@ def to_mapstack( if driver not in gis_utils.GDAL_EXT_CODE_MAP: raise ValueError(f"Extension unknown for driver: {driver}") ext = gis_utils.GDAL_EXT_CODE_MAP.get(driver) - if not isdir(root): - os.makedirs(root) + os.makedirs(root, exist_ok=True) with tempfile.TemporaryDirectory() as tmpdir: if driver == "PCRaster" and _compat.HAS_PCRASTER: clone_path = gis_utils.write_clone( @@ -3235,8 +3233,7 @@ def to_mapstack( if "/" in var: # variables with in subfolders folders = "/".join(var.split("/")[:-1]) - if not isdir(join(root, folders)): - os.makedirs(join(root, folders)) + os.makedirs(join(root, folders), exist_ok=True) var0 = var.split("/")[-1] raster_path = join(root, folders, f"{prefix}{var0}{postfix}.{ext}") else: diff --git a/hydromt/utils.py b/hydromt/utils.py index 3537811ec..5bb224753 100644 --- a/hydromt/utils.py +++ b/hydromt/utils.py @@ -1,5 +1,4 @@ """Utility functions for hydromt that have no other home.""" -import os import numpy as np @@ -44,12 +43,6 @@ def partition_dictionaries(left, right): return common, left_less_right, right_less_left -def create_folder(path): - """Create a folder if there is none.""" - if not os.path.exists(path): - os.makedirs(path) - - def elevation2rgb(val): """Convert elevation to rgb tuple.""" val += 32768 From f948149d01604bec32299a2017913dba7695068d Mon Sep 17 00:00:00 2001 From: Dirk Eilander Date: Thu, 21 Sep 2023 02:14:43 +0200 Subject: [PATCH 23/35] simplify slippy tiles --- hydromt/raster.py | 672 ++++++++++--------------------------------- hydromt/utils.py | 11 +- tests/test_raster.py | 38 ++- 3 files changed, 170 insertions(+), 551 deletions(-) diff --git a/hydromt/raster.py b/hydromt/raster.py index 542176baf..83aa646b2 100644 --- a/hydromt/raster.py +++ b/hydromt/raster.py @@ -7,7 +7,6 @@ """Extension for xarray to provide rasterio capabilities to xarray datasets/arrays.""" from __future__ import annotations -import gc import logging import math import os @@ -40,7 +39,7 @@ from shapely.geometry import LineString, Polygon, box from . import _compat, gis_utils -from .utils import elevation2rgb, rgb2elevation +from .utils import elevation2rgba, rgba2elevation logger = logging.getLogger(__name__) XDIMS = ("x", "longitude", "lon", "long") @@ -215,64 +214,6 @@ def tile_window_xyz(shape, px): yield (l, u, w, h) -def tile_window_osm( - zl, - y_ext, - px_size, - minx, - miny, - maxx, - maxy, -): - """Generate tiles in the xyz structure (according to OSM). - - Returns: transform, bounds, x and y of the tile. - """ - # Calculate the dx, px parameters - dx = (y_ext * 2) / (2**zl) - px = dx / px_size - - # Origin displacement - odx = np.floor(abs(-y_ext - minx) / dx) - ody = np.floor(abs(y_ext - maxy) / dx) - - # Set the new origin - minx = -y_ext + odx * dx - maxy = y_ext - ody * dx - - # Create window generator - lu = product( - np.arange(round(minx, 2), round(maxx, 2), dx), - np.arange(round(maxy, 2), round(miny, 2), -dx), - ) - # Loop through the windows and return the transform, bounds, x and y number - for l, u in lu: - col = int(round(odx + (l - minx) / dx, 2)) - row = int(round(ody + (maxy - u) / dx, 2)) - bounds = [ - l, - u - dx, - l + dx, - u, - ] - yield Affine(px, 0, l, 0, -px, u), bounds, col, row - - -def tile_window_png(zl, minx, maxx, miny, maxy): - """Create tile windows according to OSM using mercantile. - - Returns: x, y and bounds of the tile. - """ - ul = mct.tile(minx, maxy, zl) - lr = mct.tile(maxx, miny, zl) - - x_range = range(ul.x, lr.x + 1) - y_range = range(ul.y, lr.y + 1) - - for _x, _y in product(x_range, y_range): - yield _x, _y, mct.xy_bounds(_x, _y, zl) - - class XGeoBase(object): """Base class for the GIS extensions for xarray.""" @@ -2306,262 +2247,16 @@ def to_xyz_tiles( with open(join(root, f"{mName}.yml"), "w") as f: yaml.dump({mName: yml}, f, default_flow_style=False, sort_keys=False) - def to_osm_tiles( - self, - root: str, - zoom_levels: list, - bbox: tuple = (), - method: str = "nearest", - min_lvl: int = 5, - driver="GTiff", - **kwargs, - ): - """Export rasterdataset to tiles in an osm structure. - - Parameters - ---------- - root : str - Path where the database will be saved - Database yml will be put one directory above - zoom_levels : list - Zoom levels to be put in the database - bbox : tuple, optional - A region where the tiles should be generated for - (if absent, the entire input raster dataset will be tiled) - method : str, optional - How to resample the data when downscaling. - E.g. 'nearest' for resampling with the nearest value - min_lvl : int, optional - Minimal starting level from to where to downscale from. - Increase if memory issues are encountered. - driver : str, optional - GDAL driver (e.g., 'GTiff' for geotif files), or 'netcdf4' for netcdf files. - **kwargs - Key-word arguments to write raster files - """ - # Ensure root directory exists and set name for datacatalog and stuff - os.makedirs(root, exist_ok=True) - m_name = os.path.normpath(os.path.basename(root)) - - # Set variables and flags - vrt_fn = None - - # Dimension count check before continuing - if self._obj.ndim != 2: - raise ValueError( - f"Only 2d datasets are accepted: {self._obj.ndim} \ - dimensions are present" - ) - - # 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] - - # Fixed pixel size for XYZ tiles - px_size = 256 - - # Check whether the given resampling method is a valid one - resampling = getattr(Resampling, method, None) - if resampling is None: - raise ValueError(f"Resampling method unknown: {method}.") - - # Object to local variable, also transpose it and extract some meta - obj = self._obj.copy() - obj = obj.transpose(self.y_dim, self.x_dim) - nodata = obj.raster.nodata - - # Set bounds based on given bbox or extent of object - if bbox: - minx, miny, maxx, maxy = bbox - else: - minx, miny, maxx, maxy = obj.raster.transform_bounds("EPSG:3857") - - # Make sure it lies within the extent of Pseudo Mercator - minx, miny = map( - max, - zip((minx, miny), [-y_ext_pm] * 2), - ) - maxx, maxy = map( - min, - zip((maxx, maxy), [y_ext_pm] * 2), - ) - - # Calculate the transform in 3857 for the given object - # (Mostly just for resolution) - _obj_tr = rasterio.warp.calculate_default_transform( - obj.raster.crs, - "EPSG:3857", - *obj.shape, - *obj.raster.bounds, - )[0] - - # Determine the max number of zoom levels with the resolution - dres = _obj_tr[0] - nzl = int(np.ceil((np.log10((y_ext_pm * 2) / (dres * px_size)) / np.log10(2)))) - - # Set the zoom levels to be created - zoom_levels.sort() - if zoom_levels[-1] > min_lvl: - min_lvl = zoom_levels[-1] - zoom_levels = list(range(zoom_levels[0], min_lvl + 1, 1)) - zoom_levels.reverse() - - # Dict to save pixels size per zl, for the DataCatalog - zls_px_size = {} - - # Loop through the zoomlevels - for zl in zoom_levels: - dx = (y_ext_pm * 2) / (2**zl) - px = dx / px_size - - # Simple warning for max zl exceedence - if max(zoom_levels) > nzl: - logger.warning( - f"Pixels at zoomlevel {zl} are smaller than the original data" - ) - - # read data from previous zoomlevel is there is one - if vrt_fn is not None: - obj = xr.open_dataarray(vrt_fn, engine="rasterio").squeeze( - "band", drop=True - ) - - # Acquire some meta from the object - # Whether its the previous zl's vrt, or the original object - obj_bounds = obj.raster.bounds - obj_res = obj.raster.res[0] - obj_shape = obj.raster.shape - - # Get the transformer to warp coords to Pseudo Mercator - t_r = pyproj.Transformer.from_crs( - 3857, obj.raster.crs.to_epsg(), always_xy=True - ) - - # Ensure the directory is there per zoomlevel - # Also create the textfile needed for vrt creation - sd = join(root, f"{zl}") - os.makedirs(sd, exist_ok=True) - txt_path = join(sd, "filelist.txt") - file = open(txt_path, "w") - - # Loop through the windows - for transform, bounds, col, row in tile_window_osm( - zl, y_ext_pm, px_size, minx, miny, maxx, maxy - ): - # Ensure directory per x is there - ssd = join(sd, f"{col}") - os.makedirs(ssd, exist_ok=True) - - # Transform bounding box of tile to 3857 - ul = t_r.transform(bounds[0], bounds[3]) - lr = t_r.transform(bounds[2], bounds[1]) - out_bounds = [ul[0], lr[1], lr[0], ul[1]] - - # Clip tile bounds to the bounds of the total object bounds - temp_bounds = [ - max(out_bounds[0], obj_bounds[0]), - max(out_bounds[1], obj_bounds[1]), - min(out_bounds[2], obj_bounds[2]), - min(out_bounds[3], obj_bounds[3]), - ] - - # This honestly is some weird stuff I wrote, but it works... - # Its needed because I index by number of rows and cols - _window = [ - abs(x - y) / obj_res for x, y in zip(temp_bounds, obj_bounds) - ] - u = math.floor(round(_window[3], 2)) - l = math.floor(round(_window[0], 2)) - h = obj_shape[0] - u - math.floor(round(_window[1], 2)) - w = obj_shape[1] - l - math.floor(round(_window[2], 2)) - temp = obj[u : u + h, l : l + w] - - # Create empty tile - tile = np.full((px_size, px_size), nodata, dtype=np.float64) - - temp_transform = Affine( - obj_res, - 0, - obj_bounds[0] + l * obj_res, - 0, - -obj_res, - obj_bounds[3] - u * obj_res, - ) - # Warp the data and place it into the empty tile - rasterio.warp.reproject( - temp.values, - tile, - src_transform=temp_transform, - src_crs=obj.raster.crs, - dst_transform=transform, - dst_crs="EPSG:3857", - dst_nodata=nodata, - resampling=resampling, - ) - # Convert the tile (array) to a xr.DataArray in order to save properly - temp = xr.DataArray( - name=temp.name, - dims=("y", "x"), - coords={ - "y": np.arange( - bounds[3] + 0.5 * transform[4], bounds[1], transform[4] - )[:px_size], - "x": np.arange( - bounds[0] + 0.5 * transform[0], - bounds[2] + 0.5 * transform[0], - transform[0], - )[:px_size], - }, - data=tile.copy(), - ) - - del tile - - # Set some metadata - temp.raster.set_crs(3857) - temp.raster.set_nodata(nodata) - - # Save the data to the drive according to the driver - if driver == "netcdf4": - path = join(ssd, f"{row}.nc") - temp = temp.raster.gdal_compliant() - temp.to_netcdf(path, engine="netcdf4", **kwargs) - elif driver in gis_utils.GDAL_EXT_CODE_MAP: - ext = gis_utils.GDAL_EXT_CODE_MAP.get(driver) - path = join(ssd, f"{row}.{ext}") - temp.raster.to_raster(path, driver=driver, **kwargs) - else: - raise ValueError(f"Unkown file driver {driver}") - file.write(f"{path}\n") - - del temp - gc.collect() - - file.close() - # Create a vrt using GDAL - vrt_fn = join(root, f"{m_name}_zl{zl}.vrt") - gis_utils.create_vrt(vrt_fn, file_list_path=txt_path) - zls_px_size.update({zl: round(float(px), 2)}) - del obj - - # Write a quick yaml for the database - yml = { - "crs": 3857, - "data_type": "RasterDataset", - "driver": "raster", - "path": f"{m_name}_zl{{zoom_level}}.vrt", - "zoom_levels": zls_px_size, - } - with open(join(root, f"{m_name}.yml"), "w") as f: - yaml.dump({m_name: yml}, f, default_flow_style=False, sort_keys=False) - - def to_webviewer_tiles( + def to_slippy_tiles( self, root: Path | str, - bbox: tuple | list = None, - method: str = "nearest", + reproj_method: str = "nearest", + min_lvl: int = None, + max_lvl: int = None, + driver="png", + **kwargs, ): - """Produce png's in XYZ structure (EPSG:3857). + """Produce tiles in /zoom/x/y. structure (EPSG:3857). Generally meant for webviewers. @@ -2569,234 +2264,159 @@ def to_webviewer_tiles( ---------- root : Path | str Path where the database will be saved - bbox : tuple | list, optional - A region where the tiles should be generated for - (if absent, the entire input raster dataset will be tiled) - method : str, optional + reproj_method : str, optional How to resample the data when downscaling. E.g. 'nearest' for resampling with the nearest value (This is only used for the first/ highest zoomlevel) - - Raises - ------ - ValueError - 2d arrays only + min_lvl, max_lvl : int, optional + The minimum and maximum zoomlevel to be produced. + If None, the zoomlevels will be determined based on the data resolution + driver : str, optional + file output driver, one of 'png', 'netcdf4' or 'GTiff' + **kwargs + Key-word arguments to write file + for netcdf4, these are passed to ~:py:meth:xarray.DataArray.to_netcdf: + for GTiff, these are passed to ~:py:meth:hydromt.RasterDataArray.to_raster: + for png, these are passed to ~:py:meth:PIL.Image.Image.save: """ - # Ensure the root directory is there - os.makedirs(root, exist_ok=True) - - # Set some flags and internal variables - _from_png = False - _count = 0 - _prev_zl = -1 - - # 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] - - # Fixed pixel size for XYZ tiles + # TODO add option to use cmap instead of elevation2rgba + # Fixed pixel size and CRS for XYZ tiles px_size = 256 + crs = CRS.from_epsg(3857) + ext = {"png": "png", "netcdf4": "nc", "gtiff": "tif"}.get(driver.lower(), None) + if ext is None: + raise ValueError(f"Unkown file driver {driver}, use png, netcdf4 or GTiff") + # for now we assume output has float32 dtype + kwargs0 = { + "netcdf4": {"engine": "netcdf4"}, + "gtiff": {"driver": "GTiff", "compress": "deflate", "dtype": "float32"}, + }.get(driver.lower(), {}) + kwargs = {**kwargs0, **kwargs} # Object to local variable, also transpose it and extract some meta - assert self._obj.ndim == 2, "Only 2d datasets are accepted..." - obj = self._obj.copy() - obj = obj.transpose(self.y_dim, self.x_dim) - obj_res = obj.raster.res[0] - obj_bounds = obj.raster.bounds - - # Calculate the transform and bounds in 3857 for the given object - # (Mostly just for the resolution) - # Bound is used for determining the lowest zoom level - tr_3857 = rasterio.warp.calculate_default_transform( - obj.raster.crs, - "EPSG:3857", - *obj.shape, - *obj.raster.bounds, - )[0] - bounds_3857 = rasterio.warp.transform_bounds( - obj.raster.crs, - "EPSG:3857", - *obj.raster.bounds, - ) - # bounds_3857 = [ - # min(max(_n, -y_ext_pm), y_ext_pm) for _n in bounds_3857 - # ] - - # Determine the max number of zoom levels with the resolution - dres = tr_3857[0] - nzl = int( - math.ceil((math.log10((y_ext_pm * 2) / (dres * px_size)) / math.log10(2))) - ) - - # Check whether the given resampling method is a valid one - resampling = getattr(Resampling, method, None) - if resampling is None: - raise ValueError(f"Resampling method unknown: {method}.") - - # Set bounds based on given bbox or extent of centroids from object - if bbox: - minx, miny, maxx, maxy = rasterio.warp.transform_bounds( - obj.crs, "EPSG:4326", *bbox - ) - else: - minx = obj_bounds[0] + 0.5 * obj_res - miny = obj_bounds[1] + 0.5 * obj_res - maxx = obj_bounds[2] - 0.5 * obj_res - maxy = obj_bounds[3] - 0.5 * obj_res - # For work to follow - minx, miny, maxx, maxy = rasterio.warp.transform_bounds( - obj.raster.crs, "EPSG:4326", minx, miny, maxx, maxy + if not self._obj.ndim == 2: + ValueError("Only 2d datasets are accepted...") + # make sure the data is N-S oriented, with y-axis as first dimension + # and nodata values set to nan + obj = self.mask_nodata().transpose(self.y_dim, self.x_dim) + if obj.raster.res[1] < 0: + obj = obj.raster.flipud() + + # Calculate min/max zoomlevel based on the resolution of the data + bounds_wgsg84 = self.transform_bounds("EPSG:4326") + if min_lvl is None: + min_lat = min(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(min_lat, *obj.raster.res)[0] + C = 2 * np.pi * 6378137 # circumference of the earth + min_lvl = int( + np.ceil(np.log2(C * np.cos(np.deg2rad(min_lat)) / dx_m / px_size)) ) - - # Make sure it lies within the extent of Pseudo Mercator - minx, miny = map( - max, - zip((minx, miny), (-180, -y_ext)), - ) - maxx, maxy = map( - min, - zip((maxx, maxy), (180, y_ext)), - ) - - # For determining the most course zoom level - bounds_3857 = rasterio.warp.transform_bounds( - "EPSG:4326", "EPSG:3857", minx, miny, maxx, maxy - ) - - # Determine the lowest zoom level (meaning only one tile, remains) - xsize = bounds_3857[2] - bounds_3857[0] - ysize = bounds_3857[3] - bounds_3857[1] - xlvl = math.floor(math.log((y_ext_pm * 2) / xsize) / math.log(2)) - ylvl = math.floor(math.log((y_ext_pm * 2) / ysize) / math.log(2)) - minl = int(min(xlvl, ylvl)) - - # Shift for all the tiles from the previous level - shift = ( - (0, 0), - (0, 1), - (1, 0), - (1, 1), - ) + if max_lvl is None: + max_lvl = mct.bounding_tile(*bounds_wgsg84).z # Loop through the zoom levels - logger.info(f"Producing tiles from zoomlevel {nzl} to {minl}") - for zl in range(nzl, minl - 1, -1): - dst_res = (y_ext_pm * 2) / (px_size * 2**zl) # resolution per zoom level - # Create the sub directory for the zoom levels - sd = Path(root, f"{zl}") - os.makedirs(sd, exist_ok=True) - - # The highest zoomlevel has to be done from the original data - if not _from_png: - # Go through the zoomlevels - for x, y, bbox in tile_window_png(zl, minx, maxx, miny, maxy): - # Ensure the directory for every x location is there - ssd = Path(sd, f"{x}") + zoom_levels = {} + logger.info(f"Producing tiles from zoomlevel {min_lvl} to {max_lvl}") + for zl in range(min_lvl, max_lvl - 1, -1): + fns = [] + # Go through the zoomlevels + for i, tile in enumerate(mct.tiles(*bounds_wgsg84, zl, truncate=True)): + if i == 0: + ssd = Path(root, str(zl), f"{tile.x}") os.makedirs(ssd, exist_ok=True) - - # Get tile bounds an dst transfrom in 3857 - tile_bounds = rasterio.warp.transform_bounds( - "EPSG:3857", - obj.raster.crs, - *bbox, + bounds0 = mct.xy_bounds(tile) + zoom_levels[zl] = abs(bounds0[2] - bounds0[0]) / px_size + if zl == min_lvl: + # For the first zoomlevel, we can just clip the data + # does this need a try/except? + src_tile = obj.raster.clip_bbox( + mct.xy_bounds(tile), crs=crs, buffer=2, align=True ) - dst_transform = Affine(dst_res, 0, bbox[0], 0, -dst_res, bbox[3]) - - # Clip tile bounds with the boundary of the object - temp_bounds = [ - max(tile_bounds[0], obj_bounds[0]), - max(tile_bounds[1], obj_bounds[1]), - min(tile_bounds[2], obj_bounds[2]), - min(tile_bounds[3], obj_bounds[3]), - ] - # Select that data - z = obj.sel( - x=slice(temp_bounds[0] - obj_res, temp_bounds[2] + obj_res), - y=slice(temp_bounds[3] + obj_res, temp_bounds[1] - obj_res), - ) - z_bounds = z.raster.bounds - tile = np.full((px_size, px_size), np.nan, dtype=np.float64) - z_transform = Affine( - obj_res, - 0, - z_bounds[0], - 0, - -obj_res, - z_bounds[3], - ) - # Warp that data to Pseudo Mercator - rasterio.warp.reproject( - z.values, - tile, - src_transform=z_transform, - src_crs=obj.raster.crs, - dst_transform=dst_transform, - dst_crs="EPSG:3857", - dst_nodata=np.nan, - resampling=resampling, - ) - # Create RGB bands from the warped data - # and add transparency to nan values - rgb = elevation2rgb(tile.copy()) - _alpha = np.full((px_size, px_size), 255, dtype=np.uint8) - _alpha[np.where(np.isnan(tile))] = 0 - rgb = np.stack(rgb + (_alpha,), axis=2) - rgb_array = Image.fromarray(rgb) - rgb_array.save(Path(ssd, f"{y}.png")) - # Set the flag to True, as now png's can be used from here on out - _from_png = True - _prev_zl = zl - - # Use png's from the harddrive from the previous zoom level - else: - for x, y, _ in tile_window_png(zl, minx, maxx, miny, maxy): - # Ensure directory - ssd = Path(sd, f"{x}") - os.makedirs(ssd, exist_ok=True) + else: + # Every tile from this level has 4 child tiles on the previous lvl + children = [] + for child in mct.children(tile): + fn = Path(root, str(child.z), str(child.x), f"{child.y}.{ext}") + # Check if the file is really there, if not: it was not written + if fn.exists(): + if driver == "netcdf4": + data = xr.open_dataarray(fn).values + elif driver == "GTiff": + data = xr.open_dataarray( + fn, decode_coords=False, engine="rasterio" + ) + data = data.squeeze(drop=True).values + elif driver == "png": + data = rgba2elevation(np.array(Image.open(fn))) + else: + data = np.full((px_size, px_size), np.nan, dtype=np.float64) + # Add the data to the temporary array + children.append(data) # Create a temporary array, 4 times the size of a tile + # order of children: top-left, top-right, bottom-right, bottom-left temp = np.full((px_size * 2, px_size * 2), np.nan, dtype=np.float64) - # Every tile from this level has 4 on the previous - # Go though the shift to acquire all of them - for sh in shift: - loc = (x * 2 + sh[0], y * 2 + sh[1]) - _im_file = Path( - root, f"{_prev_zl}", f"{loc[0]}", f"{loc[1]}.png" - ) - # Check if the file is really there, if not: it was not written - if not _im_file.exists(): - continue - # Open the file and load the data as rgba bands - _im = Image.open(_im_file) - data_raw = np.array(_im.resize((px_size, px_size))) - data = data_raw[:, :, :3] - _elev = rgb2elevation( - *np.split(data, data.shape[2], axis=2) - ).reshape(px_size, px_size) - # Use the alpha band to set data to nodata where the alpha is 0 - # (This means that these cells were - # also nodata the previous time around) - _elev[np.where(data_raw[:, :, 3] == 0)] = np.nan - del data_raw - # Add it to the temporary array - temp[ - px_size * sh[1] : px_size * (1 + sh[1]), - px_size * sh[0] : px_size * (1 + sh[0]), - ] = _elev - # Downscale the array to the size of a tile - temp = np.nanmean( - np.nanmean(temp.reshape((256, 2, 256, 2)), axis=3), axis=1 + temp[:px_size, :px_size] = children[0] + temp[:px_size:, px_size:] = children[1] + temp[px_size:, px_size:] = children[2] + temp[px_size:, :px_size] = children[3] + # create a dataarray from the temporary array + src_transform = rasterio.transform.from_bounds( + *mct.xy_bounds(tile), px_size, px_size + ) + src_tile = RasterDataArray.from_numpy( + temp, src_transform, crs=crs, nodata=np.nan ) - # Do the rgba funzies again - rgb = elevation2rgb(temp) - _alpha = np.full((px_size, px_size), 255, dtype=np.uint8) - _alpha[np.where(np.isnan(temp))] = 0 - del temp - rgb = np.stack(rgb + (_alpha,), axis=2) - rgb_array = Image.fromarray(rgb) - # Save it to the drive - rgb_array.save(Path(ssd, f"{y}.png")) - # On to the next zoomlevel - _prev_zl = zl + + # reproject the data to the tile + dst_transform = rasterio.transform.from_bounds( + *mct.xy_bounds(tile), px_size, px_size + ) + dst_tile = src_tile.raster.reproject( + dst_crs=crs, + dst_transform=dst_transform, + dst_width=px_size, + dst_height=px_size, + method=reproj_method, + ) + # write the data to file + fn_out = Path(ssd, f"{tile.y}.{ext}") + fns.append(fn_out) + if driver.lower() == "png": + # Create RGBA bands from the data and save it as png + rgba = elevation2rgba(dst_tile.values) + Image.fromarray(rgba).save(fn_out, **kwargs) + elif driver.lower() == "netcdf4": + # write the data to netcdf + dst_tile = dst_tile.raster.gdal_compliant() + dst_tile.to_netcdf(fn_out, **kwargs) + elif driver.lower() == "gtiff": + # write the data to geotiff + dst_tile.raster.to_raster(fn_out, **kwargs) + if driver.lower() != "png": + # Write files to txt and create a vrt using GDAL + txt_fn = Path(root, f"filelist_lvl{zl}.txt") + vrt_fn = Path(root, f"lvl{zl}.vrt") + with open(txt_fn, "w") as f: + for fn in fns: + f.write(f"{fn}\n") + gis_utils.create_vrt(vrt_fn, file_list_path=txt_fn) + + if driver.lower() != "png": + # Write a quick yaml for the database + # zoom level : resolution in meters + yml = { + "crs": 3857, + "data_type": "RasterDataset", + "driver": "raster", + "path": "lvl{zoom_level}.vrt", + "zoom_levels": zoom_levels, + } + name = os.path.basename(root) + with open(join(root, "data.yml"), "w") as f: + yaml.dump({name: yml}, f, default_flow_style=False, sort_keys=False) def to_raster( self, diff --git a/hydromt/utils.py b/hydromt/utils.py index 5bb224753..06739d5ae 100644 --- a/hydromt/utils.py +++ b/hydromt/utils.py @@ -43,20 +43,21 @@ def partition_dictionaries(left, right): return common, left_less_right, right_less_left -def elevation2rgb(val): +def elevation2rgba(val): """Convert elevation to rgb tuple.""" val += 32768 r = np.floor(val / 256).astype(np.uint8) g = np.floor(val % 256).astype(np.uint8) b = np.floor((val - np.floor(val)) * 256).astype(np.uint8) + a = np.where(np.isnan(val), 0, 255).astype(np.uint8) + return np.stack((r, g, b, a), axis=2) - return (r, g, b) - -def rgb2elevation(r, g, b): +def rgba2elevation(rgba: np.ndarray): """Convert rgb tuple to elevation.""" + r, g, b, a = np.split(rgba, 4, axis=2) val = (r * 256 + g + b / 256) - 32768 - return val + return np.where(a == 0, np.nan, val).squeeze() def _dict_pprint(d): diff --git a/tests/test_raster.py b/tests/test_raster.py index a82cb172f..d5c22a164 100644 --- a/tests/test_raster.py +++ b/tests/test_raster.py @@ -531,31 +531,29 @@ def test_to_xyz_tiles(tmpdir, rioda_large): assert [round(_n, 2) for _n in _test_r.raster.bounds] == test_bounds -def test_to_osm_tiles(tmpdir, rioda_large): - path = str(tmpdir) - rioda_large.raster.to_osm_tiles( - os.path.join(path, "dummy_osm"), - zoom_levels=[5], - min_lvl=8, - ) - with open(os.path.join(path, "dummy_osm", "5", "filelist.txt"), "r") as f: - assert len(f.readlines()) == 1 - with open(os.path.join(path, "dummy_osm", "8", "filelist.txt"), "r") as f: - assert len(f.readlines()) == 12 - - test_bounds = [0.0, -313086.07, 313086.07, -0.0] - _test_r = open_raster(os.path.join(path, "dummy_osm", "7", "64", "64.tif")) - assert [round(_n, 2) for _n in _test_r.raster.bounds] == test_bounds +# def test_to_osm_tiles(tmpdir, rioda_large): +# path = str(tmpdir) +# rioda_large.raster.to_osm_tiles( +# os.path.join(path, "dummy_osm"), +# zoom_levels=[5], +# min_lvl=8, +# ) +# with open(os.path.join(path, "dummy_osm", "5", "filelist.txt"), "r") as f: +# assert len(f.readlines()) == 1 +# with open(os.path.join(path, "dummy_osm", "8", "filelist.txt"), "r") as f: +# assert len(f.readlines()) == 12 + +# test_bounds = [0.0, -313086.07, 313086.07, -0.0] +# _test_r = open_raster(os.path.join(path, "dummy_osm", "7", "64", "64.tif")) +# assert [round(_n, 2) for _n in _test_r.raster.bounds] == test_bounds def test_to_png_tiles(tmpdir, rioda_large): - path = str(tmpdir) - rioda_large.raster.to_webviewer_tiles( - os.path.join(path, "dummy_png"), + rioda_large.raster.to_slippy_tiles( + os.path.join(tmpdir, "dummy_png"), ) - _zl = os.listdir( - os.path.join(path, "dummy_png"), + os.path.join(tmpdir, "dummy_png"), ) _zl = [int(_n) for _n in _zl] assert len(_zl) == 4 From e865cc0a8f42efc654b010f5421a7efc80bf93ba Mon Sep 17 00:00:00 2001 From: Dirk Eilander Date: Thu, 21 Sep 2023 09:30:34 +0200 Subject: [PATCH 24/35] add optional dependencies --- docs/api.rst | 3 +- examples/working_with_tiled_raster_data.ipynb | 45 +++--- hydromt/raster.py | 143 +++++++++++------- pyproject.toml | 7 +- 4 files changed, 115 insertions(+), 83 deletions(-) diff --git a/docs/api.rst b/docs/api.rst index 7343f55b5..ace34f702 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -480,8 +480,7 @@ Raster writing methods :toctree: _generated DataArray.raster.to_xyz_tiles - DataArray.raster.to_osm_tiles - DataArray.raster.to_webviewer_tiles + DataArray.raster.to_slippy_tiles DataArray.raster.to_raster Dataset.raster.to_mapstack diff --git a/examples/working_with_tiled_raster_data.ipynb b/examples/working_with_tiled_raster_data.ipynb index be41b6a52..f7cff62e1 100644 --- a/examples/working_with_tiled_raster_data.ipynb +++ b/examples/working_with_tiled_raster_data.ipynb @@ -27,7 +27,7 @@ "from hydromt.log import setuplog\n", "from os.path import join\n", "\n", - "logger = setuplog('tiling', log_level=20)\n", + "logger = setuplog(\"tiling\", log_level=20)\n", "\n", "# get some elevation data from the data catalog\n", "data_lib = \"artifact_data\"\n", @@ -99,7 +99,7 @@ "### Tiling raster with OSM stucture\n", "\n", "Now let's have a look at tiling according to the OSM structure\n", - "an xarray.DataArray is simple written to a tile database in OSM structure via .raster.to_osm_tiles" + "an xarray.DataArray is simple written to a tile database in OSM structure via ``.raster.to_slippy_tiles``" ] }, { @@ -111,12 +111,8 @@ "# Write the database in XYZ stucture\n", "name_osm = f\"{source}_osm\"\n", "root_osm = join(\"tmpdir\", name_osm)\n", - "zoom_levels = [0, 1, 2, 5] # note that you have to present all levels\n", - "da0.raster.to_osm_tiles(\n", - " root=root_osm,\n", - " zoom_levels=zoom_levels,\n", - " min_lvl=9, # minimal level from where to upscale from (prevent memory issues)\n", - " driver=\"GTiff\", # try also 'netcdf4'\n", + "da0.raster.to_slippy_tiles(\n", + " root=root_osm, driver=\"GTiff\", reproj_method=\"average\" # try also 'netcdf4'\n", ")" ] }, @@ -152,10 +148,15 @@ "outputs": [], "source": [ "# Write the data in OSM structure, but to images!\n", - "name_png = f\"{source}_png\"\n", + "from matplotlib import pyplot as plt\n", + "\n", + "\n", + "name_png = f\"{source}_png_cmap\"\n", "root_png = join(\"tmpdir\", name_png)\n", - "da0.raster.to_webviewer_tiles(\n", - " root = root_png\n", + "da0.raster.to_slippy_tiles(\n", + " root=root_png,\n", + " cmap=\"terrain\",\n", + " norm=plt.Normalize(vmin=0, vmax=2000),\n", ")" ] }, @@ -178,10 +179,10 @@ "\n", "# Create a figure to show the image\n", "fig = plt.figure()\n", - "fig.set_size_inches(2.5,2.5)\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", + "ax.set_position([0, 0, 1, 1])\n", + "ax.axis(\"off\")\n", "\n", "# Show the image\n", "im = Image.open(join(root_png, \"9\", \"273\", \"182.png\"))\n", @@ -221,7 +222,9 @@ "from hydromt import DataCatalog\n", "\n", "# Load the yml into a DataCatalog\n", - "data_catalog = DataCatalog([join(root, f\"{name}.yml\"), join(root_osm, f\"{name_osm}.yml\")], logger=logger)\n", + "data_catalog = DataCatalog(\n", + " [join(root, f\"{name}.yml\"), join(root_osm, f\"{name_osm}.yml\")], logger=logger\n", + ")\n", "\n", "# View the structure of the DataCatalog\n", "data_catalog[name]" @@ -266,11 +269,11 @@ "outputs": [], "source": [ "# Request a raster from the Datacatolog based on zoom resolution & unit\n", - "da = data_catalog.get_rasterdataset(name, zoom_level=(1/600, 'degree'))\n", - "da = data_catalog.get_rasterdataset(name_osm, zoom_level=(1/600, 'degree'))\n", + "da = data_catalog.get_rasterdataset(name, zoom_level=(1 / 600, \"degree\"))\n", + "da = data_catalog.get_rasterdataset(name_osm, zoom_level=(1 / 600, \"degree\"))\n", "\n", - "da = data_catalog.get_rasterdataset(name, zoom_level=(1e3, 'degree'))\n", - "da = data_catalog.get_rasterdataset(name_osm, zoom_level=(1e4, 'meter'))\n" + "da = data_catalog.get_rasterdataset(name, zoom_level=(1e3, \"degree\"))\n", + "da = data_catalog.get_rasterdataset(name_osm, zoom_level=(1e4, \"meter\"))" ] }, { @@ -331,7 +334,7 @@ "metadata": {}, "outputs": [], "source": [ - "# if we run the same request again we will use the cached files (and download none) \n", + "# if we run the same request again we will use the cached files (and download none)\n", "da0 = data_catalog.get_rasterdataset(name, bbox=[11.6, 45.3, 12.0, 46.0])" ] } @@ -352,7 +355,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.5" + "version": "3.11.4" }, "orig_nbformat": 4, "vscode": { diff --git a/hydromt/raster.py b/hydromt/raster.py index 83aa646b2..b7ceb0e7b 100644 --- a/hydromt/raster.py +++ b/hydromt/raster.py @@ -18,7 +18,6 @@ import dask import geopandas as gpd -import mercantile as mct import numpy as np import pandas as pd import pyproj @@ -29,7 +28,6 @@ import xarray as xr import yaml from affine import Affine -from PIL import Image from pyproj import CRS from rasterio import features from rasterio.enums import MergeAlg, Resampling @@ -2250,10 +2248,12 @@ def to_xyz_tiles( def to_slippy_tiles( self, root: Path | str, - reproj_method: str = "nearest", + reproj_method: str = "average", min_lvl: int = None, max_lvl: int = None, driver="png", + cmap=None, + norm=None, **kwargs, ): """Produce tiles in /zoom/x/y. structure (EPSG:3857). @@ -2279,114 +2279,143 @@ def to_slippy_tiles( for GTiff, these are passed to ~:py:meth:hydromt.RasterDataArray.to_raster: for png, these are passed to ~:py:meth:PIL.Image.Image.save: """ - # TODO add option to use cmap instead of elevation2rgba + # for now these are optional dependencies + try: + import matplotlib.pyplot as plt + import mercantile as mct + from PIL import Image + + except ImportError: + raise ImportError("matplotlib, pillow and mercantile are required") + # Fixed pixel size and CRS for XYZ tiles - px_size = 256 + pxs = 256 crs = CRS.from_epsg(3857) ext = {"png": "png", "netcdf4": "nc", "gtiff": "tif"}.get(driver.lower(), None) if ext is None: raise ValueError(f"Unkown file driver {driver}, use png, netcdf4 or GTiff") - # for now we assume output has float32 dtype - kwargs0 = { - "netcdf4": {"engine": "netcdf4"}, - "gtiff": {"driver": "GTiff", "compress": "deflate", "dtype": "float32"}, - }.get(driver.lower(), {}) - kwargs = {**kwargs0, **kwargs} # Object to local variable, also transpose it and extract some meta - if not self._obj.ndim == 2: - ValueError("Only 2d datasets are accepted...") + if not self._obj.ndim == 2 or not isinstance(self._obj, xr.DataArray): + ValueError("Only 2d DataArrays are accepted.") # make sure the data is N-S oriented, with y-axis as first dimension # and nodata values set to nan obj = self.mask_nodata().transpose(self.y_dim, self.x_dim) if obj.raster.res[1] < 0: obj = obj.raster.flipud() + # make sure dataarray has a name + name = obj.name or "data" + obj.name = name + + # colormap output + if cmap is not None and driver != "png": + raise ValueError("Colormap is only supported for png output") + if isinstance(cmap, str): + cmap = plt.get_cmap(cmap) + if cmap is not None and norm is None: + norm = plt.Normalize( + vmin=obj.min().load().item(), vmax=obj.max().load().item() + ) - # Calculate min/max zoomlevel based on the resolution of the data - bounds_wgsg84 = self.transform_bounds("EPSG:4326") - if min_lvl is None: + # for now we assume output has float32 dtype + kwargs0 = { + "netcdf4": {"encoding": {name: {"dtype": "float32", "zlib": True}}}, + "gtiff": {"driver": "GTiff", "compress": "deflate", "dtype": "float32"}, + }.get(driver.lower(), {}) + kwargs = {**kwargs0, **kwargs} + + # 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 min_lat = min(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(min_lat, *obj.raster.res)[0] C = 2 * np.pi * 6378137 # circumference of the earth - min_lvl = int( - np.ceil(np.log2(C * np.cos(np.deg2rad(min_lat)) / dx_m / px_size)) + max_lvl = int( + np.ceil(np.log2(C * np.cos(np.deg2rad(min_lat)) / dx_m / pxs)) ) - if max_lvl is None: - max_lvl = mct.bounding_tile(*bounds_wgsg84).z + if min_lvl is None: # calculate min zoomlevel based on the data extent + min_lvl = mct.bounding_tile(*bounds_wgsg84).z # Loop through the zoom levels zoom_levels = {} logger.info(f"Producing tiles from zoomlevel {min_lvl} to {max_lvl}") - for zl in range(min_lvl, max_lvl - 1, -1): + 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)): + ssd = Path(root, str(zl), f"{tile.x}") + os.makedirs(ssd, exist_ok=True) + tile_bounds = mct.xy_bounds(tile) if i == 0: - ssd = Path(root, str(zl), f"{tile.x}") - os.makedirs(ssd, exist_ok=True) - bounds0 = mct.xy_bounds(tile) - zoom_levels[zl] = abs(bounds0[2] - bounds0[0]) / px_size - if zl == min_lvl: + zoom_levels[zl] = abs(tile_bounds[2] - tile_bounds[0]) / pxs + if zl == max_lvl: # For the first zoomlevel, we can just clip the data # does this need a try/except? src_tile = obj.raster.clip_bbox( - mct.xy_bounds(tile), crs=crs, buffer=2, align=True + tile_bounds, crs=crs, buffer=2, align=True ) else: # Every tile from this level has 4 child tiles on the previous lvl - children = [] - for child in mct.children(tile): + # Create a temporary array, 4 times the size of a tile + temp = np.full((pxs * 2, pxs * 2), np.nan, dtype=np.float64) + for ic, child in enumerate(mct.children(tile)): fn = Path(root, str(child.z), str(child.x), f"{child.y}.{ext}") # Check if the file is really there, if not: it was not written - if fn.exists(): - if driver == "netcdf4": - data = xr.open_dataarray(fn).values - elif driver == "GTiff": - data = xr.open_dataarray( - fn, decode_coords=False, engine="rasterio" - ) - data = data.squeeze(drop=True).values - elif driver == "png": + if not fn.exists(): + continue + # order: top-left, top-right, bottom-right, bottom-left + yslice = slice(0, pxs) if ic in [0, 1] else slice(pxs, None) + xslice = slice(0, pxs) if ic in [0, 3] else slice(pxs, None) + if driver == "netcdf4": + with xr.open_dataset(fn) as ds: + temp[yslice, xslice] = ds[name].values + elif driver == "GTiff": + with rioxarray.open_rasterio( + fn, parse_coordinates=False + ) as da: + temp[yslice, xslice] = da.squeeze(drop=True).values + elif driver == "png": + if cmap is not None: + fn_bin = str(fn).replace(f".{ext}", ".bin") + with open(fn_bin, "r") as f: + data = np.fromfile(f, "f4").reshape((pxs, pxs)) + os.remove(fn_bin) # clean up + else: data = rgba2elevation(np.array(Image.open(fn))) - else: - data = np.full((px_size, px_size), np.nan, dtype=np.float64) - # Add the data to the temporary array - children.append(data) - # Create a temporary array, 4 times the size of a tile - # order of children: top-left, top-right, bottom-right, bottom-left - temp = np.full((px_size * 2, px_size * 2), np.nan, dtype=np.float64) - temp[:px_size, :px_size] = children[0] - temp[:px_size:, px_size:] = children[1] - temp[px_size:, px_size:] = children[2] - temp[px_size:, :px_size] = children[3] + temp[yslice, xslice] = data # create a dataarray from the temporary array src_transform = rasterio.transform.from_bounds( - *mct.xy_bounds(tile), px_size, px_size + *tile_bounds, pxs * 2, pxs * 2 ) src_tile = RasterDataArray.from_numpy( temp, src_transform, crs=crs, nodata=np.nan ) + src_tile.name = name # reproject the data to the tile - dst_transform = rasterio.transform.from_bounds( - *mct.xy_bounds(tile), px_size, px_size - ) + dst_transform = rasterio.transform.from_bounds(*tile_bounds, pxs, pxs) dst_tile = src_tile.raster.reproject( dst_crs=crs, dst_transform=dst_transform, - dst_width=px_size, - dst_height=px_size, + dst_width=pxs, + dst_height=pxs, method=reproj_method, ) # write the data to file fn_out = Path(ssd, f"{tile.y}.{ext}") fns.append(fn_out) if driver.lower() == "png": - # Create RGBA bands from the data and save it as png - rgba = elevation2rgba(dst_tile.values) + if cmap is not None and zl != min_lvl: + # create temp bin file with data for upsampling + fn_bin = str(fn_out).replace(f".{ext}", ".bin") + dst_tile.values.astype("f4").tofile(fn_bin) + rgba = cmap(norm(dst_tile.values), bytes=True) + else: + # Create RGBA bands from the data and save it as png + rgba = elevation2rgba(dst_tile.values) Image.fromarray(rgba).save(fn_out, **kwargs) elif driver.lower() == "netcdf4": # write the data to netcdf diff --git a/pyproject.toml b/pyproject.toml index 5f87c8e24..eb0cc07bb 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -17,7 +17,6 @@ dependencies = [ "entrypoints", # Provide access for Plugins "geopandas>=0.10", # pandas but geo, wraps fiona and shapely "gdal>=3.1", # enable geospacial data manipulation, both raster and victor - "mercantile", # for easilty tiling in XYZ structure (OSM) "numba", # speed up computations (used in e.g. stats) "numpy>=1.20", # pin necessary to ensure compatability with C headers "netcdf4", # netcfd IO @@ -51,11 +50,13 @@ dynamic = ['version', 'description'] [project.optional-dependencies] io = [ - "s3fs", # S3 file system + "fsspec==2023.5.0", # enable cloud file systems likg GCS and S3 "gcsfs", # google cloud file system + "mercantile", # tile handling "openpyxl", # excel IO + "pillow", # image IO "requests", # donwload stuff - "fsspec==2023.5.0", # enable cloud file systems likg GCS and S3 + "s3fs", # S3 file system ] extra = [ "xugrid", # xarray wrapper for mesh processing From f603d8e120127d28ba1af9613617eaeace54f171 Mon Sep 17 00:00:00 2001 From: Dirk Eilander Date: Fri, 22 Sep 2023 09:16:44 +0200 Subject: [PATCH 25/35] test and small bugfixes; update changelog --- docs/changelog.rst | 15 ++++--- hydromt/raster.py | 12 +++--- tests/test_raster.py | 98 ++++++++++++++++++++++++++++++++------------ 3 files changed, 85 insertions(+), 40 deletions(-) diff --git a/docs/changelog.rst b/docs/changelog.rst index 2cc6a9914..a766ba00d 100644 --- a/docs/changelog.rst +++ b/docs/changelog.rst @@ -11,10 +11,8 @@ Unreleased Added ----- -- New raster method ``to_osm_tiles``: tiling of a raster dataset according to the osm structure. -- New raster method ``to_webviewer_tiles``: tiling of a raster dataset for webviewers (png files) -- docs now include a dropdown for selecting older versions of the docs. (#457) -- Support for loading the same data source but from different places (e.g. local & aws) +- docs now include a dropdown for selecting older versions of the docs. (PR #457) +- Support for loading the same data source but from different providers (e.g., local & aws) (PR #438) - Add support for reading and writing tabular data in ``parquet`` format. (PR #445) - Add support for reading model configs in ``TOML`` format. (PR #444) - new ``force-overwrite`` option in ``hydromt update`` CLI to force overwritting updated netcdf files. (PR #460) @@ -22,21 +20,22 @@ Added - Adapters can now clip data that is passed through a python object the same way as through the data catalog. (PR #481) - Model objects now have a _MODEL_VERSION attribute that plugins can use for compatibility purposes (PR # 495) - Model class now has methods for getting, setting, reading and writing arbitrary tabular data. (PR #502) -- Relevant data adapters now have functionality for reporting and detecting the spatial and temporal extent they cover (# 503) +- Relevant data adapters now have functionality for reporting and detecting the spatial and temporal extent they cover (PR #503) - Data catalogs have a ``hydromt_version`` meta key that is used to determine compatibility between the catalog and the installed hydromt version. (PR #506) - Allow the root of a data catalog to point to an archive, this will be extracted to the ~/.hydromt_data folder. (PR #512) - Support for reading overviews from (Cloud Optimized) GeoTIFFs using the zoom_level argument of ``DataCatalog.get_rasterdataset``. (PR #514) - Support for writing overviews to (Cloud Optimized) GeoTIFFs in the ``raster.to_raster`` method. (PR #514) -- Added documentation for how to start your own plugin (#446) +- Added documentation for how to start your own plugin (PR #446) +- New raster method ``to_slippy_tiles``: tiling of a raster dataset according to the slippy tile structure for e.g., webviewers (PR #440). +- Support for http and other *filesystems* in path of data source (PR #515). Changed ------- -- Updated ``MeshModel`` and related methods to support multigrids instead of one single 2D grid. PR #412 +- Updated ``MeshModel`` and related methods to support multigrids instead of one single 2D grid. (PR #412) - possibility to ``load`` the data in the model read_ functions for netcdf files (default for read_grid in r+ mode). (PR #460) - Internal model components (e.g. `Models._maps`, `GridModel._grid``) are now initialized with None and should not be accessed directly, call the corresponding model property (e.g. `Model.maps`, `GridModel.grid`) instead. (PR #473) - Use the Model.data_catalog to read the model region if defined by a geom or grid. (PR #479) -- Support for http and other *filesystems* in path of data source (PR #515). Fixed ----- diff --git a/hydromt/raster.py b/hydromt/raster.py index 226ac03ae..e085e5a79 100644 --- a/hydromt/raster.py +++ b/hydromt/raster.py @@ -2311,8 +2311,8 @@ def to_slippy_tiles( raise ValueError(f"Unkown file driver {driver}, use png, netcdf4 or GTiff") # Object to local variable, also transpose it and extract some meta - if not self._obj.ndim == 2 or not isinstance(self._obj, xr.DataArray): - ValueError("Only 2d DataArrays are accepted.") + if self._obj.ndim != 2: + raise ValueError("Only 2d DataArrays are accepted.") # make sure the data is N-S oriented, with y-axis as first dimension # and nodata values set to nan obj = self.mask_nodata().transpose(self.y_dim, self.x_dim) @@ -2342,14 +2342,14 @@ def to_slippy_tiles( # 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 - min_lat = min(np.abs(bounds_wgsg84[1]), np.abs(bounds_wgsg84[3])) + 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(min_lat, *obj.raster.res)[0] + dx_m = gis_utils.cellres(max_lat, *obj.raster.res)[0] C = 2 * np.pi * 6378137 # circumference of the earth max_lvl = int( - np.ceil(np.log2(C * np.cos(np.deg2rad(min_lat)) / dx_m / pxs)) + np.ceil(np.log2(C * np.cos(np.deg2rad(max_lat)) / dx_m / pxs)) ) if min_lvl is None: # calculate min zoomlevel based on the data extent min_lvl = mct.bounding_tile(*bounds_wgsg84).z @@ -2441,7 +2441,7 @@ def to_slippy_tiles( dst_tile.raster.to_raster(fn_out, **kwargs) if driver.lower() != "png": # Write files to txt and create a vrt using GDAL - txt_fn = Path(root, f"filelist_lvl{zl}.txt") + txt_fn = Path(root, str(zl), "filelist.txt") vrt_fn = Path(root, f"lvl{zl}.vrt") with open(txt_fn, "w") as f: for fn in fns: diff --git a/tests/test_raster.py b/tests/test_raster.py index d5c22a164..e6037f6b7 100644 --- a/tests/test_raster.py +++ b/tests/test_raster.py @@ -2,6 +2,7 @@ """Tests for the hydromt.raster submodule.""" import os +from os.path import isfile, join import dask import geopandas as gpd @@ -517,45 +518,90 @@ def test_to_xyz_tiles(tmpdir, rioda_large): # NOTE: this method does not work in debug mode because of os.subprocess path = str(tmpdir) rioda_large.raster.to_xyz_tiles( - os.path.join(path, "dummy_xyz"), + join(path, "dummy_xyz"), tile_size=256, zoom_levels=[0, 2], ) - with open(os.path.join(path, "dummy_xyz", "0", "filelist.txt"), "r") as f: + with open(join(path, "dummy_xyz", "0", "filelist.txt"), "r") as f: assert len(f.readlines()) == 16 - with open(os.path.join(path, "dummy_xyz", "2", "filelist.txt"), "r") as f: + with open(join(path, "dummy_xyz", "2", "filelist.txt"), "r") as f: assert len(f.readlines()) == 1 test_bounds = [2.13, -2.13, 3.2, -1.07] - _test_r = open_raster(os.path.join(path, "dummy_xyz", "0", "2", "1.tif")) + _test_r = open_raster(join(path, "dummy_xyz", "0", "2", "1.tif")) assert [round(_n, 2) for _n in _test_r.raster.bounds] == test_bounds -# def test_to_osm_tiles(tmpdir, rioda_large): -# path = str(tmpdir) -# rioda_large.raster.to_osm_tiles( -# os.path.join(path, "dummy_osm"), -# zoom_levels=[5], -# min_lvl=8, -# ) -# with open(os.path.join(path, "dummy_osm", "5", "filelist.txt"), "r") as f: -# assert len(f.readlines()) == 1 -# with open(os.path.join(path, "dummy_osm", "8", "filelist.txt"), "r") as f: -# assert len(f.readlines()) == 12 +def test_to_slippy_tiles(tmpdir, rioda_large): + from PIL import Image -# test_bounds = [0.0, -313086.07, 313086.07, -0.0] -# _test_r = open_raster(os.path.join(path, "dummy_osm", "7", "64", "64.tif")) -# assert [round(_n, 2) for _n in _test_r.raster.bounds] == test_bounds + # for tile at zl 7, x 64, y 64 + test_bounds = [0.0, -313086.07, 313086.07, -0.0] + # populate with random data + np.random.seed(0) + rioda_large[:] = np.random.random(rioda_large.shape) - -def test_to_png_tiles(tmpdir, rioda_large): - rioda_large.raster.to_slippy_tiles( - os.path.join(tmpdir, "dummy_png"), - ) - _zl = os.listdir( - os.path.join(tmpdir, "dummy_png"), - ) + # png + png_dir = join(tmpdir, "tiles_png") + rioda_large.raster.to_slippy_tiles(png_dir) + _zl = os.listdir(png_dir) _zl = [int(_n) for _n in _zl] assert len(_zl) == 4 assert min(_zl) == 6 assert max(_zl) == 9 + fn = join(png_dir, "7", "64", "64.png") + im = np.array(Image.open(fn)) + assert im.shape == (256, 256, 4) + assert all(im[0, 0, :] == [128, 0, 131, 255]) + + # test with cmap + png_dir = join(tmpdir, "tiles_png_cmap") + rioda_large.raster.to_slippy_tiles(png_dir, cmap="viridis", min_lvl=7, max_lvl=7) + fn = join(png_dir, "7", "64", "64.png") + im = np.array(Image.open(fn)) + assert im.shape == (256, 256, 4) + assert all(im[0, 0, :] == [128, 0, 132, 255]) + + # gtiff + tif_dir = join(tmpdir, "tiles_tif") + rioda_large.raster.to_slippy_tiles( + tif_dir, + driver="GTiff", + min_lvl=5, + max_lvl=8, + ) + with open(join(tif_dir, "5", "filelist.txt"), "r") as f: + assert len(f.readlines()) == 1 + with open(join(tif_dir, "8", "filelist.txt"), "r") as f: + assert len(f.readlines()) == 12 + _test_r = open_raster(join(tif_dir, "7", "64", "64.tif")) + assert [round(_n, 2) for _n in _test_r.raster.bounds] == test_bounds + assert all([isfile(join(tif_dir, f"lvl{zl}.vrt")) for zl in range(5, 9)]) + _test_vrt = open_raster(join(tif_dir, "lvl7.vrt")) + assert isinstance(_test_vrt, xr.DataArray) + + # nc + nc_dir = join(tmpdir, "tiles_nc") + rioda_large.raster.to_slippy_tiles( + nc_dir, + driver="netcdf4", + min_lvl=5, + max_lvl=8, + ) + with open(join(nc_dir, "5", "filelist.txt"), "r") as f: + assert len(f.readlines()) == 1 + with open(join(nc_dir, "8", "filelist.txt"), "r") as f: + assert len(f.readlines()) == 12 + _test_r = xr.open_dataset(join(nc_dir, "7", "64", "64.nc")) + assert [round(_n, 2) for _n in _test_r.raster.bounds] == test_bounds + assert all([isfile(join(nc_dir, f"lvl{zl}.vrt")) for zl in range(5, 9)]) + _test_vrt = open_raster(join(nc_dir, "lvl7.vrt")) + assert isinstance(_test_vrt, xr.DataArray) + + # test all errors in to_slippy_tiles + with pytest.raises(ValueError, match="Unkown file driver"): + rioda_large.raster.to_slippy_tiles(str(tmpdir), driver="unsupported") + with pytest.raises(ValueError, match="Only 2d DataArrays"): + rioda_large.expand_dims("t").raster.to_slippy_tiles(str(tmpdir)) + with pytest.raises(ValueError, match="Colormap is only supported for png"): + rioda_large.raster.to_slippy_tiles(str(tmpdir), cmap="viridis", driver="GTiff") From b9abc009aad756a3b449bddbf5f3a0e19777ca2f Mon Sep 17 00:00:00 2001 From: Dirk Eilander Date: Fri, 22 Sep 2023 09:34:08 +0200 Subject: [PATCH 26/35] move mercantile back to deps for docs; improve test coverage --- hydromt/raster.py | 13 +++++++------ pyproject.toml | 2 +- tests/test_raster.py | 2 +- 3 files changed, 9 insertions(+), 8 deletions(-) diff --git a/hydromt/raster.py b/hydromt/raster.py index e085e5a79..933f79e69 100644 --- a/hydromt/raster.py +++ b/hydromt/raster.py @@ -18,6 +18,7 @@ import dask import geopandas as gpd +import mercantile as mct import numpy as np import pandas as pd import pyproj @@ -2295,13 +2296,13 @@ def to_slippy_tiles( for png, these are passed to ~:py:meth:PIL.Image.Image.save: """ # for now these are optional dependencies - try: - import matplotlib.pyplot as plt - import mercantile as mct - from PIL import Image + if driver.lower() == "png": + try: + import matplotlib.pyplot as plt + from PIL import Image - except ImportError: - raise ImportError("matplotlib, pillow and mercantile are required") + except ImportError: + raise ImportError("matplotlib and pillow are required for png output") # Fixed pixel size and CRS for XYZ tiles pxs = 256 diff --git a/pyproject.toml b/pyproject.toml index af92132e6..af6c05749 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -18,6 +18,7 @@ dependencies = [ "fsspec", # general file systems utilities "geopandas>=0.10", # pandas but geo, wraps fiona and shapely "gdal>=3.1", # enable geospacial data manipulation, both raster and victor + "mercantile", # tile handling "numba", # speed up computations (used in e.g. stats) "numpy>=1.20", # pin necessary to ensure compatability with C headers "netcdf4", # netcfd IO @@ -52,7 +53,6 @@ dynamic = ['version', 'description'] [project.optional-dependencies] io = [ "gcsfs", # google cloud file system - "mercantile", # tile handling "openpyxl", # excel IO "pillow", # image IO "requests", # donwload stuff diff --git a/tests/test_raster.py b/tests/test_raster.py index e6037f6b7..6d5debdbd 100644 --- a/tests/test_raster.py +++ b/tests/test_raster.py @@ -556,7 +556,7 @@ def test_to_slippy_tiles(tmpdir, rioda_large): # test with cmap png_dir = join(tmpdir, "tiles_png_cmap") - rioda_large.raster.to_slippy_tiles(png_dir, cmap="viridis", min_lvl=7, max_lvl=7) + rioda_large.raster.to_slippy_tiles(png_dir, cmap="viridis", min_lvl=6, max_lvl=7) fn = join(png_dir, "7", "64", "64.png") im = np.array(Image.open(fn)) assert im.shape == (256, 256, 4) From eef06e576eb04b746de9ce458ea083bfde149e0a Mon Sep 17 00:00:00 2001 From: Dirk Eilander Date: Fri, 22 Sep 2023 09:41:35 +0200 Subject: [PATCH 27/35] fix test after including extra zoom lvl --- tests/test_raster.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_raster.py b/tests/test_raster.py index 6d5debdbd..a4ed0b70b 100644 --- a/tests/test_raster.py +++ b/tests/test_raster.py @@ -560,7 +560,7 @@ def test_to_slippy_tiles(tmpdir, rioda_large): fn = join(png_dir, "7", "64", "64.png") im = np.array(Image.open(fn)) assert im.shape == (256, 256, 4) - assert all(im[0, 0, :] == [128, 0, 132, 255]) + assert all(im[0, 0, :] == [31, 148, 139, 255]) # gtiff tif_dir = join(tmpdir, "tiles_tif") From 0c65f5df1fee7372c61aab131679537cd40fe545 Mon Sep 17 00:00:00 2001 From: Dirk Eilander Date: Fri, 22 Sep 2023 09:48:28 +0200 Subject: [PATCH 28/35] remove unecessary flipud --- hydromt/raster.py | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/hydromt/raster.py b/hydromt/raster.py index 5ad189e26..1302ea1a3 100644 --- a/hydromt/raster.py +++ b/hydromt/raster.py @@ -2317,11 +2317,8 @@ def to_slippy_tiles( # Object to local variable, also transpose it and extract some meta if self._obj.ndim != 2: raise ValueError("Only 2d DataArrays are accepted.") - # make sure the data is N-S oriented, with y-axis as first dimension - # and nodata values set to nan + # make sure the y-axis as first dimension and nodata values set to nan obj = self.mask_nodata().transpose(self.y_dim, self.x_dim) - if obj.raster.res[1] < 0: - obj = obj.raster.flipud() # make sure dataarray has a name name = obj.name or "data" obj.name = name @@ -2368,7 +2365,7 @@ def to_slippy_tiles( ssd = Path(root, str(zl), f"{tile.x}") os.makedirs(ssd, exist_ok=True) tile_bounds = mct.xy_bounds(tile) - if i == 0: + if i == 0: # zoom level : resolution in meters zoom_levels[zl] = abs(tile_bounds[2] - tile_bounds[0]) / pxs if zl == max_lvl: # For the first zoomlevel, we can just clip the data @@ -2414,7 +2411,7 @@ def to_slippy_tiles( ) src_tile.name = name - # reproject the data to the tile + # reproject the data to the tile / coares resolution dst_transform = rasterio.transform.from_bounds(*tile_bounds, pxs, pxs) dst_tile = src_tile.raster.reproject( dst_crs=crs, @@ -2443,8 +2440,9 @@ def to_slippy_tiles( elif driver.lower() == "gtiff": # write the data to geotiff dst_tile.raster.to_raster(fn_out, **kwargs) + + # Write files to txt and create a vrt using GDAL if driver.lower() != "png": - # Write files to txt and create a vrt using GDAL txt_fn = Path(root, str(zl), "filelist.txt") vrt_fn = Path(root, f"lvl{zl}.vrt") with open(txt_fn, "w") as f: @@ -2452,9 +2450,8 @@ def to_slippy_tiles( f.write(f"{fn}\n") gis_utils.create_vrt(vrt_fn, file_list_path=txt_fn) + # Write a quick yaml for the database if driver.lower() != "png": - # Write a quick yaml for the database - # zoom level : resolution in meters yml = { "crs": 3857, "data_type": "RasterDataset", From 1f1a0ed957b3f450c6fadd9c8278228a3691b3a3 Mon Sep 17 00:00:00 2001 From: Sam Vente Date: Fri, 22 Sep 2023 10:24:23 +0200 Subject: [PATCH 29/35] fix stack level --- hydromt/data_catalog.py | 1 + 1 file changed, 1 insertion(+) diff --git a/hydromt/data_catalog.py b/hydromt/data_catalog.py index 521c2932d..31c53b051 100644 --- a/hydromt/data_catalog.py +++ b/hydromt/data_catalog.py @@ -428,6 +428,7 @@ def __contains__(self, key: str) -> bool: "Directly checking for containement is deprecated. " " Use 'contains_source' instead.", DeprecationWarning, + stacklevel=2, ) return self.contains_source(key) From d1f432ebec47c6212a8629d9a1fa0fbe4ac91dd5 Mon Sep 17 00:00:00 2001 From: Brendan Date: Fri, 22 Sep 2023 14:31:27 +0200 Subject: [PATCH 30/35] Updated version --- examples/working_with_tiled_raster_data.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/working_with_tiled_raster_data.ipynb b/examples/working_with_tiled_raster_data.ipynb index f7cff62e1..777a57a6d 100644 --- a/examples/working_with_tiled_raster_data.ipynb +++ b/examples/working_with_tiled_raster_data.ipynb @@ -355,7 +355,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.4" + "version": "3.11.5" }, "orig_nbformat": 4, "vscode": { From 56c4b4ef2e911da330039337b9f79c87c33a32e7 Mon Sep 17 00:00:00 2001 From: Brendan Date: Fri, 22 Sep 2023 14:31:45 +0200 Subject: [PATCH 31/35] Revamped zoom level determination --- hydromt/raster.py | 63 ++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 51 insertions(+), 12 deletions(-) diff --git a/hydromt/raster.py b/hydromt/raster.py index e085e5a79..31b58e7cc 100644 --- a/hydromt/raster.py +++ b/hydromt/raster.py @@ -2305,6 +2305,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: @@ -2316,11 +2320,13 @@ def to_slippy_tiles( # make sure the data is N-S oriented, with y-axis as first dimension # and nodata values set to nan obj = self.mask_nodata().transpose(self.y_dim, self.x_dim) - if obj.raster.res[1] < 0: - obj = obj.raster.flipud() + # if obj.raster.res[1] < 0: + # obj = obj.raster.flipud() # 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": @@ -2339,20 +2345,53 @@ 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 = {} @@ -2360,7 +2399,7 @@ def to_slippy_tiles( 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) From e6e1162bc84eae6a19294b6f1b51b0d13ad15cfd Mon Sep 17 00:00:00 2001 From: Brendan Date: Fri, 22 Sep 2023 16:00:28 +0200 Subject: [PATCH 32/35] Added extra example for the oil painting --- examples/working_with_tiled_raster_data.ipynb | 58 +++++++++++++++++-- 1 file changed, 54 insertions(+), 4 deletions(-) diff --git a/examples/working_with_tiled_raster_data.ipynb b/examples/working_with_tiled_raster_data.ipynb index 777a57a6d..052d73eb3 100644 --- a/examples/working_with_tiled_raster_data.ipynb +++ b/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)" ] }, { From a0ba1245f7105d785dc0e95929369c2b11724cb1 Mon Sep 17 00:00:00 2001 From: Brendan Date: Fri, 22 Sep 2023 16:00:44 +0200 Subject: [PATCH 33/35] some docs for the colormap --- hydromt/raster.py | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/hydromt/raster.py b/hydromt/raster.py index 31b58e7cc..739323eae 100644 --- a/hydromt/raster.py +++ b/hydromt/raster.py @@ -2267,8 +2267,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. structure (EPSG:3857). @@ -2288,6 +2288,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: @@ -2345,7 +2352,7 @@ def to_slippy_tiles( }.get(driver.lower(), {}) kwargs = {**kwargs0, **kwargs} - # Setting up information for zoomlevel calculation and + # 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")) From 81355481b250594cb732469a262ddcc4a683d58a Mon Sep 17 00:00:00 2001 From: Brendan Date: Fri, 22 Sep 2023 16:45:18 +0200 Subject: [PATCH 34/35] Restored yaml filename to folder; pinned zl in example --- examples/working_with_tiled_raster_data.ipynb | 2 +- hydromt/raster.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/working_with_tiled_raster_data.ipynb b/examples/working_with_tiled_raster_data.ipynb index 052d73eb3..8cc01afd9 100644 --- a/examples/working_with_tiled_raster_data.ipynb +++ b/examples/working_with_tiled_raster_data.ipynb @@ -308,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" ] }, diff --git a/hydromt/raster.py b/hydromt/raster.py index 73e7cc8f9..4c10d96c5 100644 --- a/hydromt/raster.py +++ b/hydromt/raster.py @@ -2506,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( From 2155a2f9eb55f48440857366071f1360326c988f Mon Sep 17 00:00:00 2001 From: Brendan Date: Fri, 22 Sep 2023 18:34:38 +0200 Subject: [PATCH 35/35] Some last simplifications --- hydromt/raster.py | 43 +++++++++++++++++-------------------------- 1 file changed, 17 insertions(+), 26 deletions(-) diff --git a/hydromt/raster.py b/hydromt/raster.py index 4c10d96c5..11bff366a 100644 --- a/hydromt/raster.py +++ b/hydromt/raster.py @@ -2294,7 +2294,7 @@ def to_slippy_tiles( 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. + via that string or as a Colormap object from matplotlib itself. norm : object, optional A matplotlib Normalize object that defines a range between a maximum and minimum value @@ -2353,31 +2353,7 @@ 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 - + # Setting up information for determination of tile windows # This section is for dealing with rounding errors obj_bounds_cor = [ obj_bounds[0] + 0.5 * obj_res, @@ -2394,6 +2370,21 @@ def to_slippy_tiles( # Calculate min/max zoomlevel based if max_lvl is None: # calculate max zoomlevel close to native resolution # Determine the max number of zoom levels with the resolution + # This section is purely for the resolution + obj_clipped_to_pseudo = obj.raster.clip_bbox( + (-180, -y_ext, 180, y_ext), + 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 + + # Calculate the maximum zoom level dres = tr_3857[0] max_lvl = int( math.ceil((math.log10((y_ext_pm * 2) / (dres * pxs)) / math.log10(2)))