MSG Dust RGB adding coastilnes and grid to the image

[sofmaqueo]

Hi
I am trying to plot the MSG' Dust RGB composite

So far I am using Scene and I manage to get an image of the composite but I want to add coastilines and a grid as overlays, but I have not managed to do that.

First, how it is working:

file_name = glob.glob('/path/MSG-SEVIRI/data/Feb/18/MSG4-SEVI-MSG15-0100-NA-20220218120010.124000000Z-20220218121300-4852904-3.nat')
scn = Scene(reader="seviri_l1b_native", 
             filenames=file_name)
area = scn_cropped['dust'].attrs['area']
scn_resample_nc = scn.resample(area)
scn_resample_nc.show('dust')

Then, I tried to convert the scene file to a numpy array to plot with cartopy:

image = np.asarray(scn_resample_nc["dust"]).transpose(1,2,0)
image = np.nan_to_num(image)
crs = scn_resample_nc["dust"].attrs["area"].to_cartopy_crs()

And for the plot:

fig = plt.figure(figsize=(16, 10))
ax = fig.add_subplot(1, 1, 1, projection=crs)

ax.coastlines(resolution="10m", color="white")

gl = ax.gridlines(draw_labels=True, linestyle='--', xlocs=range(-30,40,1), ylocs=range(0,40,1))
gl.top_labels=False
gl.right_labels=False
gl.xformatter=LONGITUDE_FORMATTER
gl.yformatter=LATITUDE_FORMATTER
gl.xlabel_style={'size':14}
gl.ylabel_style={'size':14}

ax.imshow(image, transform=crs, extent=crs.bounds, origin="upper")

The coastilines, and gridlines are fine but I want the colors of the first graph to be on the second, somehow I haven't managed to do that.
Can someone help me? I think it is something about how the colors are being adjudicated to the numbers in the imshow function but I do not know how to change that...

Any help would be highly appreciated!

Cheers,
Sofía

[djhoese]

The .show() method of the Scene is doing one extra step internally that you're missing in your own work. It is "enhancing" the image. You can perform this step manually by doing:

from satpy.writers import get_enhanced_image

data_arr = get_enhanced_image(scn["dust"]).data

Note the .data at the end. Then data_arr will be a xarray DataArray. Additionally, you could simplify your code by doing:

data_arr.plot.imshow(ax=ax, rgb="bands", transform=crs, extent=crs.bounds, origin="upper")

I don't remember if the origin= is needed with this. This is a shortcut that uses xarray to do semi-standard things for you (like your transpose call which could be removed if you used this).

[sofmaqueo]

Thank you so much! It works like a charm :)