Some composite are produced with the wrong colors

[guidocioni]

Describe the bug
I'm using SEVIRI L1.5 Rapid scan data.
When creating some of the composites the resulting data has wrong RGB channels data, resulting in the wrong image being produced.

To Reproduce

channel='airmass'
scn = Scene(filenames = {reader:[file]})
scn.load([channel])

area = create_area_def('my_area',
                        {'proj': 'latlong', 'lon_0': 0},
                        area_extent=[-2, 36, 22, 55],
                        resolution=0.005,
                        units='degrees',
                        description='test')

scn = scn.resample(area)

### FIRST METHOD (using matplotlib and cartopy) 

image = np.asarray(scn[channel]).transpose(1,2,0)
image = np.interp(image, (np.nanpercentile(image,1), np.nanpercentile(image,99)), (0, 1))
crs = scn[channel].attrs["area"].to_cartopy_crs()


fig = plt.figure(figsize=(20,20))
ax = fig.add_subplot(1, 1, 1, projection=crs)
ax.coastlines(resolution="50m", color="white")

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

plt.savefig('test_2.png', dpi=120, bbox_inches='tight')

### SECOND METHOD

overlay = {
    'coast_dir': './gshhg-shp-2.3.7/',
    'color': (255, ),
    'resolution': 'i',
    'width': 5
}
scn.save_dataset(channel, filename='test.png', overlay=overlay)

Expected behavior
The images produced by the two methods should coincide.

Actual results
Image obtained with matplotlib

Image obtained with second method

Some of the other composites are fine, e.g. for natural color
Image obtained with matplotlib

Image obtained with second method

The single channels image are also fine

Environment Info:

• OS: MacOS
• Satpy Version: 0.35.0
• PyResample Version: 1.23.0

[djhoese]

The RGB data needs to be "enhanced" before being plotted with matplotlib. See the mention of get_enhanced_image in the tutorial:

https://github.com/pytroll/tutorial-satpy-half-day/blob/644a0b9419c628c00f66a94064d7cb70f11b10d6/notebooks/05_composites.ipynb

[guidocioni]

yes, you're right! I didn't see that tutorial.
Now I get exactly the same image with the two methods.

Any clue why cartopy does not draw some of the borders? it only happens at the fine resolution.

BTW here is the (much cleaner) code for reference which plots both normal channel and composites.

channel='airmass'
scn = Scene(filenames = {reader:[file]})
scn.load([channel])

area = create_area_def('my_area',
                        {'proj': 'latlong', 'lon_0': 0},
                        area_extent=[-2, 36, 22, 55],
                        resolution=0.005,
                        units='degrees',
                        description='test')

scn = scn.resample(area)
image = get_enhanced_image(scn[channel]).data
crs = scn[channel].attrs["area"].to_cartopy_crs()

fig = plt.figure(figsize=(20,20))
ax = fig.add_subplot(1, 1, 1, projection=crs)
ax.coastlines(resolution="10m", color="white")

if len(scn[channel].shape) > 2:
    image.plot.imshow(vmin=0, vmax=1, add_colorbar=False, rgb='bands', ax=ax)
else:
    image[0].plot.imshow(vmin=0, vmax=1, cmap='Greys_r', add_colorbar=False, ax=ax)

[djhoese]

Probably not a problem, but try renaming scn = scn.resample(area) to a new name like new_scn = scn.resample(area). Then make sure you update the code below it to use that new variable. This is more of a best practice thing than fixing your specific problem.

Very odd with the coastlines. I'm not sure why it would not be able to do some of those. It is important to note that pycoast (the pytroll package drawing the coastlines in your code) and cartopy use two different sets of shapefiles. Pycoast uses GSHHG and cartopy uses NaturalEarth. However, the polygons missing are pretty easy ones to draw (not over the dateline, not going outside the extents of the image, etc). 🤔 Or maybe they are going outside of the image extents. You could try looking at forcing the x/y extents to be larger just to see if it has an effect on where the coastlines show up.

[djhoese]

I'm going to close this issue as I think the main satpy stuff is figured out, but feel free to keep commenting here and we can continue the discussion/debugging.

[guidocioni]

Hey @djhoese , thank you a lot for the support, the coastlines are definitely not satpy issue; as a matter of fact I've seen this problem popping up randomly in the past and I'm not sure where it's coming from.

Anyway I'm exploiting your kindness to do some additional debugging :)

I'm trying to write a processing script to operationally plot SEVIRI rapid scan native images. I'm using the EUMETSAT APIs to download the latest scans and then spawn some parallel processes with matplotlib to do all the plotting. I'm quite close but I still don't get how can I get the same quality output as when using the Scene.show() method.

What I want to achieve is to be able to plot

• with cartopy (so that I can manipulate better the overlay afterwards. Also the rapid scan image is upside down)
• a local domain that I choose with a lat-lon box
• data at the highest resolution available
• without resampling (if possible!): as it is operational the least amount of cpu time I use the better!

