ColumnDataSource - Image

[lcboutin]

Hi,

I’m trying to display a series of geo TIFF images, which would be components of a ColumnDataSource. I am struggling on how to pass the source so that it can retrieve the right image.

Step 1) Create several rasters (TIF) from my original geopandas dataframe
Step 2) Create ColumnDataSource
Step 3) Show different rasters using a slider *Need the CDS and CDSView to work and pass the information to image.

##################################################

import pandas as pd
import geopandas as gpd
import rasterio
from geocube.api.core import make_geocube
import os


from bokeh.plotting import figure, show, output_file, save
from bokeh.models import ColumnDataSource, HoverTool, ColorMapper, ColorBar,Ticker,OpenURL, TapTool, Circle, LogTicker, LabelSet,FixedTicker, Button, Dropdown, CheckboxButtonGroup, RadioButtonGroup, GeoJSONDataSource, Patches 
from bokeh.models import widgets, CustomJS, Slider, Label, LinearColorMapper, LogColorMapper, ContinuousTicker, CategoricalColorMapper, axes,CDSView,GroupFilter, BooleanFilter
from bokeh.models import Range1d, CDSView, Div, Paragraph, LogAxis, Select
from bokeh.models import MultiLine, Scatter, Rect, Line, LegendItem, Legend, PrintfTickFormatter, CheckboxGroup, CustomJSFilter
from bokeh.layouts import gridplot
from bokeh.models.glyphs import Quad
from bokeh.layouts import layout, column, row, gridplot
from bokeh.tile_providers import get_provider, Vendors
from bokeh.palettes import Spectral8, magma
from bokeh.transform import linear_cmap, factor_cmap, factor_mark
from bokeh.models import ColorBar

#import pyodbc
import bokeh as bk

import numpy as np

####################################################
#%%   STEP 1  - Data Generation
####################################################

#Note: rasters don't have the same size
data = [[1,1,1,2],[1,2,1,2],[1,2,2,2],[1,1,2,2],\
        [2,1,1,4],[2,2,1,4],[2,2,2,4],[2,1,2,4],[2,1,3,4],[2,3,3,4]]

df = pd.DataFrame(data,columns=['rasterID','X','Y','col'])

adf = df[df['rasterID']==1]
gdf1 = gpd.GeoDataFrame(adf, geometry=gpd.points_from_xy(adf.X, adf.Y), crs="EPSG:3857")
out_grid = make_geocube(
                vector_data=gdf1,
                measurements=['col'],
                resolution=(-0.5, 0.5))
#Save my TIFF1
out_grid['col'].rio.to_raster("Raster1.TIF")
    
bdf = df[df['rasterID']==2]
gdf2 = gpd.GeoDataFrame(bdf, geometry=gpd.points_from_xy(bdf.X, bdf.Y), crs="EPSG:3857")
out_grid = make_geocube(
                vector_data=gdf2,
                measurements=['col'],
                resolution=(-0.5, 0.5),
            )
#Save my TIFF2
out_grid['col'].rio.to_raster("Raster2.TIF")

####################################################
#%%   STEP 2 - Data Source
####################################################
DataSource = []
raster1 = rasterio.open('Raster1.TIF')
raster2 = rasterio.open('Raster2.TIF')
DataSource.append(['1', raster1])
DataSource.append(['2', raster2])
df_DS = pd.DataFrame(DataSource,columns=['ID','raster'])

####################################################
#%%  STEP 3 - Bokeh
####################################################

pv = figure(tools=['wheel_zoom','box_zoom','reset','pan','lasso_select','box_select']
            ,match_aspect=True,x_axis_type='mercator',y_axis_type='mercator'
            ,width=1000,height=700)

#Specified the map extent to start with.
pv.x_range = Range1d(start=1,end=5)
pv.y_range = Range1d(start=1,end=5)


#This is where I would like to use the DS_source and DS_view, but can't figure it out. Just ploting first raster instead.
#DS_source = ColumnDataSource(data=df_DS)
#DS_view = CDSView(source=DS_source,filters=[GroupFilter(column_name='ID',group='1')])
info = df_DS.iloc[1]['raster']
DS_rend = pv.image(image = [np.flipud(info.read(1))],
                       x=info.bounds[0],
                       y=info.bounds[1],
                       dw=info.bounds[2]-info.bounds[0],
                       dh=info.bounds[3]-info.bounds[1],
                       legend_label="Value",
                       palette = 'RdYlGn11',
                       alpha=0.5
                      )
DS_hover = HoverTool(renderers=[DS_rend],tooltips=[('Value','@image{0.0}')])
pv.add_tools(DS_hover)
pv.legend.click_policy='hide'

