Add lat lon to Seviri plots

[SamothKoc]

Hello, I'm new to Satpy and I'm facing a bit of trouble. I've tried various approaches and couldn't find a solution through Google either. I'm able to work with natural-color images without any issues. I want to integrate the image into Cartopy to add latitude and longitude gridlines. My next step is to incorporate model data, which has worked well for other satellites.

Even with natural-color images, everything is fine. However, when I use imshow() on the Seviri dust product, the image suddenly turns blue. When I use the show() method itself and work with the temporarily generated image, it appears correct. Is there a trick to make imshow() display the normal RGB dust color scheme? It seems that imshow() is altering the color scheme between 0 and 1. .show(): imshow:

####DUST
composite_id = ['dust']
scn.load(composite_id, upper_right_corner="NE")
scn_cropped = scn.resample(area)
#scn_cropped.show('dust')

scn_cropped['dust'].attrs['area']
area = scn_cropped['dust'].attrs['area']

scn_resample_nc = scn.resample(area)
#scn_resample_nc.show('dust')  -> Correct image without grid and no automatic option

#Cartopy
image = np.asarray(scn_resample_nc["dust"]).transpose(1,2,0)
image.shape

image = np.interp(image, (np.percentile(image,1), np.percentile(image,99)), (0, 1))
crs = scn_resample_nc["dust"].attrs["area"].to_cartopy_crs()

# Initiatie a subplot and axes with the CRS information previously defined
fig = plt.figure(figsize=(16, 10))
ax = fig.add_subplot(1, 1, 1, projection=crs)

# Add coastline features to the plot
ax.coastlines(resolution="10m", color="white")

# Define a grid to be added to the plot
gl = ax.gridlines(draw_labels=True, linestyle='--', xlocs=range(-20,20,5), ylocs=range(0,50,5))
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}

# Bounderies
countries = cfeature.NaturalEarthFeature(
    category='cultural',
    name='admin_0_countries',
    scale='50m',
    facecolor='none',
    edgecolor='white'
)
ax.add_feature(countries)

# In the end, we can plot our image data...
ax.imshow(image, transform=crs, extent=crs.bounds, origin="upper")

# Define a title for the plot
plt.title("RGB Natural color composite recorded by MSG at " + scn_resample_nc['dust'].attrs["start_time"].strftime("%Y-%m-%d %H:%M"), fontsize=20, pad=20.0)

plt.savefig('plot.png', dpi=300, bbox_inches='tight', transparent=True)

I hope someone can assist me with this.

[djhoese]

You'll need to "enhance" your RGB image so that it gets each band stretched to the expected limits. Satpy's .show does this for you. Satpy also does this as part of the writing step in .save_datasets. The simple workaround is to change the code like this:

from satpy.writers import get_enhanced_image
image = get_enhanced_image(scn_resample_nc["dust"]).data.data.compute().transpose(1, 2, 0)

instead of what you're doing now which is this:

image = np.asarray(scn_resample_nc["dust"]).transpose(1,2,0)

Note that get_enhanced_image is returning a special XRImage object. The first .data is the underlying enhanced DataArray object. The second .data is the dask array inside the DataArray. .compute() on the dask array returns a numpy array and then you can do the transpose as you did before.

[djhoese]

Now for my question, what does this line do:

image = np.interp(image, (np.percentile(image,1), np.percentile(image,99)), (0, 1))

[SamothKoc]

Well... imshow scales the color range between 0 and 1, so you can adjust the scaling to customize the color scale. That's why I thought... In summary, this line performs a contrast adjustment where the pixel values of the image are scaled to increase contrast and range between 0 and 1 to highlight and customize specific features in the image.... that was my first idea :D but it doesn't workd :D...

[djhoese]

Ah, if I'm understanding you correctly that shouldn't be needed once get_enhanced_image is used. That normalization to a 0 to 1 range is done by the "enhancement".

[SamothKoc]

Yes :D that is it... A further question... is there a way to get the pixel coordinate to combine model data with the current Seviri plt?

Greetings

[djhoese]

Just to make sure we're talking about the same thing, what do you mean by "pixel coordinate"?

Normally I would recommend resampling the SEVIRI data and the model data to the same projection/grid (an AreaDefinition in satpy/pyresample terms).

[SamothKoc]

Pyresample was the key :) thanks.