Here is the code that I'm using (only pasting the minimum required, although the rest handling the file downloading and parallelization is quite interesting, I could share it in the future).

# channel is either HRV or a composite name
scn = Scene(filenames={"seviri_l1b_native": [fname]})
scn.load([channel])
if channel not in scn: # means we're plotting composite so we have to resample
    area = create_area_def('my_area',
                           {'proj': 'latlong', 'lon_0': 0},
                           area_extent=area,
                           resolution=0.01, # chosen arbitrarily... don't know what to put
                           units='degrees',
                           description='test')
    scn = scn.resample(area)
else: # normal channel, no resample needed, just crop
    scn = scn.crop(ll_bbox=tuple(area))
#
image = get_enhanced_image(scn[channel]).data
crs = scn[channel].attrs["area"].to_cartopy_crs()

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

if len(scn[channel].shape) > 2:
    image.plot.imshow(vmin=0, vmax=1, add_colorbar=False,
                      rgb='bands', ax=ax)
else:
    image[0].plot.imshow(vmin=0, vmax=1, cmap='Greys_r',
                         add_colorbar=False, ax=ax)

I cannot seem to get the same quality as in the .show() method. Here is a comparison of the same scene: the image with the coastlines is obviosuly the one obtained with cartopy.

Do you think the problem is related to the way imshow does the interpolation? How can I be sure that I'm plotting the data version with the highes possible resolution every time, regardless of whether it is a composite or a normal channel?

[djhoese]

a local domain that I choose with a lat-lon box

You are already using the Scene.crop method so that should be fine. Instead of doing regular resampling you can do Scene.resample(resampler='native') for geostationary data and it should allow the composites to be generated without having to do the more expensive real resampling.

data at the highest resolution available

As long as you provide all the resolutions as input then Satpy should default to loading the highest resolution it can find.

without resampling (if possible!)

I'm not an expert on the SEVIRI readers, but maybe @pnuu or @ameraner can point to some of the options with those readers. I think there is one to force the data to be flipped inside the reader so that it doesn't get produced upside down.

[guidocioni]

a local domain that I choose with a lat-lon box

You are already using the Scene.crop method so that should be fine. Instead of doing regular resampling you can do Scene.resample(resampler='native') for geostationary data and it should allow the composites to be generated without having to do the more expensive real resampling.

I get a ValueError: Expand factor must be a whole number if I try to use that resampler.

data at the highest resolution available

As long as you provide all the resolutions as input then Satpy should default to loading the highest resolution it can find.

without resampling (if possible!)

I'm not an expert on the SEVIRI readers, but maybe @pnuu or @ameraner can point to some of the options with those readers. I think there is one to force the data to be flipped inside the reader so that it doesn't get produced upside down.

Again, I'm just trying to do as little processing as possible and to plot the data as-is without resampling or interpolation, just with the geospatial metadata extracted from the native file.
Maybe the differences in the images are only due to the way imshow process the data; I don't know what is Scene.show() using under the hood to produce the plot.

[djhoese]

Scene.show() it using PIL (Pillow) .show() with a PNG version of the raw data array. It isn't doing any interpolation, but the pixels aren't interpolated or morphed to fit a particular size. What I mean by that is that your plot isn't the same size ratio as the original array so pixels are probably stretched vertically or horizontally.

What I suggested about the keyword arguments to the reader would not be any "real" additional processing. It would be mostly just flipping the data around. You are probably using the HRV channel which may require an additional keyword argument to produce a "padded" version of the HRV channel (HRV is made up of two separate images, but a single full disk is needed to use it without resampling). I think if you do Scene.load(..., upper_right_corner='NE') then the flipping might be fixed. If you create the Scene with reader_kwargs={"fill_disk": True} when using the seviri_l1b_native reader then that should fix the native resampler issue...I think.

[guidocioni]

Scene.show() it using PIL (Pillow) .show() with a PNG version of the raw data array. It isn't doing any interpolation, but the pixels aren't interpolated or morphed to fit a particular size. What I mean by that is that your plot isn't the same size ratio as the original array so pixels are probably stretched vertically or horizontally.

What I suggested about the keyword arguments to the reader would not be any "real" additional processing. It would be mostly just flipping the data around. You are probably using the HRV channel which may require an additional keyword argument to produce a "padded" version of the HRV channel (HRV is made up of two separate images, but a single full disk is needed to use it without resampling). I think if you do Scene.load(..., upper_right_corner='NE') then the flipping might be fixed. If you create the Scene with reader_kwargs={"fill_disk": True} when using the seviri_l1b_native reader then that should fix the native resampler issue...I think.

Thanks I'll definitely try that!
Maybe if I don't need to put many overlays it would make more sense to use the show method and forget about plotting with cartopy..may be even more efficient