can't create ABI geo_color composite

[wbresky]

Describe the bug A clear and concise description of what the bug is. I am using the satpy ABI reader to create color composites of the CONUS sector. I can successfully create true_color and natural_color composites but when I try to create a geo_color composite it fails. To Reproduce

import numpy as np
from satpy import Scene
from glob import glob
import matplotlib.pyplot as plt
import cartopy
from satpy.utils import debug_on; debug_on()

cmp='geo_color'

# Get available G16 full disk file names
filenames = glob("./OR_ABI-L1b-RadF-M6C*_G16*.nc")

# create a "scene" object using ABI L1b reader
FD_scene = Scene(reader="abi_l1b",filenames=filenames)

# load ABI channels needed for color composites
FD_scene.load(['C01','C02','C03','C05','C07','C13'])

# resample to CONUS sector and create composite
new_scn = FD_scene.resample('goes_east_abi_c_500m')
new_scn.load([cmp])

# references for code below:
#       https://fire.trainhub.eumetsat.int/docs/figure1_MODIS_L1B.html
#       https://python-kurs.github.io/sommersemester_2019/units/S01E07.html

# transpose the color composite values to a shape that can be interpreted by the imshow method: (M,N,3)
image = np.asarray(new_scn[cmp]).transpose(1,2,0)

# scale the values to the range between 0 and 1, clipping the lower and upper percentiles
image = np.nan_to_num(image)
image = np.interp(image, (np.percentile(image,1), np.percentile(image,99)), (0, 1))

# convert scene to a cartopy CRS projection
crs = new_scn[cmp].attrs["area"].to_cartopy_crs()

# define axis and plot with matplotlib
ax = plt.axes(projection=crs)
ax.coastlines()
ax.gridlines()
ax.set_global()
plt.imshow(image, transform=crs, extent=crs.bounds, origin='upper')
plt.show()

Expected behavior

To be able to display a geo_color composite image of ABI CONUS sector

Actual results [DEBUG: 2024-11-08 19:15:50 : satpy.scene] Delaying generation of DataID(name='geo_color_low_clouds') because of incompatible areas

The error above causes a cascade of errors resulting in the error message:

Traceback (most recent call last): File "/data/home004/wayne.bresky/miniconda3/envs/test_env/lib/python3.11/site-packages/satpy/dataset/data_dict.py", line 169, in getitem return super(DatasetDict, self).getitem(item) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ KeyError: 'geo_color'

Screenshots If applicable, add screenshots to help explain your problem.

Environment Info:

• OS: Linux
• Satpy Version: 0.52.1
• PyResample Version: 1.31.0
• Readers and writers dependencies (when relevant): [run from satpy.utils import check_satpy; check_satpy()]

Additional context

I can create other color composites (true_color, natural_color, etc) with the same ABI data and python program but the more complicated geo_color composite is failing even though I have resampled the l1b data to the same projection (i.e., CONUS 500m resolution).

[pnuu]

See also https://github.com/pytroll/satpy/issues/2733

My workaround is to pre-resample the background images to the target area/projection and use those in my composite configs.

[wbresky]

Thanks for the quick response. I did see this bug report but I wasn't sure if it applied to my case.

I'm fairly new to satpy. Are you suggesting I create local copies of the static background images resampled to the same CONUS projection and use those images instead?

Could you possibly send me an example from your code?

Thanks.

Wayne

[strandgren]

In your example it should work if you simply move the loading of geo_color to the first load call, i.e. before the resampling. geo_color uses two static data files, which by default are on a global grid, causing the area incompatibility error you get.

If you're epeatedly generating the composite for the same area, it's indeed a good idea to resample the two static images to your target grid and use them instead. If you ping me on Monday I can provide an example, unless Panu is faster :)

[pnuu]

This is what I used to resample the required background images:

target_area_def_name = "EPSG_3035_1km"
fname = "gshhs_land_water_mask_3km_i.tif"
scn = Scene(reader='generic_image', filenames=[fname])
scn.load(["image"])
scn2 = scn.resample(target_area_def_name, resampler="gradient_search", method="nearest_neighbour")
scn2.save_dataset("image", filename="gshhs_land_water_mask_EPSG_3035_1km.tif", enhance=False, dtype="uint8")

And same for BlackMarble_2016_3km_geo.tif.

The original URIs are here and here.

So after resampling I've placed these images in the cache directory (or some https location), copied the composite definitions to my $SATPY_CONFIG_PATH/composites/<instrument>.yaml file, replaced the url (and hash, you can also just remove the hash lines) and then just used geo_color composite in the normal way.