Plotting Scene with Cartopy gives correct borders but red background

[christianversloot]

First things first: thanks for making such a well-supported library for satellite image processing.

Secondly: I'm aware that my question is Cartopy-oriented, but I'm opening an issue here because it is grounded in code one of Satpy's maintainers wrote in another issue (#1714). More specifically: #1714 (comment).

The problem: I am using Satpy to resample GOES-16 imagery using an area definition spanning the Texas/Oklahoma/New Mexico border area, using the true_color composite. I am trying to add borders with Cartopy. After resampling, I do get what I want:

Then, when using the resampled Scene as a Cartopy CRS, I do see the borders - but the background is red:

I have tried to tinker a lot, using other examples as well, but I have no clue what I'm doing wrong. The code posted in the other issue (although slightly older than a year) seems to have worked, so why isn't mine :) What am I missing?

Thanks for your help!

import os
from satpy import Scene
from glob import glob
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from pyresample import create_area_def

# Print channel files used
BASE_DIR = os.path.join(os.getcwd(), '23g16')
fn = glob(os.path.join(BASE_DIR, '*s20221432301174*.nc'))
for f in fn:
	print(f)

# Create scene and define composite
scn = Scene(filenames=fn, reader='abi_l1b')
composite = 'true_color'
scn.load([composite])

# Create area definition
ext = [-110.629439, 31.594328, -95.933875, 37.384890]
width = 3840
height = 2160
area_def = create_area_def('TX/OK/NM',
                           {'proj': 'longlat', 'datum': 'WGS84'},
                            area_extent=ext,
                            shape=(height, width),
                            units='degrees',
                            description='TX/OK/NM border')

# Resample to AreaDefinition
new_scn = scn.resample(area_def)
new_scn.load([composite])

# Save without Cartopy
new_scn.save_dataset(composite, filename='no_cartopy.png')

# Plot composite
CRS = new_scn[composite].attrs['area'].to_cartopy_crs()
fig = plt.figure(figsize=(30,25))
ax = fig.add_subplot(1, 1, 1, projection=CRS)
new_scn[composite].plot.imshow(rgb='bands', transform=CRS, origin='upper')
ax.add_feature(ccrs.cartopy.feature.STATES, linewidth=0.25)

# Save plot
plt.savefig('example.png')

[djhoese]

Try passing ax=ax to the .plot.imshow(...) call. Just a wild guess.

[christianversloot]

Thanks for your quick response! Unfortunately:

new_scn[composite].plot.imshow(rgb='bands', ax=ax, transform=CRS, origin='upper')

When I see the CRS description in the plot title, the code seems to be unknown. Although it seems odd to be related (as the borders are displayed correctly), it may be doing something to the Scene?

[mraspaud]

I just tried something similar with MSG data, and it looks like the data passed to matplotlib is to big. (the enhancements we usually apply to imagery before saving is not applied when calling .plot.imshow as far as I know). So maybe
(new_scn[composite] / 255).plot.imshow(rgb='bands', ax=ax, transform=CRS, origin='upper') would show something more sensible?

[mraspaud]

Ah, I got something better with this:

img = get_enhanced_image(new_scn[composite])
img.data.plot.imshow(rgb="bands", transform=CRS, origin="upper")
plt.show()

Not as straightforward as I would have wished, but seems to work.

[christianversloot]

Thanks Martin! This gives me a gray image:

import os
from satpy import Scene
from glob import glob
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from pyresample import create_area_def
from satpy.writers import get_enhanced_image

# Print channel files used
BASE_DIR = os.path.join(os.getcwd(), '23g16')
fn = glob(os.path.join(BASE_DIR, '*s20221432301174*.nc'))
for f in fn:
	print(f)

# Create scene and define composite
scn = Scene(filenames=fn, reader='abi_l1b')
composite = 'true_color'
scn.load([composite])

# Create area definition
ext = [-110.629439, 31.594328, -95.933875, 37.384890]
width = 3840
height = 2160
area_def = create_area_def('TX/OK/NM',
                           {'proj': 'longlat', 'datum': 'WGS84'},
                            area_extent=ext,
                            shape=(height, width),
                            units='degrees',
                            description='TX/OK/NM border')

# Resample to AreaDefinition
new_scn = scn.resample(area_def)
new_scn.load([composite])

# Save without Cartopy
# new_scn.save_dataset(composite, filename='no_cartopy.png')

# Plot composite
CRS = new_scn[composite].attrs['area'].to_cartopy_crs()
fig = plt.figure(figsize=(30,25))
ax = fig.add_subplot(1, 1, 1, projection=CRS)
img = get_enhanced_image(new_scn[composite])
img.data.plot.imshow(rgb="bands", transform=CRS, origin="upper")
ax.add_feature(ccrs.cartopy.feature.STATES, linewidth=0.25)

# Save plot
plt.savefig('example2.png')

[djhoese]

@christianversloot Your script with ax=ax in the imshow call works for me:

import os
from satpy import Scene
from glob import glob
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from pyresample import create_area_def

# Print channel files used
fn = glob("/data/satellite/abi/20170920/*.nc")
print(fn)

# Create scene and define composite
scn = Scene(filenames=fn, reader='abi_l1b')
composite = 'true_color'
scn.load([composite])

# Create area definition
ext = [-110.629439, 31.594328, -95.933875, 37.384890]
width = 3840
height = 2160
area_def = create_area_def('TX/OK/NM',
                           {'proj': 'longlat', 'datum': 'WGS84'},
                            area_extent=ext,
                            shape=(height, width),
                            units='degrees',
                            description='TX/OK/NM border')

# Resample to AreaDefinition
new_scn = scn.resample(area_def)
new_scn.load([composite])

# Save without Cartopy
new_scn.save_dataset(composite, filename='no_cartopy.png')

# Plot composite
CRS = new_scn[composite].attrs['area'].to_cartopy_crs()
fig = plt.figure(figsize=(30,25))
ax = fig.add_subplot(1, 1, 1, projection=CRS)
from satpy.writers import get_enhanced_image
img = get_enhanced_image(new_scn[composite])
img.data.plot.imshow(rgb='bands', transform=CRS, origin='upper', ax=ax)
ax.add_feature(ccrs.cartopy.feature.STATES, linewidth=0.25)

# Save plot
plt.savefig('example.png')

Note: This was with some old GOES-16 data I had in the "TEST" position but both your data and mine should cover Texas just fine. I also should note your second new_scn.load should not be needed.

[christianversloot]

I have it working now! For reference of anyone else running into this problem: ensure that your version of cartopy is up to date, which itself requires some hassle as it requires PROJ >8.0.0 but apt has ~6.3.7 as latest version. Building PROJ from source (https://stackoverflow.com/a/72442088) and then installing cartopy and then ensuring that the rest (satpy, pyresample,...) is up to date fixed it for me.

Thanks both for your help!

[djhoese]

Glad you got it working. My recommendation for anyone working with cartopy would be to use conda and install packages from conda-forge. Mixing apt and PROJ/cartopy is going to be difficult.