Problem on visualizing Himawari-8 NetCDF data using cartopy

I'm trying to plot Himawari-8 AHI Data (downloaded from THREDDs) using Python and following this simple notebook (Go to Himawari-8 Section). Everything went alright when first plotting it through imshow.

from netCDF4 import Dataset
import matplotlib.pyplot as plt
import numpy as np
import cartopy.crs as ccrs

# Data paths
h8b1_path = '0800_20160814080000-P1S-ABOM_OBS_B01-PRJ_GEOS141_2000-HIMAWARI8-AHI.nc'
h8b2_path = '0800_20160814080000-P1S-ABOM_OBS_B02-PRJ_GEOS141_2000-HIMAWARI8-AHI.nc'
h8b3_path = '0800_20160814080000-P1S-ABOM_OBS_B03-PRJ_GEOS141_2000-HIMAWARI8-AHI.nc'

# Open data
h8b1 = Dataset(h8b1_path)
h8b2 = Dataset(h8b2_path)
h8b3 = Dataset(h8b3_path)

for item in h8b1.dimensions:
    print("dimensions: {}, name: {}\n".format(h8b1.dimensions[item].name, 
                                              h8b1.dimensions[item].size))

vars = h8b1.variables.keys()
for item in vars:
    print('Variable:\t{}'.format(item))
    print('Dimensions:\t{}'.format(h8b1[item].dimensions))
    print('Shape:\t{}\n'.format(h8b1[item].shape))

b1 = h8b1.variables['channel_0001_scaled_radiance'][0,:,:]
b2 = h8b2.variables['channel_0002_scaled_radiance'][0,:,:]
b3 = h8b3.variables['channel_0003_scaled_radiance'][0,:,:]

x = h8b1.variables['x'][:]
y = h8b1.variables['y'][:]

b1 = h8b1.variables['channel_0001_scaled_radiance'][0,:,:]
b2 = h8b2.variables['channel_0002_scaled_radiance'][0,:,:]
b3 = h8b3.variables['channel_0003_scaled_radiance'][0,:,:]

x = h8b1.variables['x'][:]
y = h8b1.variables['y'][:]

fig = plt.figure(figsize=(12, 12))
spec = fig.add_gridspec(2, 4, hspace=0, wspace=0.1)

# Doing something here as to why I used gridspec.

ax0 = fig.add_subplot(spec[0, :])
ax0.imshow(rgb)

plt.show()

However, when using/applying a CRS from Cartopy and replot the 'Himawari subset', it just shows a blank figure.

plt.figure(figsize=(12,12))

# Or perhaps the projection scheme being used???

# Define the projection information of the image
img_proj = ccrs.Geostationary(central_longitude=h8b1['geostationary'].longitude_of_projection_origin, 
                                  satellite_height=h8b1['geostationary'].satellite_height)

# The extents of the image we are plotting
img_extent = (X[0], X[-1], Y[-1], Y[0])

# Setup the axes projection
ax = plt.axes(projection=ccrs.Miller(central_longitude=float(h8b1['geostationary'].longitude_of_projection_origin)))

# Extent of the axes in the plot
ax.set_extent([117, 125, 13, 20]) # Area of Luzon Landmass

# Plot image
ax.imshow(rgb, transform=img_proj, extent=img_extent, origin='upper')

Is there a way to fix this out? Attached here is the sample data as well as in the code attached.

Solution

You don't show how you create the rgb variable that you're trying to plot, which makes it difficult to know what is happening. Keep in mind that matplotlib is flexible regarding the range of values when using an RGB(A) definition, but this flexibility might bite you as well. Integers are assumed to be in the 0-255 range, and floats in 0-1. So if your RGB (array?) is floats from 0-255, it would perhaps all clip to white.

Something that's certainly going wrong is a mismatch in units. The data projection assumes meters, but the NetCDF data provides x/y in kilometers. This would cause the data to be plotted as a tiny image near the origin of the projection, perhaps also outside the map extent you set.

The example below shows how you can resolve this. I've changed a few things which are subjective to some extent.

Using xarray to make reading the data slightly easier, but NetCDF4 should work fine as well if you want to avoid the dependency.
Using pcolorfast over imshow for plotting. I tend to avoid on-the-fly warming of images with imshow (different data vs map projection) with Cartopy, especially with projections like Geostationary that contain points that aren't even on earth. This is also why I read only a subset (crudely with .isel), which speeds things up a little but also avoids those "singularities" when reprojecting.

import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from pathlib import Path
import numpy as np

# I/O
data_dir = Path(r"C:\Temp")

with xr.open_mfdataset((
    data_dir / '0800_20160814080000-P1S-ABOM_OBS_B01-PRJ_GEOS141_2000-HIMAWARI8-AHI.nc',
    data_dir / '0800_20160814080000-P1S-ABOM_OBS_B02-PRJ_GEOS141_2000-HIMAWARI8-AHI.nc',
    data_dir / '0800_20160814080000-P1S-ABOM_OBS_B03-PRJ_GEOS141_2000-HIMAWARI8-AHI.nc',
)) as ds:

    ds = ds.isel(x=slice(1000, 4000), y=slice(1000, 3000), time=0)

    central_longitude = ds['geostationary'].longitude_of_projection_origin
    satellite_height = ds['geostationary'].satellite_height

    # convert km to m
    mapx = ds['x'].to_numpy() * 1000.
    mapy = ds['y'].to_numpy() * 1000.

    rgb = np.stack((
        ds['channel_0001_scaled_radiance'].to_numpy(),
        ds['channel_0002_scaled_radiance'].to_numpy(),
        ds['channel_0003_scaled_radiance'].to_numpy(),
    ), axis=-1)


# stretch RGB values a little
rgb_stretched = np.clip((rgb/0.7)**0.7, 0, 1)

# define projections
data_proj = ccrs.Geostationary(
    central_longitude=central_longitude, 
    satellite_height=satellite_height,
)
map_proj =  ccrs.Miller(central_longitude=central_longitude)

# plotting
fig, ax = plt.subplots(
    figsize=(6, 6), facecolor="w", dpi=100, layout="compressed",
    subplot_kw=dict(projection=map_proj),
)

pcm = ax.pcolorfast(mapx, mapy, rgb_stretched, transform=data_proj)

ax.coastlines(color="w")
ax.set_extent([115, 127, 11, 22])

I've scaled the raw RGB a little just for visualization purposes, and also set the extent slightly wider than you to show a bit more context.

For the relatively tiny subset that you have I probably wouldn't change the map projection to something different than the data projection. That just introduces an extra resampling without much added benefit (to me), it's perfectly valid to do so of course.