#Simple slider to change the raster
slider = Slider(start=0, end=1, width=400
                  , value=0, step=1, title="ID: Raster1"
                  , sizing_mode="fixed")
slider.show_value = True
slider.tooltips = True
slider_dict = {'0':'Raster1', '1':'Raster2'}
cb = bk.models.CustomJS(args=dict(
        slider=slider
        ,slider_dict=slider_dict
        #,DS_source=DS_source
        ),
    code="""
    console.log(cb_obj.value);
    //SLIDER
    var id = slider.value;
    slider.title = 'ID: ' + slider_dict[id]
    //DS_source.filters[0].group = slider_dict[id]
    """)
slider.js_on_change('value',cb)
lo = layout([[slider]
            ,[pv]]
            ,sizing_mode="fixed")
save(lo,r'Bokeh_discourse.20240609.html')

[Bryan]

Hi @lcboutin please edit your post to use code formatting so that the code is intelligible (either with the </> icon on the editing toolbar, or triple backtick ``` fences around the code blocks)

[lcboutin]

Corrected now.

[gmerritt123]

Do your rasters all have the same dimensions (i.e. x, y, dw, dh)? There are a couple ways I can think of to structure your input but what’s optimal probably depends on this.

[lcboutin]

Hi Gaelen,

Unfortunately, they all have different size, since the “plume extent” varies overtime. Hence, why I thin I can’t only store the np.array of the raster, but benefit from storing the meta data from the raster.

[gmerritt123]

Right, so there’s a bit more data prep you need to (you’re kinda already half doing it I think but tough to tell exactly as I’m not a rasterio expert).

Ultimately, to get it “bokeh ready”, you need to extract the 2D arrays from each raster along with their georeferencing info (i.e. w, h, dw, dh). Warning: you’ll need to pre-warp your rasters into ‘epsg:3857’ (web merc) to be able to plot with esri basemap tiles etc

Then with that all prepped, you can assemble a ColumnDataSource that holds all that info and you tell the image renderer to plot using those respective fields. Next, to display only one raster at a time (i.e. based on a slider value), you can use IndexFilter and CDSView to show only one at a time. This makes the CustomJS super super minimal

Commented/functional example:

from bokeh.plotting import figure, show
from bokeh.layouts import column
from bokeh.models import ColumnDataSource, Image, LinearColorMapper, CustomJS, Slider, ColorBar, IndexFilter, CDSView
import numpy as np

x = np.linspace(0, 50, 300)
y = np.linspace(0, 50, 300)
xx, yy = np.meshgrid(x, y)
#make a bunch of 2d arrays (equivalent to the data in your tifs)
d1 = xx**2 * yy**2 #say this is the plume at t = 0
d2 = xx * yy**2 #say this is the plume at t = 1
d3 = xx**2 * yy**0.5 #say this is the plume at t = 2

#setup a data structure containing all the rasters and their coord info
#you'll have to extract this info from the rasterio objects (it's GetGeoTransform etc in gdal but I don't know rasterio well)
src_dict={'x':[0,20,30], #lower left coord for d1, d2, d3 
          'y':[10,15,30], #ditto
          'dw':[25,30,50], #width of d1,d2,d3
          'dh':[5,10,15],  #height
          'array':[d1,d2,d3] #this is a 3D "stack" of the raster data
          }

src = ColumnDataSource(data=src_dict)

#create a filter that will display only one at a time
filt = IndexFilter(indices=[0]) #start w d1 (i.e. only index 0)


p = figure(width=800,height=600,match_aspect=True,x_range=[0,50],y_range=[0,50])

# create the image renderer to follow not only the 'array' field for its data, but also the x,y, dh, and dw fields
r= p.image(image='array',x='x',y='y',dw = 'dw',dh='dh'
       ,source=src, view = CDSView(filter=filt) #assign the renderer's view to follow the index filter  
       ,palette='Sunset11',alpha=0.5)

r.glyph.color_mapper.low = 0
r.glyph.color_mapper.high = 1200
cbar = ColorBar(color_mapper=r.glyph.color_mapper)
p.add_layout(cbar,'right')


sl = Slider(start=0,end=len(src_dict['x'])-1,value=0,step=1)

cb = CustomJS(args=dict(sl=sl,filt=filt,src=src)
             ,code='''
             filt.indices = [sl.value] //set the filter indices based on the slider value             
             src.change.emit() //triggers a redraw with the updated filter
             ''')
sl.js_on_change('value',cb)

show(column([p,sl]))

[lcboutin]

OK, let me try this and see if it will work. Thanks, really appreciate the help. I’ll reformat the data dictionary to align with your example.

[lcboutin]

Well, it worked like a charm! I reduced my HTML export from 360MB to 45MB using rasters instead of scatter. Thanks a bunch Gaelen, that is fantastic!!