Hi forum,
I am working with large raster files (each one about 200k by 200k). I have mainly been using hvPlot to create interactive plots for certain small regions in my rasters i.e. I take a coordinate, extract about a 1k by 1k window from the raster and plot that in hvPlot. This has been great to visualise the rasters in specific regions. However, what would be perfect would be being able to visualise the entire raster in one go without having to create separate plots for each region.
Solutions in mind
-
The plot initially shows the base tiles for the entire region covered by the raster but no data is plotted yet. If you zoom in enough, then the raster is plotted for the area that you have zoomed into. The threshold would be set so that plotting the raster for that area should be manageable. Once zoomed in, if you pan around the plot of the raster should update for that region. Though, if you zoom out past the threshold then the plot of the raster disappears showing only the tiles.
-
(Preferred one) Initially the plot shows the entire region and plots an aggregated version of the original raster (I can create this separately). Then if zoom in enough it shows the original detailed raster. This is my preferred solution, though if I can get the first solution working then I think I can work out how to do this.
Overall, my goal is just to be able to visualise the entire raster in a single plot. So I am not fixed on the above, happy to go with any other solutions which achieves this.
My thoughts
I am quite new to the HoloViz libraries, so admittedly my knowledge is quite basic. From a read of the various user guides, it seems like I would need to use HoloViews and a RangeXY stream.
I saw this post which shows how a layer can be showed when zoomed in sufficiently. This seems to be very close to what I am after. Though, it seems to me the solution in that post requires to first create a plot of the entire raster (the point plot in that post) and then the RangeXY stream can be applied to filter and turn the layer on/off as required. That would not be possible in my case as creating the initial entire plot of the raster would require loading the entire raster which won’t fit in my RAM.
Below is some reproducible data and some code that I have written. I know this code doesn’t work (additionally, it creates a plot of the entire raster which is not feasible with my actual data), but hopefully it demonstrates the general idea that I am after.
# create dummy raster - in reality this would be saved on disk not completely loaded into memory
data = np.random.rand(1000, 1000)
lons = np.linspace(-100, -80, data.shape[0])
lats = np.linspace(30, 45, data.shape[1])
da = xr.DataArray(
data,
coords={'y': lats, 'x': lons},
dims=['y', 'x']
)
da = da.rio.write_crs("EPSG:4326")
# reproject to Web Mercator so can be overlayed with tiles
da = da.rio.reproject("EPSG:3857")
# create base tiles for entire region covered by the raster
tiles_osm = hv.element.tiles.OSM().opts(
opts.Tiles(width=500, height=500,
xlim=(lnglat_to_meters(-105, 30)[0], lnglat_to_meters(-75, 30)[0]),
ylim=(lnglat_to_meters(-105, 30)[1], lnglat_to_meters(-75, 45)[1]))
)
# function to only show the raster when zoomed in enough
def filterRast(plot, x_range, y_range):
if x_range is None or y_range is None:
return tiles_osm # don't think this right, not sure what should be returned if not zoomed in enough
elif (x_range[1] - x_range[0] > 2) | (y_range[1] - y_range[0] > 2):
return tiles_osm # don't think this right, not sure what should be returned if not zoomed in enough
else:
return plot[x_range, y_range]
# stream to capture the rendered range of the plot
range_stream = hv.streams.RangeXY(source=rast_plot)
streams=[range_stream]
# plot of entire raster - this won't be possible with my actual data
rast_plot = hv.QuadMesh(da, ['x', 'y']).opts(
opts.QuadMesh(tools=['hover']))
# apply the filter function to the raster plot
filter_rast_plot = rast_plot.apply(filterRast, streams=streams)
tiles_osm * filter_rast_plot
Also, if it is of any use, below is a function that extracts a window around a coordinate from a raster. I feel like this function needs to be used somewhere but not sure how it can be implemented with the streams since the apply function is performed on the plot rather than the underlying raster.
def extractRast(rast, lon, lat, size = 100):
# get index of coordinates in array
transform = da.rio.transform()
y_index, x_index = ~ transform * (lon, lat)
# extract required window from array
window = da[int(y_index - size/2):int(y_index + size/2),
int(x_index - size/2):int(x_index + size/2)]
return window
It would be much appreciated if anyone could point me in the right direction.
Many thanks!