Matplotlib set divergent colorbar to zero and change color of nan values

Question

I am plotting sea ice anomalies from the ECMWF's Copernicus data store. The values were from -1 to 1. I multiplied all the values by 100 to get them as a percent. I am having trouble getting the 0% anomaly values to be white with even negative/positive contours. I would like the min and max to be -50 and +50, with values below/above that shown in the same contour. I figured that part out with "extend='both'" for contourf. See example figure. I'd also like to change the color of nan values where there's no data to a blue color, to show areas of open water, but I'm not sure how to go about doing that. The grib files and an example figure can be found here. Below is what I have got so far:

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy
import cartopy.feature as cfeature
import xarray as xr
import datetime as datetime

ds=xr.open_mfdataset('./Data/copernicus-seaice/*.grib', combine='nested', concat_dim='time',
               engine='cfgrib')

# Create variables
lons = ds['longitude'][:]
lats = ds['latitude'][:]
ice = ds['siconc'] # Sea Ice Concentration
time = ds['time'].dt.strftime('%m-%d-%Y') # format time for figure title

# multiply values by 100 to get percentage
ice = ice * 100

for i in range(len(ds.time)):
    fig = plt.figure(figsize=(10,10))
    ax = fig.add_subplot(111, projection=ccrs.NorthPolarStereo(central_longitude=-47))
    ax.set_extent([-180, 180, 60, 90], crs=ccrs.PlateCarree())


    levels = [-50,-45,-40,-35,-30,-25,-20,-15,-10,-5,5,10,15,20,25,30,35,40,45,50]

    plot = ax.contourf(lons, lats, ice.isel(time=i),
                   levels, cmap=plt.cm.seismic_r, extend='both', 
                   transform=ccrs.PlateCarree())

    cb = fig.colorbar(plot, orientation='horizontal', shrink=0.7, pad=0.05)

    cb.set_label('Sea Ice Anomaly (%)', size=12, weight='bold')

    resol = '50m'
    land = cartopy.feature.NaturalEarthFeature('physical', 'land', \
       scale=resol, edgecolor='white', facecolor=cfeature.COLORS['land'])

    # Features
    ax.add_feature(land, facecolor='black')
    ax.add_feature(cfeature.BORDERS, edgecolor='white')

    ax.set_title('Arctic Sea Ice Anomalies: ' + time[i].values, size=14, weight='bold')