Datashader is a python library for fast rasterization and visualization of larger datasets preserving the data’s distribution. Datashader is built using Numba which a open source JIT compiler that translates subsets of python and numpy into faster machine code.
Before generating the heatmap, we will set up conda environment and install required python GIS packages for this task. Make sure conda is installed on your system. Use following commands to create a conda environment and to install python libraries. When you install geocube library using conda, it will automatically install related dependencies that includes gdal, shapely, rasterio, geopandas, xarray, rioxarray etc. After that, will install datashader.
And we will also install pygeos library which is used to speed up the vectorized operations in GeoPandas and Shapely.
|
1 2 3 4 5 |
(base) geoknight@pop-os:~$conda create -n spatial-dev.guru python=3.10 (base) geoknight@pop-os:~$conda activate spatial-dev.guru (spatial-dev.guru) geoknight@pop-os:~$conda install -c conda-forge geocube (spatial-dev.guru) geoknight@pop-os:~$conda install -c conda-forge datashader (spatial-dev.guru) geoknight@pop-os:~$conda install -c conda-forge pygeos |
In this tutorial we will create heatmap for uk accidents dataset and then generate the tiles and serve them in openlayers map. We will divide this tutorial in three steps.
1. Generate Heatmap using Datashader
2. Create tiles for Heatmap using render_tiles method of Datashader
3. Serve the heatmap tiles in OpenLayers
To get full source code for this tutorial, click here.
To check the live demo for this tutorial, click here.
1. Generate Heatmap using Datashader
We will create heatmap for uk_accidents dataset which will show the density of accidents. uk_accidents data set is in csv format. We have three csv dataset for uk_accidents for different years. Each csv file have latitude and longitude coordinates of each accident. We are reading these csvs using pandas and merging them into single pandas dataframe.
Using geopandas, we will convert the lat long of dataframe into geodataframe. This geodataframe is in 4326 Geographic coordinate system. We will reproject it to 3857 projected coordinate system. See line number from 16 to 30.
At line number 30, we have dataframe containing the projected longitude(x) and latitude(y) values. We will use these coordinates to generate the heatmap.
To generate heatmap for the coordinates, we will use datashader.
1. At line number 32, we are creating an empty canvas using datashader. Canvas is like a plane where you will project your coordinates.
2. At line number 33, we are adding projected x and y values to canvas. This will return a rasterized xarray object.
3. At line number 34, we are using shade function from datashader to aggregate and rasterize the the points. I am using ‘log’ interpolation method to aggregate the points. You can use other interpolation method for testing.
4. At line number 35, we are using spread function from datashader. spread is used to make sparse heatmap more visible.
That’s it, you have successfully generated heatmap. You can run the code in jupyter notebook for better visualization.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 |
import colorcet import datashader as ds, pandas as pd, geopandas as gpd # Read UK Accidents Point csvs uk_accidents_1 = pd.read_csv('csvs/uk_accidents_2005_to_2007.tar.xz', low_memory=False) uk_accidents_2 = pd.read_csv('csvs/uk_accidents_2009_to_2011.tar.xz', low_memory=False) uk_accidents_3 = pd.read_csv('csvs/uk_accidents_2012_to_2014.tar.xz', low_memory=False) uk_accidents = pd.concat([uk_accidents_1, uk_accidents_2, uk_accidents_3]) # Convert Long Lat into numeric type uk_accidents['Longitude'] = pd.to_numeric(uk_accidents['Longitude']) uk_accidents['Latitude'] = pd.to_numeric(uk_accidents['Latitude']) # Convert Long Lat into Point Geometry uk_accidents = gpd.GeoDataFrame(geometry = gpd.points_from_xy(x=uk_accidents['Longitude'], y=uk_accidents['Latitude'])) uk_accidents = uk_accidents.set_crs('EPSG:4326') # Reprojecting to 3857 coordinate system uk_accidents = uk_accidents.to_crs('EPSG:3857') uk_accidents['value'] = 1 uk_accidents = uk_accidents[uk_accidents.is_valid] uk_accidents = uk_accidents[~uk_accidents.is_empty] # Get x, y coordinates uk_accidents['x'] = uk_accidents.geometry.x uk_accidents['y'] = uk_accidents.geometry.y df = uk_accidents[['x', 'y']] # Generate Heatmap cvs = ds.Canvas(plot_width=1000, plot_height=1000) agg = cvs.points(df, 'x', 'y') img = ds.tf.shade(agg, cmap=colorcet.fire, how='log') img = ds.tf.spread(img, px=1, shape='circle', how='add', mask=None, name=None) |
2. Create tiles for Heatmap using render_tiles method of Datashader
In the first step, we have only visualized heatmap for uk_accidents dataset. If you want to create the tiles for the heatmap, you have to use render_tiles method. This method, will generate the heatmap tiles in 3857 projected coordinate system. That was the reason, we first projected the uk_accidents coordinates from 4326 to 3857.
In the render_tiles method at line number 71, we pass many parameters to generate tiles.
1. First parameter is map_extent. This parameter tells the datashader for which map extent it will create the tiles.
2. levels parameter takes a range object which specifies the range of zoom levels. The datashader will create tiles for those zoom levels in range object.
3. output_path define the storage path for created tiles.
4. load_data_func parameter load the dataframe from function. See line number 39.
5. rasterize_func parameter takes a function argument that returns rasterized xarray object. See line number 43.
6. shader_func parameter takes a function argument which applies interpolation process on the dataset to generate heatmap. See line number 50.
7. post_render method takes a function argument. This function generate tiles images. See line number 59.
Running render_tiles will create the heatmap tiles and store them in tiles directory. You can specify any output_path. In our case it is tiles directory where generated tiles will be stored. tiles directory will look like below structure. In our code, we specified 11 zoom levels for tile creation. That is why we have 11 directory naming from 0 to 10 which contains tiles images for respective zoom levels.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 |
import pandas as pd from datashader.tiles import render_tiles import geopandas as gpd import datashader as ds, pandas as pd, colorcet import datashader.transfer_functions as tf from PIL import ImageDraw def read_uk_accidents(): # Read UK Accidents Point csvs uk_accidents_1 = pd.read_csv('csvs/uk_accidents_2005_to_2007.tar.xz', low_memory=False) uk_accidents_2 = pd.read_csv('csvs/uk_accidents_2009_to_2011.tar.xz', low_memory=False) uk_accidents_3 = pd.read_csv('csvs/uk_accidents_2012_to_2014.tar.xz', low_memory=False) uk_accidents = pd.concat([uk_accidents_1, uk_accidents_2, uk_accidents_3]) # Convert Long Lat into numeric type uk_accidents['Longitude'] = pd.to_numeric(uk_accidents['Longitude']) uk_accidents['Latitude'] = pd.to_numeric(uk_accidents['Latitude']) # Convert Long Lat into Point Geometry uk_accidents = gpd.GeoDataFrame(geometry = gpd.points_from_xy(x=uk_accidents['Longitude'], y=uk_accidents['Latitude'])) uk_accidents = uk_accidents.set_crs('EPSG:4326') # Reprojecting to 3857 coordinate system uk_accidents = uk_accidents.to_crs('EPSG:3857') uk_accidents['value'] = 1 uk_accidents = uk_accidents[uk_accidents.is_valid] uk_accidents = uk_accidents[~uk_accidents.is_empty] # Get x, y coordinates uk_accidents['x'] = uk_accidents.geometry.x uk_accidents['y'] = uk_accidents.geometry.y df = uk_accidents[['x', 'y']] return [df, tuple(uk_accidents.total_bounds)] def load_data_func(x_range, y_range): global df return df def rasterize_func(df, x_range, y_range, height, width): # aggregate cvs = ds.Canvas(x_range=x_range, y_range=y_range, plot_height=height, plot_width=width) agg = cvs.points(df, 'x', 'y') return agg def shader_func(agg, span=None): # shader func img = tf.shade(agg, cmap=colorcet.fire, how='log') img = tf.spread(img, px=1, shape='circle', how='add', mask=None, name=None) img = tf.set_background(img, None) return img def post_render_func(img, **kwargs): # Create tiles draw = ImageDraw.Draw(img) draw.text((5, 5), '', fill='rgb(255, 255, 255)') return img if __name__ == "__main__": global df df, map_extent = read_uk_accidents() render_tiles(map_extent, levels=range(11), output_path='tiles', load_data_func=load_data_func, rasterize_func=rasterize_func, shader_func=shader_func, post_render_func=post_render_func) |
3. Serve the generated heatmap tiles in OpenLayers.
In previous step, we created heatmap tiles for 11 zoom levels and stored them in tiles directory. We will use openlayers to render them on map.
In openlayers, we will use ol.source.XYZ with ol.layer.Tile to serve tiles. See line number 35-42. At line number 37-38, we need to pass the url path of tiles that we created.
In our case, it is tile directory, so the url path for this will be /tiles/{z}/{x}/{y}.png. Our html file and tiles directory are in same location. Run html file where our openlayers code resides using any server. Then load the file in web browser. It will render the map and load the tiles from tile directory in your map.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 |
<html> <head> <title>OpenLayers zoom in and zoom out by box selection</title> <script src="dependencies/jsts.js"></script> <script src="dependencies/ol.js"></script> <link rel="stylesheet" href="dependencies/ol.css"> <style> .ol-dragzoom { background-color: rgba(255,255,255,0.4); border-color: rgb(0, 0, 0); border-width: 2; } html, body { height: 100%; margin: 0; padding: 0; } </style> </head> <body> <div id = "map"> <div id = "popup"></div> </div> </body> <script> //Define Map var map = new ol.Map({ target: 'map', layers: [ new ol.layer.Tile({ source: new ol.source.OSM() }), new ol.layer.Tile({ source: new ol.source.XYZ({ url: '/tiles/{z}/{x}/{y}.png', }), minZoom: 1, // maxZoom: 10 }) ], view: new ol.View({ center: ol.proj.fromLonLat([-0.1275, 51.507222]), zoom: 5 }) }); </script> </html> |
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 email at contact@spatial-dev.guru.
Cool! BTW, that colormap is best for a dark map layer, so I’d suggest switching to such a map or switching to a colormap that starts from a lighter color if you want to stay with that map layer. Right now if there are any isolated peaks the user won’t notice them as they will be drawn in yellow on top of a similarly colored map.