Search code examples
python-xarraycartopynetcdf4python-siphon

Sea Surface Temperature NetCDF4 Input z must be 2D, not 3D


I am trying to plot recent sea surface temperature data on a map. I got it working when I downloaded the netCDF4 file, but when I try to access the file from https://www.ncei.noaa.gov/thredds/, I get a TypeError saying, "Input z must be 2D, not 3D". The downloaded file I used originally was from: https://psl.noaa.gov/ Here's what I have:

from netCDF4 import Dataset as netcdf_dataset
from datetime import datetime, timedelta

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from mpl_toolkits.axes_grid1.inset_locator import inset_axes

from siphon.catalog import TDSCatalog
from xarray.backends import NetCDF4DataStore
import xarray as xr

# Catalog only has data from day before so must use timedelta(days=1)
today = datetime.utcnow()
date = datetime(today.year, today.month, today.day, 12) - timedelta(days=1)


base_url = 'https://www.ncei.noaa.gov/thredds/catalog/OisstBase/NetCDF/V2.1/AVHRR/'
cat = TDSCatalog(f'{base_url}{date:%Y%m}/catalog.xml')
ncss = cat.datasets[f'oisst-avhrr-v02r01.{date:%Y%m%d}_preliminary.nc'].subset()

query = ncss.query()
query.time(date)
query.lonlat_box(north=31, south=20, east=283, west=262)
query.accept('netcdf')
query.variables('sst')

data = ncss.get_data(query)
ds = xr.open_dataset(NetCDF4DataStore(data))

sst = ds.variables['sst'][0, :,:]
lats = ds.variables['lat'][:]
lons = ds.variables['lon'][:]

fig = plt.figure(figsize=(15,10))
ax = fig.add_subplot(111, projection=ccrs.PlateCarree())

# TypeError occurs here
sst_contour = ax.contourf(lons, lats, sst, levels=np.arange(0,34,2), cmap='turbo',
                          vmin=0, vmax=38, transform=ccrs.PlateCarree())
isotherm_contour = ax.contour(lons, lats, sst, levels=[27], colors='black', linestyle='--',
                              transform=ccrs.PlateCarree())

This is what my working version looks like using the downloaded file from https://psl.noaa.gov/. Any input on what I'm doing wrong here would be greatly appreciated: enter image description here


Solution

  • I did not consider zlev when pulling the data out of the dataset. Changing

    sst = ds.variables['sst'][0,:,:]
    

    to the following fixed my error

    sst = ds.variables['sst'][0,0,:,:]
    

    enter image description here enter image description here