Interactive Satellite Map using Python

Question

I'm trying to overlay a Lambert Conformal Conical satellite image onto a Holoviews interactive map. I can map the satellite image just fine, but I can't figure out how to translate this map onto a Holoviews map properly. Below is reproducible code where I grab data using the Unidata Siphon library.

Import Packages

from datetime import datetime
import matplotlib.pyplot as plt
from netCDF4 import Dataset
from siphon.catalog import TDSCatalog
import holoviews as hv
import geoviews as gv
import geoviews.feature as gf
from cartopy import crs
from cartopy import feature as cf

hv.extension('bokeh')

Grab data and create figure

date=datetime.utcnow()
idx=-2
regionstr = 'CONUS'
channelnum = 13
datestr = str(date.year) + "%02d"%date.month + "%02d"%date.day
channelstr = 'Channel' + "%02d"%channelnum

cat = TDSCatalog('http://thredds-test.unidata.ucar.edu/thredds/catalog/satellite/goes16/GOES16/' + regionstr +
    '/' + channelstr + '/' + datestr + '/catalog.xml')
ds = cat.datasets[idx].remote_access(service='OPENDAP')
x = ds.variables['x'][:]
y = ds.variables['y'][:]
z = ds.variables['Sectorized_CMI'][:]

proj_var = ds.variables[ds.variables['Sectorized_CMI'].grid_mapping]

# Create a Globe specifying a spherical earth with the correct radius
globe = ccrs.Globe(ellipse='sphere', semimajor_axis=proj_var.semi_major,
                   semiminor_axis=proj_var.semi_minor)

proj = ccrs.LambertConformal(central_longitude=proj_var.longitude_of_central_meridian,
                             central_latitude=proj_var.latitude_of_projection_origin,
                             standard_parallels=[proj_var.standard_parallel],
                             globe=globe)


fig = plt.figure(figsize=(14, 10))
ax = fig.add_subplot(1, 1, 1, projection=proj)
ax.coastlines(resolution='50m', color='lightblue')
ax.add_feature(cf.STATES, linestyle=':', edgecolor='lightblue')
ax.add_feature(cf.BORDERS, linewidth=1, edgecolor='lightblue')

for im in ax.images:
    im.remove()
im = ax.imshow(z, extent=(x.min(), x.max(), y.min(), y.max()), origin='upper', cmap='jet')
plt.colorbar(im)

Now plot an interactive image using Holoviews (which uses Bokeh backend)

goes = hv.Dataset((x, y, z),['Lat', 'Lon'], 'ABI13')
%opts Image (cmap='jet') [width=1000 height=800 xaxis='bottom' yaxis='left' colorbar=True toolbar='above' projection=proj] 
goes.to.image()* gf.coastline().options(projection=crs.LambertConformal(central_longitude=proj_var.longitude_of_central_meridian,central_latitude=proj_var.latitude_of_projection_origin,standard_parallels=[proj_var.standard_parallel],globe=globe))

I must not be translating it properly, though I've found documentation on Holoviews with regards to Lambert Conformal Conical projection to be sparse. I'm open to using any other interactive map package. My main desire is to be able to plot relatively quickly, get state/country lines on the image properly, and be able to zoom in. I've tried folium but also fell into projection issues there.

Answer 1

So I think the main thing to understand, which is explained here: is how projections are declared. The elements (e.g. Image, Points etc.) in GeoViews have a parameter called crs which declares the coordinate system the data is in, while the projection plot option declares what to project the data to for display.

In your case I think you want to display the image in the same coordinate system it is already in (Lambert Conformal), so technically you don't have to declare the coordinate system (crs) on the element at all and can just use a hv.Image (which is entirely unaware of projections).

As far as I can tell your code should already work as expected if you are using GeoViews 1.5, but here is what I would do:

# Apply mask
masked = np.ma.filled(z, np.NaN)

# Declare image from data
goes = hv.Image((x, y, masked),['Lat', 'Lon'], 'ABI13')

# Declare some options
options = dict(width=1000, height=800, yaxis='left', colorbar=True,
               toolbar='above', cmap='jet', projection=proj)

# Display plot
gf.ocean * gf.land * goes.options(**options) * gf.coastline.options(show_bounds=False)

Note how we declare the projection on the Image but not crs. If you do want to display the data in a different projection it is defined in, you do have to declare the crs and use a gv.Image. In this case I'd recommend using project_image with the fast option enabled (which might introduce some artifacts but is much faster):

# Apply mask
masked = np.ma.filled(z, np.NaN)

# Declare the gv.Image with the crs
goes = gv.Image((x, y, masked), ['Lat', 'Lon'], 'ABI13', crs=proj)

# Now regrid the data and apply the reprojection
projected_goes = gv.operation.project_image(goes, fast=False, projection=ccrs.GOOGLE_MERCATOR)

# Declare some options
options = dict(width=1000, height=800, yaxis='left', colorbar=True,
               toolbar='above', cmap='jet')

# Display plot
projected_goes.options(**options) * gv.tile_sources.ESRI.options(show_bounds=False)

Another final tip, when you plot with bokeh all the data you are plotting will be sent to the browser, so when dealing with images any larger than you are already using I'd recommend using the holoviews' regrid operation which uses datashader to dynamically resize the image as you are zooming. To use that simply apply the operation to the image like so:

from holoviews.operation.datashader import regrid
regridded = regrid(goes)