RGB image gets padded along dimension when plotting

[ricardobarroslourenco]

I have developed the following function to plot a RGB composite, from data arrays in a dataset:

def lambert_azimuthal_rgb_plot(xarray_dataset, title='Empty Title', latitude_cutoff=50, img_multiplier=1):
    # Extract the red, green, and blue data arrays
    red = xarray_dataset['red'].values
    green = xarray_dataset['green'].values
    blue = xarray_dataset['blue'].values

    # Combine red, green, and blue arrays into an RGB image
    rgb = np.dstack((red, green, blue))

    # Create a mask to filter out NaN values
    mask = np.isnan(rgb).any(axis=2)
    
    # Set NaN values to transparent (0, 0, 0, 0) in RGBA
    rgba = np.zeros((rgb.shape[0], rgb.shape[1], 4))
    rgba[~mask, :3] = rgb[~mask]
    rgba[~mask, 3] = 1  # Set alpha channel to 1 for non-NaN values

    # Define the coordinate reference system (CRS) for the data
    data_crs = ccrs.PlateCarree()

    # Plot the RGB image using cartopy with Lambert Azimuthal Equal Area projection
    fig, ax = plt.subplots(figsize=(10, 10),
                           subplot_kw={'projection': ccrs.LambertAzimuthalEqualArea(central_latitude=90,
                                                                                    central_longitude=0)})
    ax.set_title(title)

    # Add 50m resolution coastlines and country borders
    coastlines_50m = cfeature.NaturalEarthFeature('physical', 'coastline', '50m', edgecolor='black', facecolor='none')
    borders_50m = cfeature.NaturalEarthFeature('cultural', 'admin_0_boundary_lines_land', '50m', edgecolor='black', facecolor='none')
    ax.add_feature(coastlines_50m, linewidth=1)
    ax.add_feature(borders_50m, linestyle=':', linewidth=1)

    # Calculate the extent of the data from xarray metadata
    extent = [xarray_dataset['x'].min().item(), 
              xarray_dataset['x'].max().item(), 
              xarray_dataset['y'].min().item(), 
              xarray_dataset['y'].max().item()]

    # Display the RGBA image
    ax.imshow(rgba,
              origin='lower',  # Using 'lower' as appropriate for this data
              extent=extent,
              transform=data_crs,
              interpolation='none'
             )

    # Set the extent of the map to include the area of interest
    ax.set_extent([-180, 180, latitude_cutoff - 3, 90], crs=data_crs)

    # Adding gridlines
    gl = ax.gridlines(draw_labels=True, color='gray', alpha=0.5, linestyle='--')
    gl.xlocator = plt.FixedLocator(np.arange(-180, 181, 10))
    gl.ylocator = plt.FixedLocator(np.arange(50, 91, 10))
    gl.xlabel_style = {'size': 10, 'color': 'gray'}
    gl.ylabel_style = {'size': 10, 'color': 'gray'}

    plt.show()

When doing a basic call, lambert_azimuthal_rgb_plot(xr_from_df_arctic, title='RGB Map', latitude_cutoff=50, img_multiplier=1) , it generates the following image:

It is almost correct, but data is misaligned along the x axis (longitude - roughly by 10 degrees). Any hint on what may be going on? When I do plot the individual rgb channels, they are all properly aligned.

---

Edit 2024-Jul-14:

Jupyter Notebook Reproducible example: https://gist.github.com/ricardobarroslourenco/bfb3e4bc52118f347ccd16633a00f0b0

NetCDF4 file for the notebook: https://www.dropbox.com/scl/fi/fvxz028eeb11064oj71t9/reproducible_rgb.nc4?rlkey=uclhs2hdxv9m3c7xj7nkiu4so&st=yvu6z6k4&dl=0

[TomNicholas]

I'm not sure that this question has anything to do with xarray itself. The first step of your example extracts the numpy arrays from the xarray objects, and from that point everything else is either matplotlib or numpy. You might have more luck asking on the matplotlib issue tracker.

[ricardobarroslourenco]

Hi @TomNicholas! Thanks for pointing this out.

I am trying to do a basic ternary composite, so I posted it here in case someone else has had the same issue. As I wasn't able to do it natively with xr.Dataset.plot I exported, but I would be happy to do this within Xarray.

[mathause]

It's also difficult to help without an easy way to reproduce the example.

Is extent given for the center coordinates but should be the edges? Does it plot corectly with pcolormesh?

[ricardobarroslourenco]

Hi @mathause ! I have included an edit with the reproducible example (Jupyter Notebook file and netcdf) :)

Regarding the extent, I believe that when setting with origin lower , it already considers the edges for anchoring the image limits ( https://matplotlib.org/stable/users/explain/artists/imshow_extent.html#default-extent ). Regarding your pcolormesh question, it works fine (refer to the individual channel plots).

[ricardobarroslourenco]

@mathause I have discovered the issue. I was plotting with imshow instead of pcolormesh in this snippet. As I have masked out data, the x and y coordinates were dropped. However, it seems that xarray.DataArray.plot uses pcolormesh by default, so it is robust to these changes, while imshow is not. I will be closing it here, but putting a note on my end on using xarray with an irregular coordinate grid on matplotlib.