Adding features to gridded dataset map balloons output file size and slows performance

Software Version Info
bokeh==3.6.2
Cartopy==0.24.1
geoviews==1.14.0
holoviews==1.20.0
ipython==8.32.0
jupyterlab==4.3.5
panel==1.6.0
python==3.13.0
xarray==2025.1.2

Description of expected behavior and the observed behavior

I have been experimenting with GeoViews with netCDF gridded data following the Gridded Datasets I and Gridded Datasets II pages of the user guide. Specifically, I want to flip through timesteps and have the map update quickly as I scroll. I am able to make basic interactive output when I only include the netCDF data. However, adding coastlines, country boundaries and other features rapidly increases the output file size and slows performance.

I am not sure if this is a user issue or if this is a bug. My question is essentially whether there is a better way to add coastlines and other similar features to a map that does not lead to these much larger files and slower interactive behavior.

Expected Behavior

Expected that adding map features (coastlines, state outlines etc.) should only have a marginal impact on file size and performance

Observed Behavior

Adding map features increases the output file size by a factor of 5 (2.2 MB → 11.0 MB) and greatly slows down the interactive map. In other words, following the example below, the output files have the following sizes

"gv_mwe_1.html" : 2.2 MB
"gv_mwe_2.html" : 355 KB
"gv_mwe_3.html" : 11.0 MB

Minimal Working Example

Note that code developed in Jupyter Notebook, but does not appear to be a Jupyter-specific issue

# ============================================================================
# Import Statements
# ============================================================================

import geoviews as gv
import xarray as xr
import cartopy.crs as ccrs
from IPython.display import IFrame

gv.extension('bokeh')
gv.output(widget_location='bottom')

# ============================================================================
# Prep sample data
# ============================================================================

# Load xarray sample air temperature data
ds = xr.tutorial.load_dataset("air_temperature")

# Select the first 24 timesteps of 6-hourly data
ds = ds.isel(time=slice(None,24))

# wrap the xarray dataset as a geoviews Dataset and specify the native CRS
ds_gv = gv.Dataset(ds,crs=ccrs.PlateCarree())

# ============================================================================
# Create layout and map feature
# ============================================================================

proj = ccrs.PlateCarree()

# ----------------------------------------------------------------------------
# QuadMesh to display xarray / netcdf data
# ----------------------------------------------------------------------------
qmesh = ds_gv.to(gv.QuadMesh,kdims=['lon','lat'],vdims=['air','lat','lon'])

qmesh_layout = qmesh.opts(backend='bokeh',colorbar=True,projection=proj,)

# ----------------------------------------------------------------------------
# gv.features for coastlines, state boundaries etc.
# ----------------------------------------------------------------------------

map_features =  (
    gv.feature.coastline().opts(line_color='black') * \
    gv.feature.borders() * \
    gv.feature.states() * \
    gv.feature.lakes().opts(line_color='black',fill_color='none')
).opts(global_extent=True,projection=proj)

# ============================================================================
# Write Output Files
# ============================================================================

# ----------------------------------------------------------------------------
# File 1: Qmesh Only
# ----------------------------------------------------------------------------
f1 = "./gv_mwe_1.html"
output1 = qmesh_layout
gv.save(output1,f1)

# ----------------------------------------------------------------------------
# File 2: Map Features Only
# ----------------------------------------------------------------------------
f2 = "./gv_mwe_2.html"
output2 = map_features
gv.save(output2,f2)

# ----------------------------------------------------------------------------
# File 3: Qmesh and Map Features
# ----------------------------------------------------------------------------
f3 = "./gv_mwe_3.html"
output3 = (qmesh_layout * map_features).opts(global_extent=False)
gv.save(output3,f3)

# ============================================================================
# Display Outputs (Jupyter)
# ============================================================================

IFrame(f1,width=400,height=400) # <--- Works fine
IFrame(f2,width=400,height=400) # <--- Works fine
IFrame(f3,width=400,height=400) # <--- Extremely slow and unresponsive

You can crop the features to the bounds you care about How to show cropped geometries / features upon zoom.

Alternatively, you can use tiles instead.

Thanks @ahuang11, your solution of cropping the features improved the performance. For my case I specified a fixed scale and didn’t worry about changing the scale upon zooming. The final file (including features and gridded data) is still significantly larger than the file with only gridded data (4.4 MB vs 2.2 MB), but that is fine for my applications.

For future reference, here’s an updated example of how I accomplished this based on the link @ahuang11 provided. The only change is that now I create the map_features variable from cropped gpd GeoDataFrames:

# ============================================================================
# Import Statements
# ============================================================================

import geoviews as gv
import xarray as xr
import cartopy.crs as ccrs
from cartopy import feature as cf
import geopandas as gpd

gv.extension('bokeh')
gv.output(widget_location='bottom')

# ============================================================================
# Prep sample data
# ============================================================================

# Load xarray sample air temperature data
ds = xr.tutorial.load_dataset("air_temperature")

# Select the first 24 timesteps of 6-hourly data
ds = ds.isel(time=slice(None,24))

# wrap the xarray dataset as a geoviews Dataset and specify the native CRS
ds_gv = gv.Dataset(ds,crs=ccrs.PlateCarree())

# ============================================================================
# Create layout and map feature
# ============================================================================

proj = ccrs.PlateCarree()

# ----------------------------------------------------------------------------
# QuadMesh to display xarray / netcdf data
# ----------------------------------------------------------------------------
qmesh = ds_gv.to(gv.QuadMesh,kdims=['lon','lat'],vdims=['air','lat','lon'])

qmesh_layout = qmesh.opts(backend='bokeh',colorbar=True,projection=proj,)

# ----------------------------------------------------------------------------
# cartopy.features for coastlines, state boundaries
# ----------------------------------------------------------------------------

# Get data bounding box based on the dataset (convert longitudes
# from [0,360] to [-180,180] by subtracting 360
bbox = [float(x) for x in \
        [ds.lon.min()-360,ds.lat.min(),ds.lon.max()-360,ds.lat.max()]
       ]

# Specify the shapefile scale, one of 110m, 50m, 10m
scale = "110m"

states = gpd.clip(
    gpd.GeoDataFrame(
        geometry=gpd.GeoSeries(cf.STATES.with_scale(scale).geometries())
    ),
    bbox,
)

coastline = gpd.clip(
    gpd.GeoDataFrame(
        geometry=gpd.GeoSeries(cf.COASTLINE.with_scale(scale).geometries())
    ),
    bbox,
)

map_features = gv.Path(coastline).opts(color='black') * gv.Path(states).opts(color='black')


# ============================================================================
# Write Output Files
# ============================================================================

# ----------------------------------------------------------------------------
# File 1: Qmesh Only
# ----------------------------------------------------------------------------
f1 = "./gv_mwe_1.html"
output1 = qmesh_layout
gv.save(output1,f1)

# ----------------------------------------------------------------------------
# File 2: Map Features Only
# ----------------------------------------------------------------------------
f2 = "./gv_mwe_2.html"
output2 = map_features
gv.save(output2,f2)

# ----------------------------------------------------------------------------
# File 3: Qmesh and Map Features
# ----------------------------------------------------------------------------
f3 = "./gv_mwe_3.html"
output3 = (qmesh_layout * map_features).opts(global_extent=False)
gv.save(output3,f3)
1 Like

Note the related issue on the GeoViews github page.