Create your own layers in JupyterGIS
Introduction
JupyterGIS is a JupyterLab extension that makes it possible to work with geospatial data. Users can interact with a map directly from the JupyterLab UI, or through a notebook. For instance, one can open a jGIS file from the file browser and it will show the layers that it contains, just like in QGIS. But the same file can also be opened programmatically in a notebook, and changes to the map will be synchronized between the two views. That allows for real-time collaboration and it is great for visualizing globally accessible data, such as pre-generated raster data from OpenStreetMap or satellite data from NASA.
Research is done in notebooks
But what if you are more of a scientist who processes and analyses data? A typical workflow usually consists of a notebook where you load raw data, do all the science, and show the results on a map. The latter often consists of specific parts carefully chosen to illustrate your results, because you cannot show it all. It is good if you want to share your findings with people who are not as versed into the scientific side as you are, but it also comes with drawbacks. Maybe they don’t want to open the notebook in the first place, with all the details of the intermediary steps, and instead interact with the final map. And you’re giving them the opposite: the complex details of your research, without the ability to explore the results on their own. They fall into the category of people who use JupyterLab as an application to visualize rich documents, while you use it to do research with notebooks. What’s great about JupyterGIS is that it acts like a bridge between both kinds of users: it allows them to collaborate on a map either from JupyterLab’s UI or from a notebook.
The idea is to use a notebook to compute the raster tiles on-the-fly, as they are requested by a user exploring a map outside the notebook. As far as the frontend is concerned (the code that shows the map in the browser), tiles are just provided by a server at a specific URL pattern, such as https://tile.openstreetmap.org/{z}/{x}/{y}.png. When working locally, we can run a server to serve tiles at e.g. http://127.0.0.1:8000/{z}/{x}/{y}.png in our notebook, which acts as a tile provider, and point the frontend to that. But there is still a missing piece in our data flow, which consists of going from geographical data, represented by a 2D array of numerical values, to PNG images ready to be sent over the wire. And there happens to be a great Python library for that, called TiTiler.
The notebook as a tile server
TiTiler is developed by Development Seed and is described as “a modern dynamic tile server built on top of FastAPI and Rasterio/GDAL”. There are two parts in it: Rasterio which does all the heavy lifting for generating tiles, and FastAPI for serving tiles through a web API.
- Rasterio can clip data to a tile area, reproject data to the projection used by the frontend (such as Web Mercator), resample, apply colormaps, process data with (custom) algorithms, etc. TiTiler also supports extensions, and one that is very interesting for scientific computing is titiler.xarray. As its name implies, it allows working with data through the Xarray library. Xarray is very popular for processing (big) multidimensional data, thanks to its integration with Dask and Zarr. If you work with satellite and earth observation data, you must be familiar with these libraries. In fact, the Pangeo project is built on top of this software stack.
- FastAPI is a very popular Python framework for building web server applications. In TiTiler, it is used to serve the tiles requested by the frontend but also to access the underlying application. For instance, the frontend could make a request to get the list of supported colormaps. Note that since FastAPI supports OpenAPI automatic UI documentation, you could use it as a frontend to query TiTiler.
JupyterGIS-tiler
JupyterGIS-tiler is a JupyterGIS extension that lets you use a notebook to process your (Xarray) data and make it available as raster tiles, so that it can be visualized in JupyterGIS. From the frontend’s point of view it is a regular raster source, so there is nothing special about it except that:
- it is only available as long as the notebook is running. Once you close the notebook (and in particular the associated kernel), the tiles won’t be served anymore and the frontend will get a 404 on tile requests.
- it is served from your local host, which is fine if you run JupyterGIS locally but problematic when deploying on the Internet. In that case, users wouldn’t be able to collaborate on your layer. So rather than letting users deal with reverse proxies and such, we take advantage of the fact that JupyterLab also has a web server, and we proxy the TiTiler server through the Jupyter server. This guarantees that if you deployed JupyterLab, then JupyterGIS-tiler is deployed for free.
Example: flow accumulation
A good example of what JupyterGIS-tiler is able to do is illustrated in the south_america.ipynb notebook. This notebook starts by downloading digital elevation and flow accumulation data from HydroSHEDS.
Flow accumulation is a product derived from flow direction, itself derived from a digital elevation model (DEM). Flow accumulation can represent the river network, as illustrated above. The specificity of this data is that it is encoded at the pixel level. Unlike a DEM where there is some kind of continuity between surrounding pixels, flow accumulation can be very high on a single pixel, and close to zero next to it. This makes it hard to visualize because naive downsampling won’t yield good results. Instead, data needs to be reprojected by taking the maximum value of surrounding pixels. Also, one needs to thicken the single-pixel flow accumulation values to get a nice representation of rivers at all zoom levels, with colors reflecting flow intensity.
Fortunately one can make use of GDAL and Python algorithms to process the data the way we want. GDAL will mainly be useful at the reprojection level, while Python can provide additional processing from libraries such as SciPy. Here is the custom algorithm we use to thicken flow accumulation and get a nice visual river network:
import numpy as np
from titiler.core.algorithm import BaseAlgorithm
from rio_tiler.models import ImageData
import scipy.ndimage
radius = 1
circle = np.zeros((2*radius+1, 2*radius+1)).astype("uint8")
y, x = np.ogrid[-radius:radius+1, -radius:radius+1]
index = x**2 + y**2 <= radius**2
circle[index] = 1
class AlgoAcc(BaseAlgorithm):
def __call__(self, img: ImageData) -> ImageData:
data = np.log(img.data[0])
data = scipy.ndimage.maximum_filter(data, footprint=circle)
return ImageData(
data,
img.mask,
assets=img.assets,
crs=img.crs,
bounds=img.bounds,
)This algorithm can then be passed to TiTiler and used for creating tiles on the fly:
from jupytergis.tiler import GISDocument
doc = GISDocument()
await doc.add_tiler_layer(
name="Flow accumulation layer",
data_array=da_acc,
colormap_name="plasma",
rescale=np.log([vmin_acc, vmax_acc]),
algorithm=AlgoAcc,
reproject="max",
)Here is the final map with the flow accumulation layer overlayed on top of the terrain layer, where tiles are computed on the fly from our Xarray data:
Conclusion
This opens up new possibilities by connecting JupyterGIS’s interactive and collaborative visualization system with the power of the Python ecosystem for data processing. Projects like Pangeo could benefit from JupyterGIS-tiler by plugging their infrastructure to access big data and provide compute capabilities, enabling an open-source alternative to Google Earth Engine and other platforms.
About the author
David Brochart is a technical director at QuantStack. He contributes to the real-time collaboration system in JupyterLab, and is interested in big-geo-data.
Acknowledgments
The initial development of JupyterGIS at QuantStack was funded by the European Space Agency (ESA). Today, the JupyterGIS development team comprises other significant stakeholders, with significant contributions from the Eric and Wendy Schmidt Center for Data Science & Environment (DSE) at UC Berkeley.
The work by David Brochart on integrating JupyterGIS with the Xarray stack was funded by the Centre National d’Études Spatiales (CNES).
