diff --git a/cht_tiling/tiled_web_map.py b/cht_tiling/tiled_web_map.py index 9bb09b4..c996124 100644 --- a/cht_tiling/tiled_web_map.py +++ b/cht_tiling/tiled_web_map.py @@ -15,7 +15,9 @@ from botocore import UNSIGNED from botocore.client import Config -from cht_tiling.topobathy import make_topobathy_tiles +from cht_tiling.topobathy import make_topobathy_tiles_top_level, make_topobathy_tiles_lower_levels +from cht_tiling.webviewer import write_html + from cht_tiling.utils import ( get_zoom_level_for_resolution, list_files, @@ -43,6 +45,7 @@ def __init__(self, path, name, parameter="elevation"): self.name = name self.path = path self.url = None + self.npix = 256 self.parameter = "elevation" self.encoder = "terrarium" self.encoder_vmin = None @@ -109,18 +112,18 @@ def get_data(self, xl, yl, max_pixel_size, crs=None, waitbox=None): izoom = min(izoom, self.max_zoom) # Determine the indices of required tiles - ix0, it0 = xy2num(xl[0], yl[1], izoom) - ix1, it1 = xy2num(xl[1], yl[0], izoom) + ix0, iy0 = xy2num(xl[0], yl[1], izoom) + ix1, iy1 = xy2num(xl[1], yl[0], izoom) # Make sure indices are within bounds ix0 = max(0, ix0) - it0 = max(0, it0) + it0 = max(0, iy0) ix1 = min(2**izoom - 1, ix1) - it1 = min(2**izoom - 1, it1) + iy1 = min(2**izoom - 1, iy1) # Create empty array nx = (ix1 - ix0 + 1) * 256 - ny = (it1 - it0 + 1) * 256 + ny = (iy1 - iy0 + 1) * 256 z = np.empty((ny, nx)) z[:] = np.nan @@ -131,7 +134,7 @@ def get_data(self, xl, yl, max_pixel_size, crs=None, waitbox=None): if self.download: for i in range(ix0, ix1 + 1): ifolder = str(i) - for j in range(it0, it1 + 1): + for j in range(it0, iy1 + 1): png_file = os.path.join( self.path, str(izoom), ifolder, str(j) + ".png" ) @@ -175,7 +178,7 @@ def get_data(self, xl, yl, max_pixel_size, crs=None, waitbox=None): # Loop over required tiles for i in range(ix0, ix1 + 1): ifolder = str(i) - for j in range(it0, it1 + 1): + for j in range(iy0, iy1 + 1): png_file = os.path.join(self.path, str(izoom), ifolder, str(j) + ".png") if not os.path.exists(png_file): @@ -192,13 +195,13 @@ def get_data(self, xl, yl, max_pixel_size, crs=None, waitbox=None): # Fill array ii0 = (i - ix0) * 256 ii1 = ii0 + 256 - jj0 = (j - it0) * 256 + jj0 = (j - iy0) * 256 jj1 = jj0 + 256 z[jj0:jj1, ii0:ii1] = valt # Compute x and y coordinates - x0, y0 = num2xy(ix0, it1 + 1, izoom) # lower left - x1, y1 = num2xy(ix1 + 1, it0, izoom) # upper right + x0, y0 = num2xy(ix0, iy1 + 1, izoom) # lower left + x1, y1 = num2xy(ix1 + 1, iy0, izoom) # upper right # Data is stored in centres of pixels so we need to shift the coordinates dx = (x1 - x0) / nx dy = (y1 - y0) / ny @@ -208,11 +211,46 @@ def get_data(self, xl, yl, max_pixel_size, crs=None, waitbox=None): return x, y, z - def generate_topobathy_tiles(self, **kwargs): - make_topobathy_tiles(self.path, **kwargs) - if "make_availability_file" in kwargs: - if kwargs["make_availability_file"]: - self.make_availability_file() + def generate_topobathy_tiles(self, + datalist, + index_path=None, + lon_range=None, + lat_range=None, + zoom_range=None, + make_webviewer=True, + write_metadata=True, + make_availability_file=True, + make_lower_levels=True, + make_highest_level=True, + skip_existing=False, + parallel=True, + interpolation_method="linear"): + + if make_highest_level: + for data_dict in datalist: + # Can loop here around differet datasets + make_topobathy_tiles_top_level(self, + data_dict, + index_path=index_path, + lon_range=lon_range, + lat_range=lat_range, + zoom_range=zoom_range, + skip_existing=skip_existing, + parallel=parallel, + interpolation_method=interpolation_method) + + if make_lower_levels: + make_topobathy_tiles_lower_levels(self, + skip_existing=skip_existing, + parallel=parallel) + + # For anything but global datasets, make an availability file + if make_availability_file: + self.make_availability_file() + if write_metadata: + self.write_metadata() + if make_webviewer: + write_html(os.path.join(self.path, "index.html"), max_native_zoom=self.max_zoom) def check_availability(self, i, j, izoom): # Check if tile exists at all @@ -326,3 +364,50 @@ def make_availability_file(self): ds.to_netcdf(nc_file) ds.close() + + def write_metadata(self): + + metadata = {} + + metadata["longname"] = self.long_name + metadata["format"] = "tiled_web_map" + metadata["encoder"] = self.encoder + if self.encoder_vmin is not None: + metadata["encoder_vmin"] = self.encoder_vmin + if self.encoder_vmax is not None: + metadata["encoder_vmax"] = self.encoder_vmax + metadata["url"] = self.url + if self.s3_bucket is not None: + metadata["s3_bucket"] = self.s3_bucket + if self.s3_key is not None: + metadata["s3_key"] = self.s3_key + if self.s3_region is not None: + metadata["s3_region"] = self.s3_region + metadata["max_zoom"] = self.max_zoom + metadata["interpolation_method"] = "linear" + metadata["source"] = self.source + metadata["coord_ref_sys_name"] = "WGS 84 / Pseudo-Mercator" + metadata["coord_ref_sys_kind"] = "projected" + metadata["vertical_reference_level"] = self.vertical_reference_level + metadata["vertical_units"] = self.vertical_units + metadata["difference_with_msl"] = self.difference_with_msl + metadata["available_tiles"] = True + + metadata["description"] = {} + metadata["description"]["title"] = "LONG_NAME" + metadata["description"]["institution"] = "INST_NAME" + metadata["description"]["history"] = "created by : AUTHOR" + metadata["description"]["references"] = "No reference material available" + metadata["description"]["comment"] = "none" + metadata["description"]["email"] = "Your email here" + metadata["description"]["version"] = "1.0" + metadata["description"]["terms_for_use"] = "Use as you like" + metadata["description"]["disclaimer"] = ( + "These data are made available in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE." + ) + + # Write to toml file + # toml_file = os.path.join(path, name + ".tml") + toml_file = os.path.join(self.path, "metadata.tml") + with open(toml_file, "w") as f: + toml.dump(metadata, f) diff --git a/cht_tiling/topobathy.py b/cht_tiling/topobathy.py index f80f523..59ef034 100644 --- a/cht_tiling/topobathy.py +++ b/cht_tiling/topobathy.py @@ -1,9 +1,12 @@ import os import time from multiprocessing.pool import ThreadPool - import numpy as np +import xarray as xr +import rioxarray import toml + + from cht_utils.misc_tools import interp2 from pyproj import CRS, Transformer @@ -15,150 +18,21 @@ num2deg, png2elevation, ) -from cht_tiling.webviewer import write_html - - -def make_lower_level_tile( - zoom_path_i, - zoom_path_higher, - i, - j, - npix, - encoder, - encoder_vmin, - encoder_vmax, - compress_level, -): - # Get the indices of the tiles in the higher zoom level - i00, j00 = 2 * i, 2 * j # upper left - i10, j10 = 2 * i, 2 * j + 1 # lower left - i01, j01 = 2 * i + 1, 2 * j # upper right - i11, j11 = 2 * i + 1, 2 * j + 1 # lower right - - # Create empty array of NaN to store the elevation data from the higher zoom level - zg512 = np.zeros((npix * 2, npix * 2)) - zg512[:] = np.nan +# from cht_tiling.webviewer import write_html - # Create empty array of NaN of 4*npix*npix to store the 2-strid elevation data from higher zoom level - zg4 = np.zeros((4, npix, npix)) - zg4[:] = np.nan - - okay = False - - # Get the file names of the tiles in the higher zoom level - # Upper left - file_name = os.path.join(zoom_path_higher, str(i00), str(j00) + ".png") - if os.path.exists(file_name): - zgh = png2elevation( - file_name, - encoder=encoder, - encoder_vmin=encoder_vmin, - encoder_vmax=encoder_vmax, - ) - zg512[0:npix, 0:npix] = zgh - okay = True - # Lower left - file_name = os.path.join(zoom_path_higher, str(i10), str(j10) + ".png") - if os.path.exists(file_name): - zgh = png2elevation( - file_name, - encoder=encoder, - encoder_vmin=encoder_vmin, - encoder_vmax=encoder_vmax, - ) - zg512[npix:, 0:npix] = zgh - okay = True - # Upper right - file_name = os.path.join(zoom_path_higher, str(i01), str(j01) + ".png") - if os.path.exists(file_name): - zgh = png2elevation( - file_name, - encoder=encoder, - encoder_vmin=encoder_vmin, - encoder_vmax=encoder_vmax, - ) - zg512[0:npix, npix:] = zgh - okay = True - # Lower right - file_name = os.path.join(zoom_path_higher, str(i11), str(j11) + ".png") - if os.path.exists(file_name): - zgh = png2elevation( - file_name, - encoder=encoder, - encoder_vmin=encoder_vmin, - encoder_vmax=encoder_vmax, - ) - zg512[npix:, npix:] = zgh - okay = True - - if not okay: - # No tiles in higher zoom level, so continue - return - - # Compute average of 4 tiles in higher zoom level - # Data from zg512 with stride 2 - zg4[0, :, :] = zg512[0 : npix * 2 : 2, 0 : npix * 2 : 2] - zg4[1, :, :] = zg512[1 : npix * 2 : 2, 0 : npix * 2 : 2] - zg4[2, :, :] = zg512[0 : npix * 2 : 2, 1 : npix * 2 : 2] - zg4[3, :, :] = zg512[1 : npix * 2 : 2, 1 : npix * 2 : 2] - - # Compute average of 4 tiles - zg = np.nanmean(zg4, axis=0) - - # Write to terrarium png format - file_name = os.path.join(zoom_path_i, str(j) + ".png") - elevation2png( - zg, - file_name, - encoder=encoder, - encoder_vmin=encoder_vmin, - encoder_vmax=encoder_vmax, - compress_level=compress_level, - ) - - -def make_topobathy_tiles( - path, - dem_names=None, - dem_list=None, - bathymetry_database=None, - dataset=None, # Must be XArray - dataarray_name="elevation", - dataarray_x_name="lon", - dataarray_y_name="lat", +def make_topobathy_tiles_top_level( + twm, + data_dict, + index_path=None, lon_range=None, lat_range=None, - index_path=None, zoom_range=None, - dx_max_zoom=None, - z_range=None, - bathymetry_database_path="d:\\delftdashboard\\data\\bathymetry", - quiet=False, - make_webviewer=True, - write_metadata=True, - make_availability_file=True, - metadata=None, - make_lower_levels=True, - make_highest_level=True, skip_existing=False, + parallel=True, interpolation_method="linear", - encoder="terrarium", - encoder_vmin=None, - encoder_vmax=None, - compress_level=6, - name="unknown", - long_name="unknown", - url=None, - s3_bucket=None, - s3_key=None, - s3_region=None, - source="unknown", - vertical_reference_level="MSL", - vertical_units="m", - difference_with_msl=0.0, ): """ - Generates topo/bathy tiles + Generates highest zoom level topo/bathy tiles :param path: Path where topo/bathy tiles will be stored. :type path: str @@ -174,330 +48,258 @@ def make_topobathy_tiles( """ - if dem_names is not None: - dem_list = [] - for dem_name in dem_names: - dem = {} - dem["name"] = dem_name - dem["zmin"] = -10000.0 - dem["zmax"] = 10000.0 - dem_list.append(dem) - - if dem_list: - dem_type = "ddb" - - if bathymetry_database is None: - from cht_bathymetry.bathymetry_database import BathymetryDatabase - - bathymetry_database = BathymetryDatabase(None) - bathymetry_database.initialize(bathymetry_database_path) - - if lon_range is None and index_path is None: - # Try to get lon_range and lat_range from the first dataset - dataset = bathymetry_database.get_dataset(dem_list[0]["name"]) - lon_range, lat_range = dataset.get_bbox(crs=CRS.from_epsg(4326)) - - elif dataset is not None: - dem_type = "xarray" - ds_lon = dataset[dataarray_x_name].to_numpy() - ds_lat = dataset[dataarray_y_name].to_numpy() - ds_z_parameter = dataarray_name - - if "crs" not in dataset: - dataset_crs_code = 4326 + npix = 256 + + transformer_4326_to_3857 = Transformer.from_crs( + CRS.from_epsg(4326), CRS.from_epsg(3857), always_xy=True + ) + + # Set vertical range + zmin = -20000.0 + zmax = 20000.0 + if "zmin" in data_dict: + zmin = data_dict["zmin"] + if "zmax" in data_dict: + zmax = data_dict["zmax"] + z_range = [zmin, zmax] + + # First we determine lon_range, lat_range, dx and crs_data + # lon_range and lat_range are in geographic coordinates and are used to determine the bounding box of the tiles + # dx is the resolution of the dataset in meters and is used to determine the zoom levels + # crs_data is the CRS of the dataset and is used for the transformation between this crs and web mercator + + if "twm" in data_dict: + input_twm = data_dict["twm"] + data_type = "twm" + + elif "ncfile" in data_dict: + # Read the netcdf file and get the data array + ds = xr.open_dataset(data_dict["ncfile"]) + da = ds[data_dict["dataarray_name"]] + data_type = "xarray" + + elif "da" in data_dict: + # XArray dataarray as input + da = data_dict["da"] + data_type = "xarray" + + elif "tiffile" in data_dict: + da = rioxarray.open_rasterio(data_dict["tiffile"]) + data_type = "xarray" + + if data_type == "xarray": + + # Get the CRS + crs_data = da.rio.crs + + # Get the lon and lat range and dx + if crs_data.is_geographic: + if not lon_range: + lon_range = [da.x.min(), da.x.max()] + if not lat_range: + lat_range = [da.y.min(), da.y.max()] + dx = np.mean(np.diff(da.y)) * 111000.0 + else: - dataset_crs_code = dataset.crs.attrs["epsg_code"] - transformer_3857_to_crs = Transformer.from_crs( - CRS.from_epsg(3857), dataset_crs_code, always_xy=True - ) - if lon_range is None: - lon_range = [np.min(ds_lon), np.max(ds_lon)] - if lat_range is None: - lat_range = [np.min(ds_lat), np.max(ds_lat)] - dx = np.mean(np.diff(ds_lat)) * 111000.0 + if not lon_range and not lat_range: + # Get bounding box + x_range = [da.x.min(), da.x.max()] + y_range = [da.y.min(), da.y.max()] + # Get lon/lat range + transformer = Transformer.from_crs(crs_data, CRS.from_epsg(4326), always_xy=True) + lon_range = [transformer.transform(x_range[0], y_range[0])[0], transformer.transform(x_range[1], y_range[1])[0]] + lat_range = [transformer.transform(x_range[0], y_range[0])[1], transformer.transform(x_range[1], y_range[1])[1]] + dx = np.abs(np.mean(np.diff(da.y))) + + transformer_3857_to_crs = Transformer.from_crs( + CRS.from_epsg(3857), crs_data, always_xy=True + ) + # Determine zoom range if zoom_range is None: - if dx_max_zoom is None: - dx_max_zoom = dx * 0.5 - zoom_max = get_zoom_level_for_resolution(dx_max_zoom) + zoom_max = get_zoom_level_for_resolution(dx * 0.5) zoom_range = [0, zoom_max] + twm.zoom_range = [0, zoom_max] + twm.max_zoom = zoom_range[1] + + # Use highest zoom level + izoom = zoom_range[1] + zoom_path = os.path.join(twm.path, str(izoom)) + + # Determine elapsed time + t0 = time.time() + + # Create rectangular mesh with origin (0.0) and 256x256 pixels + dxy = (40075016.686 / npix) / 2**izoom + xx = np.linspace(0.0, (npix - 1) * dxy, num=npix) + yy = xx[:] + xv, yv = np.meshgrid(xx, yy) + + # Determine min and max indices for this zoom level + # If the index path is given, then get ix0, ix1, iy0 and iy1 from the existing index files + # Otherwise, use lon_range and lat_range + # This option is used when we only want to make tiles that cover a model domain + if index_path: + index_zoom_path = os.path.join(index_path, str(izoom)) + # List folders and turn names into integers + iy0 = 1e15 + iy1 = -1e15 + ix_list = [int(i) for i in os.listdir(index_zoom_path)] + ix0 = min(ix_list) + ix1 = max(ix_list) + # Now loop through the folders to get the min and max y indices + for i in range(ix0, ix1 + 1): + it_list = [ + int(j.split(".")[0]) + for j in os.listdir(os.path.join(index_zoom_path, str(i))) + ] + iy0 = min(iy0, min(it_list)) + iy1 = max(iy1, max(it_list)) + else: + # Use the lon_range and lat_range from the data array or from user input + ix0, iy0 = deg2num(lat_range[1], lon_range[0], izoom) + ix1, iy1 = deg2num(lat_range[0], lon_range[1], izoom) + + # Limit the indices + ix0 = max(0, ix0) + iy0 = max(0, iy0) + ix1 = min(2**izoom - 1, ix1) + iy1 = min(2**izoom - 1, iy1) + + # Add some stuff to options dict, which is used for parallel processing + options = {} + options["index_path"] = index_path + options["transformer_4326_to_3857"] = transformer_4326_to_3857 + options["transformer_3857_to_crs"] = transformer_3857_to_crs + options["xv"] = xv + options["yv"] = yv + options["dxy"] = dxy + options["interpolation_method"] = interpolation_method + options["z_range"] = z_range + + # Loop in x direction + for i in range(ix0, ix1 + 1): + + print(f"Processing column {i - ix0 + 1} of {ix1 - ix0 + 1}") + + zoom_path_i = os.path.join(zoom_path, str(i)) + + if not os.path.exists(zoom_path_i): + makedir(zoom_path_i) + + # Loop in y direction + if parallel: + with ThreadPool() as pool: + pool.starmap( + create_highest_zoom_level_tile, + [ + ( + zoom_path_i, + i, + j, + izoom, + twm, + da, + data_type, + options, + ) + for j in range(iy0, iy1 + 1) + ], + ) + else: + # Loop in y direction + for j in range(iy0, iy1 + 1): + # Create highest zoom level tile + create_highest_zoom_level_tile(zoom_path_i, i, j, izoom, twm, da, data_type, options) + + # If zoom_path_i is empty, then remove it again + if not os.listdir(zoom_path_i): + os.rmdir(zoom_path_i) - if not z_range: - z_range = [-20000.0, 20000.0] + t1 = time.time() + + print("Elapsed time for zoom level " + str(izoom) + ": " + str(t1 - t0)) + + # Done with highest zoom level + +def make_topobathy_tiles_lower_levels( + twm, + skip_existing=False, + parallel=True, +): npix = 256 - if make_highest_level: - transformer_4326_to_3857 = Transformer.from_crs( - CRS.from_epsg(4326), CRS.from_epsg(3857), always_xy=True - ) + # Now loop through other zoom levels starting with highest minus 1 + + for izoom in range(twm.zoom_range[1] - 1, twm.zoom_range[0] - 1, -1): - # First do highest zoom level - izoom = zoom_range[1] + print("Processing zoom level " + str(izoom)) # Determine elapsed time t0 = time.time() - if not quiet: - print("Processing zoom level " + str(izoom)) - - zoom_path = os.path.join(path, str(izoom)) - - dxy = (40075016.686 / npix) / 2**izoom - xx = np.linspace(0.0, (npix - 1) * dxy, num=npix) - yy = xx[:] - xv, yv = np.meshgrid(xx, yy) - - # Determine min and max indices for this zoom level - # If the index path is given, then get ix0, ix1, it0 and it1 from the existing index files - # Otherwise, use lon_range and lat_range - - if index_path: - index_zoom_path = os.path.join(index_path, str(izoom)) - # List folders and turn names into integers - it0 = 1e15 - it1 = -1e15 - ix_list = [int(i) for i in os.listdir(index_zoom_path)] - ix0 = min(ix_list) - ix1 = max(ix_list) - # Now loop through the folders to get the min and max y indices - for i in range(ix0, ix1 + 1): - it_list = [ - int(j.split(".")[0]) - for j in os.listdir(os.path.join(index_zoom_path, str(i))) - ] - it0 = min(it0, min(it_list)) - it1 = max(it1, max(it_list)) - else: - ix0, it0 = deg2num(lat_range[1], lon_range[0], izoom) - ix1, it1 = deg2num(lat_range[0], lon_range[1], izoom) - - ix0 = max(0, ix0) - it0 = max(0, it0) - ix1 = min(2**izoom - 1, ix1) - it1 = min(2**izoom - 1, it1) + # Rather than interpolating the data onto tiles, we will take average of 4 tiles in higher zoom level + + zoom_path = os.path.join(twm.path, str(izoom)) + zoom_path_higher = os.path.join(twm.path, str(izoom + 1)) + + # # + # if index_path: + # index_zoom_path = os.path.join(index_path, str(izoom)) + # # List folders and turn names into integers + # iy0 = 1e15 + # iy1 = -1e15 + # ix_list = [int(i) for i in os.listdir(index_zoom_path)] + # ix0 = min(ix_list) + # ix1 = max(ix_list) + # # Now loop through the folders to get the min and max y indices + # for i in range(ix0, ix1 + 1): + # iy_list = [ + # int(j.split(".")[0]) + # for j in os.listdir(os.path.join(index_zoom_path, str(i))) + # ] + # iy0 = min(iy0, min(iy_list)) + # iy1 = max(iy1, max(iy_list)) + # else: + # ix0, iy0 = deg2num(lat_range[1], lon_range[0], izoom) + # iy1, iy1 = deg2num(lat_range[0], lon_range[1], izoom) + + # ix0 = max(0, ix0) + # iy0 = max(0, iy0) + # ix1 = min(2**izoom - 1, ix1) + # iy1 = min(2**izoom - 1, iy1) + + # First determine ix0 and ix1 based on higher zoom level + # Get the folders in zoom_path_higher and turn them into integers + ix_list = [int(i) for i in os.listdir(zoom_path_higher)] + ix0_higher = min(ix_list) + ix1_higher = max(ix_list) + ix0 = int(ix0_higher / 2) + ix1 = int(ix1_higher / 2) + + # Now loop through the folders to get the min and max y indices + iy0_higher = 1e15 + iy1_higher = -1e15 + for i in os.listdir(zoom_path_higher): + iy_list = [int(j.split(".")[0]) for j in os.listdir(os.path.join(zoom_path_higher, i))] + iy0_higher = min(iy0_higher, min(iy_list)) + iy1_higher = max(iy1_higher, max(iy_list)) + iy0 = int(iy0_higher / 2) + iy1 = int(iy1_higher / 2) # Loop in x direction - filename_list = [] - i_list = [] - j_list = [] - for i in range(ix0, ix1 + 1): - print(f"Processing column {i - ix0 + 1} of {ix1 - ix0 + 1}") path_okay = False zoom_path_i = os.path.join(zoom_path, str(i)) - if not os.path.exists(zoom_path_i): - makedir(zoom_path_i) - path_okay = True - - # Loop in y direction - for j in range(it0, it1 + 1): - # Create highest zoom level tile - - file_name = os.path.join(zoom_path_i, str(j) + ".png") - if os.path.exists(file_name): - if skip_existing: - # Tile already exists - continue - else: - # Read the tile - zg0 = png2elevation( - file_name, - encoder=encoder, - encoder_vmin=encoder_vmin, - encoder_vmax=encoder_vmax, - ) - pass - else: - # Tile does not exist - zg0 = np.zeros((npix, npix)) - zg0[:] = np.nan - - # if index_path: - # # Only make tiles for which there is an index file - # index_file_name = os.path.join( - # index_path, str(izoom), str(i), str(j) + ".png" - # ) - # if not os.path.exists(index_file_name): - # continue - - # # Compute lat/lon at upper left corner of tile - # lat, lon = num2deg(i, j, izoom) - - # # Convert origin to Global Mercator - # xo, yo = transformer_4326_to_3857.transform(lon, lat) - - # # Tile grid on Global mercator - # x3857 = xo + xv[:] + 0.5 * dxy - # y3857 = yo - yv[:] - 0.5 * dxy - - # if dem_type == "ddb": - # zg = bathymetry_database.get_bathymetry_on_grid( - # x3857, y3857, CRS.from_epsg(3857), dem_list - # ) - # elif dem_type == "xarray": - # # Make grid of x3857 and y3857, and convert to crs of dataset - # # xg, yg = np.meshgrid(x3857, y3857) - # xg, yg = transformer_3857_to_crs.transform(x3857, y3857) - # # Get min and max of xg, yg - # xg_min = np.min(xg) - # xg_max = np.max(xg) - # yg_min = np.min(yg) - # yg_max = np.max(yg) - # # Add buffer to grid - # dbuff = 0.05 * max(xg_max - xg_min, yg_max - yg_min) - # xg_min = xg_min - dbuff - # xg_max = xg_max + dbuff - # yg_min = yg_min - dbuff - # yg_max = yg_max + dbuff - - # # Get the indices of the dataset that are within the xg, yg range - # i0 = np.where(ds_lon <= xg_min)[0] - # if len(i0) == 0: - # # Take first index - # i0 = 0 - # else: - # # Take last index - # i0 = i0[-1] - # i1 = np.where(ds_lon >= xg_max)[0] - # if len(i1) == 0: - # i1 = len(ds_lon) - 1 - # else: - # i1 = i1[0] - # if i1 <= i0: - # # No data for this tile - # continue - - # xd = ds_lon[i0:i1] - - # if ds_lat[0] < ds_lat[-1]: - # # South to North - # j0 = np.where(ds_lat <= yg_min)[0] - # if len(j0) == 0: - # j0 = 0 - # else: - # j0 = j0[-1] - # j1 = np.where(ds_lat >= yg_max)[0] - # if len(j1) == 0: - # j1 = len(ds_lat) - 1 - # else: - # j1 = j1[0] - # if j1 <= j0: - # # No data for this tile - # continue - # # Get the dataset within the range - # yd = ds_lat[j0:j1] - # # Get number of dimensions of dataarray - # if len(dataset[ds_z_parameter].shape) == 2: - # zd = dataset[ds_z_parameter][j0:j1, i0:i1].values[:] - # else: - # zd = np.squeeze(dataset[ds_z_parameter][0, j0:j1, i0:i1].values[:]) - # else: - # # North to South - # j0 = np.where(ds_lat <= yg_min)[0] - # if len(j0) == 0: - # # Use last index - # j0 = len(ds_lat) - 1 - # else: - # # Use first index - # j0 = j0[0] - # j1 = np.where(ds_lat >= yg_max)[0] - # if len(j1) == 0: - # j1 = 0 - # else: - # j1 = j1[-1] - # if j0 <= j1: - # # No data for this tile - # continue - # # Get the dataset within the range - # yd = np.flip(ds_lat[j1:j0]) - # if len(dataset[ds_z_parameter].shape) == 2: - # zd = np.flip(dataset[ds_z_parameter][j1:j0, i0:i1].values[:], axis=0) - # else: - # zd = np.squeeze(np.flip(dataset[ds_z_parameter][0, j1:j0, i0:i1].values[:], axis=0)) - # zg = interp2(xd, yd, zd, xg, yg, method=interpolation_method) - - # if np.isnan(zg).all(): - # # only nans in this tile - # continue - - # if np.nanmax(zg) < z_range[0] or np.nanmin(zg) > z_range[1]: - # # all values in tile outside z_range - # continue - - # # Overwrite zg with zg0 where not nan - # mask = np.isnan(zg) - # zg[mask] = zg0[mask] - - # Write to terrarium png format - elevation2png( - zg, - file_name, - compress_level=compress_level, - encoder=encoder, - encoder_vmin=encoder_vmin, - encoder_vmax=encoder_vmax, - ) - - t1 = time.time() - - if not quiet: - print("Elapsed time for zoom level " + str(izoom) + ": " + str(t1 - t0)) - - # Done with highest zoom level - - if make_lower_levels: - # Now loop through other zoom levels starting with highest minus 1 - - for izoom in range(zoom_range[1] - 1, zoom_range[0] - 1, -1): - if not quiet: - print("Processing zoom level " + str(izoom)) - - # Determine elapsed time - t0 = time.time() - - # Rather than interpolating the data onto tiles, we will take average of 4 tiles in higher zoom level - - zoom_path = os.path.join(path, str(izoom)) - zoom_path_higher = os.path.join(path, str(izoom + 1)) - - if index_path: - index_zoom_path = os.path.join(index_path, str(izoom)) - # List folders and turn names into integers - it0 = 1e15 - it1 = -1e15 - ix_list = [int(i) for i in os.listdir(index_zoom_path)] - ix0 = min(ix_list) - ix1 = max(ix_list) - # Now loop through the folders to get the min and max y indices - for i in range(ix0, ix1 + 1): - it_list = [ - int(j.split(".")[0]) - for j in os.listdir(os.path.join(index_zoom_path, str(i))) - ] - it0 = min(it0, min(it_list)) - it1 = max(it1, max(it_list)) - else: - ix0, it0 = deg2num(lat_range[1], lon_range[0], izoom) - ix1, it1 = deg2num(lat_range[0], lon_range[1], izoom) - - ix0 = max(0, ix0) - it0 = max(0, it0) - ix1 = min(2**izoom - 1, ix1) - it1 = min(2**izoom - 1, it1) - - # Loop in x direction - for i in range(ix0, ix1 + 1): - path_okay = False - zoom_path_i = os.path.join(zoom_path, str(i)) - - if not path_okay: - if not os.path.exists(zoom_path_i): - makedir(zoom_path_i) - path_okay = True + if not path_okay: + if not os.path.exists(zoom_path_i): + makedir(zoom_path_i) + path_okay = True + if parallel: # Loop in y direction with ThreadPool() as pool: pool.starmap( @@ -509,71 +311,20 @@ def make_topobathy_tiles( i, j, npix, - encoder, - encoder_vmin, - encoder_vmax, - compress_level, + twm, ) - for j in range(it0, it1 + 1) + for j in range(iy0, iy1 + 1) ], ) + else: + # Loop in y direction + for j in range(iy0, iy1 + 1): + # Create lower level tile + make_lower_level_tile(zoom_path_i, zoom_path_higher, i, j, npix, twm) - t1 = time.time() - - if not quiet: - print("Elapsed time for zoom level " + str(izoom) + ": " + str(t1 - t0)) - - if make_webviewer: - # Make webviewer - write_html(os.path.join(path, "index.html"), max_native_zoom=zoom_range[1]) - - if write_metadata: - if metadata is None: - metadata = {} - - metadata["longname"] = long_name - metadata["format"] = "tiled_web_map" - metadata["encoder"] = encoder - if encoder_vmin is not None: - metadata["encoder_vmin"] = encoder_vmin - if encoder_vmax is not None: - metadata["encoder_vmax"] = encoder_vmax - metadata["url"] = url - if s3_bucket is not None: - metadata["s3_bucket"] = s3_bucket - if s3_key is not None: - metadata["s3_key"] = s3_key - if s3_region is not None: - metadata["s3_region"] = s3_region - metadata["max_zoom"] = zoom_range[1] - metadata["interpolation_method"] = interpolation_method - metadata["source"] = source - metadata["coord_ref_sys_name"] = "WGS 84 / Pseudo-Mercator" - metadata["coord_ref_sys_kind"] = "projected" - metadata["vertical_reference_level"] = vertical_reference_level - metadata["vertical_units"] = vertical_units - metadata["difference_with_msl"] = difference_with_msl - metadata["available_tiles"] = True - - metadata["description"] = {} - metadata["description"]["title"] = "LONG_NAME" - metadata["description"]["institution"] = "INST_NAME" - metadata["description"]["history"] = "created by : AUTHOR" - metadata["description"]["references"] = "No reference material available" - metadata["description"]["comment"] = "none" - metadata["description"]["email"] = "Your email here" - metadata["description"]["version"] = "1.0" - metadata["description"]["terms_for_use"] = "Use as you like" - metadata["description"]["disclaimer"] = ( - "These data are made available in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE." - ) - - # Write to toml file - # toml_file = os.path.join(path, name + ".tml") - toml_file = os.path.join(path, "metadata.tml") - with open(toml_file, "w") as f: - toml.dump(metadata, f) + t1 = time.time() + print("Elapsed time for zoom level " + str(izoom) + ": " + str(t1 - t0)) def bbox_xy2latlon(x0, x1, y0, y1, crs): # Create a transformer @@ -585,31 +336,17 @@ def bbox_xy2latlon(x0, x1, y0, y1, crs): lon_max, lat_max = transformer.transform(x1, y1) return lon_min, lon_max, lat_min, lat_max +def create_highest_zoom_level_tile(zoom_path_i, i, j, izoom, twm, da, data_type, options): -def create_highest_zoom_level_tile(file_name, i, j, options): - encoder = options["encoder"] - encoder_vmin = options["encoder_vmin"] - encoder_vmax = options["encoder_vmax"] - skip_existing = options["skip_existing"] - index_path = options["index_path"] - izoom = options["izoom"] - npix = options["npix"] + file_name = os.path.join(zoom_path_i, str(j) + ".png") transformer_4326_to_3857 = options["transformer_4326_to_3857"] transformer_3857_to_crs = options["transformer_3857_to_crs"] - dem_type = options["dem_type"] - dem_list = options["dem_list"] - dataset = options["dataset"] - ds_lon = options["ds_lon"] - ds_lat = options["ds_lat"] - ds_z_parameter = options["ds_z_parameter"] - npix = options["npix"] xv = options["xv"] yv = options["yv"] dxy = options["dxy"] - bathymetry_database = options["bathymetry_database"] - interpolation_method = options["interpolation_method"] z_range = options["z_range"] - compress_level = options["compress_level"] + + skip_existing = False # Create highest zoom level tile if os.path.exists(file_name): @@ -620,18 +357,18 @@ def create_highest_zoom_level_tile(file_name, i, j, options): # Read the tile zg0 = png2elevation( file_name, - encoder=encoder, - encoder_vmin=encoder_vmin, - encoder_vmax=encoder_vmax, + encoder=twm.encoder, + encoder_vmin=twm.encoder_vmin, + encoder_vmax=twm.encoder_vmax, ) else: # Tile does not exist - zg0 = np.zeros((npix, npix)) + zg0 = np.zeros((twm.npix, twm.npix)) zg0[:] = np.nan - if index_path: + if options["index_path"]: # Only make tiles for which there is an index file - index_file_name = os.path.join(index_path, str(izoom), str(i), str(j) + ".png") + index_file_name = os.path.join(options["index_path"], str(izoom), str(i), str(j) + ".png") if not os.path.exists(index_file_name): return @@ -645,11 +382,12 @@ def create_highest_zoom_level_tile(file_name, i, j, options): x3857 = xo + xv[:] + 0.5 * dxy y3857 = yo - yv[:] - 0.5 * dxy - if dem_type == "ddb": - zg = bathymetry_database.get_bathymetry_on_grid( - x3857, y3857, CRS.from_epsg(3857), dem_list - ) - elif dem_type == "xarray": + if data_type == "ddb": + pass + # zg = bathymetry_database.get_bathymetry_on_grid( + # x3857, y3857, CRS.from_epsg(3857), dem_list + # ) + elif data_type == "xarray": # Make grid of x3857 and y3857, and convert to crs of dataset # xg, yg = np.meshgrid(x3857, y3857) xg, yg = transformer_3857_to_crs.transform(x3857, y3857) @@ -666,56 +404,56 @@ def create_highest_zoom_level_tile(file_name, i, j, options): yg_max = yg_max + dbuff # Get the indices of the dataset that are within the xg, yg range - i0 = np.where(ds_lon <= xg_min)[0] + i0 = np.where(da.x <= xg_min)[0] if len(i0) == 0: # Take first index i0 = 0 else: # Take last index i0 = i0[-1] - i1 = np.where(ds_lon >= xg_max)[0] + i1 = np.where(da.x >= xg_max)[0] if len(i1) == 0: - i1 = len(ds_lon) - 1 + i1 = len(da.x) - 1 else: i1 = i1[0] if i1 <= i0 + 1: # No data for this tile return - xd = ds_lon[i0:i1] + xd = da.x[i0:i1] - if ds_lat[0] < ds_lat[-1]: + if da.y[0] < da.y[-1]: # South to North - j0 = np.where(ds_lat <= yg_min)[0] + j0 = np.where(da.y <= yg_min)[0] if len(j0) == 0: j0 = 0 else: j0 = j0[-1] - j1 = np.where(ds_lat >= yg_max)[0] + j1 = np.where(da.y >= yg_max)[0] if len(j1) == 0: - j1 = len(ds_lat) - 1 + j1 = len(da.y) - 1 else: j1 = j1[0] if j1 <= j0 + 1: # No data for this tile return # Get the dataset within the range - yd = ds_lat[j0:j1] + yd = da.y[j0:j1] # Get number of dimensions of dataarray - if len(dataset[ds_z_parameter].shape) == 2: - zd = dataset[ds_z_parameter][j0:j1, i0:i1].to_numpy()[:] + if len(da.shape) == 2: + zd = da[j0:j1, i0:i1].to_numpy()[:] else: - zd = np.squeeze(dataset[ds_z_parameter][0, j0:j1, i0:i1].to_numpy()[:]) + zd = np.squeeze(da[0, j0:j1, i0:i1].to_numpy()[:]) else: # North to South - j0 = np.where(ds_lat <= yg_min)[0] + j0 = np.where(da.y <= yg_min)[0] if len(j0) == 0: # Use last index - j0 = len(ds_lat) - 1 + j0 = len(da.y) - 1 else: # Use first index j0 = j0[0] - j1 = np.where(ds_lat >= yg_max)[0] + j1 = np.where(da.y >= yg_max)[0] if len(j1) == 0: j1 = 0 else: @@ -724,27 +462,28 @@ def create_highest_zoom_level_tile(file_name, i, j, options): # No data for this tile return # Get the dataset within the range - yd = np.flip(ds_lat[j1:j0]) - if len(dataset[ds_z_parameter].shape) == 2: + yd = np.flip(da.y[j1:j0]) + if len(da.shape) == 2: zd = np.flip( - dataset[ds_z_parameter][j1:j0, i0:i1].to_numpy()[:], axis=0 + da[j1:j0, i0:i1].to_numpy()[:], axis=0 ) else: zd = np.squeeze( np.flip( - dataset[ds_z_parameter][0, j1:j0, i0:i1].to_numpy()[:], axis=0 + da[0, j1:j0, i0:i1].to_numpy()[:], axis=0 ) ) - zg = interp2(xd, yd, zd, xg, yg, method=interpolation_method) + zg = interp2(xd, yd, zd, xg, yg, method=options["interpolation_method"]) + + # Any value below zmin is set NaN + zg[np.where(zg < z_range[0])] = np.nan + # Any value above zmax is set NaN + zg[np.where(zg > z_range[1])] = np.nan if np.isnan(zg).all(): # only nans in this tile return - if np.nanmax(zg) < z_range[0] or np.nanmin(zg) > z_range[1]: - # all values in tile outside z_range - return - # Overwrite zg with zg0 where not nan mask = np.isnan(zg) zg[mask] = zg0[mask] @@ -753,8 +492,102 @@ def create_highest_zoom_level_tile(file_name, i, j, options): elevation2png( zg, file_name, - compress_level=compress_level, - encoder=encoder, - encoder_vmin=encoder_vmin, - encoder_vmax=encoder_vmax, + # compress_level=twm.compress_level, + encoder=twm.encoder, + encoder_vmin=twm.encoder_vmin, + encoder_vmax=twm.encoder_vmax, + ) + +def make_lower_level_tile( + zoom_path_i, + zoom_path_higher, + i, + j, + npix, + twm, +): + # Get the indices of the tiles in the higher zoom level + i00, j00 = 2 * i, 2 * j # upper left + i10, j10 = 2 * i, 2 * j + 1 # lower left + i01, j01 = 2 * i + 1, 2 * j # upper right + i11, j11 = 2 * i + 1, 2 * j + 1 # lower right + + # Create empty array of NaN to store the elevation data from the higher zoom level + zg512 = np.zeros((npix * 2, npix * 2)) + zg512[:] = np.nan + + # Create empty array of NaN of 4*npix*npix to store the 2-strid elevation data from higher zoom level + zg4 = np.zeros((4, npix, npix)) + zg4[:] = np.nan + + okay = False + + # Get the file names of the tiles in the higher zoom level + # Upper left + file_name = os.path.join(zoom_path_higher, str(i00), str(j00) + ".png") + if os.path.exists(file_name): + zgh = png2elevation( + file_name, + encoder=twm.encoder, + encoder_vmin=twm.encoder_vmin, + encoder_vmax=twm.encoder_vmax, + ) + zg512[0:npix, 0:npix] = zgh + okay = True + # Lower left + file_name = os.path.join(zoom_path_higher, str(i10), str(j10) + ".png") + if os.path.exists(file_name): + zgh = png2elevation( + file_name, + encoder=twm.encoder, + encoder_vmin=twm.encoder_vmin, + encoder_vmax=twm.encoder_vmax, + ) + zg512[npix:, 0:npix] = zgh + okay = True + # Upper right + file_name = os.path.join(zoom_path_higher, str(i01), str(j01) + ".png") + if os.path.exists(file_name): + zgh = png2elevation( + file_name, + encoder=twm.encoder, + encoder_vmin=twm.encoder_vmin, + encoder_vmax=twm.encoder_vmax, + ) + zg512[0:npix, npix:] = zgh + okay = True + # Lower right + file_name = os.path.join(zoom_path_higher, str(i11), str(j11) + ".png") + if os.path.exists(file_name): + zgh = png2elevation( + file_name, + encoder=twm.encoder, + encoder_vmin=twm.encoder_vmin, + encoder_vmax=twm.encoder_vmax, + ) + zg512[npix:, npix:] = zgh + okay = True + + if not okay: + # No tiles in higher zoom level, so continue + return + + # Compute average of 4 tiles in higher zoom level + # Data from zg512 with stride 2 + zg4[0, :, :] = zg512[0 : npix * 2 : 2, 0 : npix * 2 : 2] + zg4[1, :, :] = zg512[1 : npix * 2 : 2, 0 : npix * 2 : 2] + zg4[2, :, :] = zg512[0 : npix * 2 : 2, 1 : npix * 2 : 2] + zg4[3, :, :] = zg512[1 : npix * 2 : 2, 1 : npix * 2 : 2] + + # Compute average of 4 tiles + zg = np.nanmean(zg4, axis=0) + + # Write to terrarium png format + file_name = os.path.join(zoom_path_i, str(j) + ".png") + elevation2png( + zg, + file_name, + encoder=twm.encoder, + encoder_vmin=twm.encoder_vmin, + encoder_vmax=twm.encoder_vmax, ) diff --git a/cht_tiling/utils.py b/cht_tiling/utils.py index 434b003..287cb01 100644 --- a/cht_tiling/utils.py +++ b/cht_tiling/utils.py @@ -321,12 +321,12 @@ def png2elevation(png_file, encoder="terrarium", encoder_vmin=0.0, encoder_vmax= rgb = np.array(img.convert("RGB")).astype(float) elevation = (rgb[:, :, 0] * 256 + rgb[:, :, 1] + rgb[:, :, 2] / 256) - 32768.0 # where val is less than -32767, set to NaN - elevation[elevation < -32767.0] = np.nan + elevation[np.where(elevation < -32767.0)] = np.nan elif encoder == "terrarium16": rgb = np.array(img.convert("RGB")).astype(float) elevation = (rgb[:, :, 0] * 256 + rgb[:, :, 1]) - 32768.0 # where val is less than -32767, set to NaN - elevation[elevation < -32767.0] = np.nan + elevation[np.where(elevation < -32767.0)] = np.nan elif encoder == "uint8": rgb = np.array(img.convert("RGB")).astype(int) elevation = rgb[:, :, 0] @@ -334,11 +334,11 @@ def png2elevation(png_file, encoder="terrarium", encoder_vmin=0.0, encoder_vmax= elif encoder == "uint16": rgb = np.array(img.convert("RGB")).astype(int) elevation = rgb[:, :, 0] * 256 + rgb[:, :, 1] - elevation[elevation == 65535] = -1 + elevation[np.where(elevation == 65535)] = -1 elif encoder == "uint24": rgb = np.array(img.convert("RGB")).astype(int) elevation = rgb[:, :, 0] * 65536 + rgb[:, :, 1] * 256 + rgb[:, :, 2] - elevation[elevation == 16777215] = -1 + elevation[np.where(elevation == 16777215)] = -1 elif encoder == "uint32": rgb = np.array(img.convert("RGBA")).astype(int) elevation = ( @@ -347,22 +347,22 @@ def png2elevation(png_file, encoder="terrarium", encoder_vmin=0.0, encoder_vmax= + rgb[:, :, 2] * 256 + rgb[:, :, 3] ) - elevation[elevation == 4294967295] = -1 + elevation[np.where(elevation == 4294967295)] = -1 elif encoder == "float8": rgb = np.array(img.convert("RGB")).astype(float) i = rgb[:, :, 0] elevation = encoder_vmin + (encoder_vmax - encoder_vmin) * i / 254 - elevation[i == 0] = np.nan + elevation[np.where(i == 0)] = np.nan elif encoder == "float16": rgb = np.array(img.convert("RGB")).astype(float) i = rgb[:, :, 0] * 256 + rgb[:, :, 1] elevation = encoder_vmin + (encoder_vmax - encoder_vmin) * i / 65534 - elevation[i == 0] = np.nan + elevation[np.where(i == 0)] = np.nan elif encoder == "float24": rgb = np.array(img.convert("RGB")).astype(float) i = rgb[:, :, 0] * 65536 + rgb[:, :, 1] * 256 + rgb[:, :, 2] elevation = encoder_vmin + (encoder_vmax - encoder_vmin) * i / 16777214 - elevation[i == 0] = np.nan + elevation[np.where(i == 0)] = np.nan elif encoder == "float32": rgb = np.array(img.convert("RGBA")).astype(float) i = ( @@ -372,7 +372,7 @@ def png2elevation(png_file, encoder="terrarium", encoder_vmin=0.0, encoder_vmax= + rgb[:, :, 3] ) elevation = encoder_vmin + (encoder_vmax - encoder_vmin) * i / 4294967294 - elevation[i == 0] = np.nan + elevation[np.where(i == 0)] = np.nan return elevation