Creating a GeoTIFF raster XYZ tile service in python with caching capability

In this tutorial, we’ll create a tile service for GeoTIFF files using Python, FastAPI, and various libraries such as xarray, rioxarray, rasterio, and Pillow. This service will dynamically generate map tiles on-the-fly and cache them for efficient retrieval.

Prerequisites

Before starting, ensure you have the following installed:

• Python (3.7 or later)
• FastAPI (pip install fastapi)
• Uvicorn (pip install uvicorn)
• xarray (pip install xarray)
• rioxarray (pip install rioxarray)
• shapely (pip install shapely)
• rasterio (pip install rasterio)
• Pillow (pip install pillow)

Overview

1. Setup FastAPI: Initialize a FastAPI application.
2. Load GeoTIFF Data: Load the GeoTIFF file using xarray and rioxarray.
3. Calculate Tile Extent: Define a function to calculate the geographic extent of a tile.
4. Define Tile Generation Functions: Implement functions to generate map tiles, interpolation etc.
5. Serving and Caching Tiles: Create an API endpoint to serve map tiles.
6. Runnig TileServer
7. Rendering raster tiles in OpenLayers map

To get full source code for this tutorial, click here.

To check the live demo for this tutorial, click here.

Let’s dive into the implementation.

Step 1: Setup FastAPI

from fastapi import FastAPI, Response
from starlette.responses import FileResponse

app = FastAPI()

We initialize a FastAPI application to serve our tile service.

Step 2: Load GeoTIFF Data

import xarray as xr
import rioxarray as riox

# Load GeoTIFF file
ds = xr.open_dataset('raster/dem_hp.tif')
ds = ds.rio.reproject("EPSG:3857")  # Reproject to Web Mercator
ds = ds.band_data  # Extract band data
ds = ds.astype(np.int32)  # Convert data type to int32

Here, we load the GeoTIFF file using xarray and reproject it to the Web Mercator coordinate system (EPSG:3857) using rioxarray. We then extract band data and convert it to int32 for efficient processing.

Step 3: Calculate Tile Extent

We define a function ST_TileEnvelope to calculate the geographic extent of a tile based on its zoom level, x, and y coordinates. This function calculates the bounding box coordinates of the tile in the Web Mercator projection. This function is directly inspired from postgis function ST_TileEnvelope

Step 4: Define Tile Generation Functions

def normalize_array(arr, min_value, max_value):
    # Function to normalize array values
    # ...

def convert_black_to_transparent(image):
    # Function to convert black pixels to transparent
    # ...

def generate_tile(zoom, x, y, bounds, margin, min_value, max_value):
    tile_envelope = ST_TileEnvelope(zoom, x, y, bounds, margin)

    xmin, ymin = tile_envelope[0]
    xmax, ymax = tile_envelope[2]

    try:
        data = ds.rio.clip_box(minx=xmin, miny=ymin,
                               maxx=xmax, maxy=ymax, auto_expand=True)
        data = data.rio.reproject(data.rio.crs, shape=(256, 256))
        data = data.rio.pad_box(xmin, ymin, xmax, ymax, 0)
    except riox.exceptions.NoDataInBounds:
        # Return empty tile if no data found in bounds
        return None

    arr = normalize_array(data.values[0], min_value, max_value)
    arr = arr.astype(np.int32)
    img = Image.fromarray(arr)

    img = img.convert('RGBA')

    img = convert_black_to_transparent(img)

    img_bytes = io.BytesIO()
    img.save(img_bytes, format='PNG')
    img_bytes.seek(0)

    return img_bytes.getvalue()

We define helper functions to facilitate tile generation. normalize_array normalizes array values to the 0-255 range. convert_black_to_transparent converts black pixels to transparent in an image. generate_tile generates a map tile based on zoom level, tile coordinates, and other parameters.
This step is basically interpolation of raster data based on extent.

The interpolation and resizing steps in the tile generation process involve two main components:

1. Interpolation using rio.reproject():
2. Resizing based on Tile Extent using rio.pad_box()

Step 5: Serving and Caching Tiles

Here, we define an API endpoint /tiles/{z}/{x}/{y}.png to serve map tiles. When a request is made for a tile at a specific zoom level (z), x-coordinate (x), and y-coordinate (y), the server checks if the tile is available in the cache directory. If not, it generates the tile using the generate_tile function and saves it to the cache directory. Finally, it returns the tile as a PNG image.

Cache Directory structure looks like below where path of tile is maintained based on {z}/{x}/{y}.png

Step 6: Running the Application

To run the tile service, execute the following command:

Step 7: Rendering Raster Tiles in OpenLayers Map

The provided HTML code sets up an OpenLayers map with two layers:

1. OSM Layer: This layer displays a base map using OpenStreetMap tiles.
2. Custom Raster Tile Layer: This layer is configured to display raster tiles using the specified URL pattern (tiles/{z}/{x}/{y}.png), where {z}, {x}, and {y} represent the zoom level, tile’s column index, and tile’s row index respectively.

The JavaScript code initializes the map, sets its view to a specific center and zoom level, and adds the OSM and custom raster tile layers to the map. The raster tile layer will dynamically fetch tiles based on the provided URL pattern and display them on the map.

I hope this tutorial will create a good foundation for you. If you want tutorials on another GIS topic or you have any queries, please send an mail at contact@spatial-dev.guru.

We also offer freelancing services. Please email us at contact@spatial-dev.guru for any query.