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)