Issues while plotting a ternary map from Xarray with Cartopy

[rlourenco]

Hi everyone. I’m trying to create a ternary composition (I have RGB channels as Data Arrays, in a Xarray Dataset) and almost getting things well. Currently, code is like this:

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 calling the plot with lambert_azimuthal_rgb_plot(xr_from_df_arctic, title='RGB Map', latitude_cutoff=50, img_multiplier=1), the result is like this:

Things are almost okay, but the results are padded in the longitude (roughly by 10 degrees). However, when I plot the original channels ( ex.: xarray_dataset['red'] ) all is on the correct aspect.

Do you have any suggestions on what may be going on? FYI, I have cross-posted this on Cartopy gitter and at the Xarray discussion (RGB image gets padded along dimension when plotting · pydata/xarray · Discussion #9240 · GitHub ); however, @TomNicholas pointed out that it was slightly off-topic, which is why I am posting it here. Also other suggestions on plotting ternary map compositions are welcome.

[rlourenco]

Just including here an example if you want to try to replicate the problem (in the notebook I also plot the individual channels, which show up correctly):

Jupyter Notebook Reproducible example: Xarray-Cartopy RGB channels plot reproducible example · GitHub

NetCDF4 file for the notebook: Dropbox

[Michael_Sumner]

this file is not georeferenced in the normal raster way, the rectlinear coords for lon and lat have gaps and non-regular sequence

It’s not obvious if it’s really meant to be represented that way (there isn’t much metadata in there). I think you need to read it via a regridding process (like GDALwarp), or plot it as a mesh - it 's not obvious how that should be represented with these data, and how that works in Python is not my strength.

This is where the error comes from, we can’t assume this is a regular grid with predictable spacing between the min and max coordinates:

  # 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()]

This is how it looks when plotted as the (expanded) set of points in the LAEA projection:

I can show some examples of using GDALwarp (via cli or osgeo.gdal in python), but that kind of regridding process is where you might need some expertise on how this data should be being handled. (I would experiment with heuristics for the output regular grid spec in LAEA and resampling until it “looks good”, but with RGB and actual specific longlat points for these data it’s not obvious that interpolation would be a good idea. It mght be best to resolve to a regular longlat grid first, and then reproject that). It looks to me like a point dataset, perhaps one that was a full grid but has only kept certain meridians … (or something).

this is what the long and lat coordinates look like individually, and with their sequential differences showing why this can’t be treated like a regular grid. The big jump in longitude by 10 is just the most obvious step that causes a problem under our “regular grid” assumption:

[rlourenco]

Hi @Michael_Sumner ! Thanks for the thorough analysis.

In the time between my post and yours, I was able to solve the problem by both using scikit-image img_as_ubyte and manually presetting the grid definition (your point on the irregular grid was a major issue - in my case, my data preprocessing masks out data given some criteria [ex.: managed vs native vegetated areas] and this is likely the culprit).

I wasn’t able to use the presetting strategy with pcolormesh / imshow , though (probably my deficiencies in matplotlib are more evident here).

Anyway, this is how it became: