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

We initialize a FastAPI application to serve our tile service.

Step 2: Load GeoTIFF Data

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

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():
   • This step resizes the raster data to match the desired shape of the tile, typically 256×256 pixels.
   • It ensures that the raster data covers the exact spatial extent of the tile.
   • Interpolation helps maintain the resolution and accuracy of the raster data when resizing to fit the tile shape.
2. Resizing based on Tile Extent using rio.pad_box()
   • After interpolation, additional resizing is performed to ensure that the tile covers the specified geographic extent.
   • This step adjusts the size of the tile to match its bounding box coordinates, ensuring accurate representation of the geographic area.
   • Resizing based on the tile extent helps generate map tiles that accurately represent the specified area while maintaining consistency in tile dimensions.
   • This resizing prevents the tile distortion

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.