Improving resolution when setting extent

[ustropics]

I am having an issue on exactly what would give me 'true' resolution of a geostationary data set (GOES in particular) while allowing for conversion to cartopy. I have provided minimal code to reproduce below, but what would give me the highest resolution when setting a more 'zoomed' in extent? Is this where I need to use native resampling or are there some attributes I'm missing?

area_def = create_area_def(
    'composite_area',  # Name of the area definition
    proj_dict,
    units='degrees',
    width=5240, 
    height=5240,
    area_extent=[lon_min, lat_min, lon_max, lat_max]
)

scn.load([recipe])
new_scn = scn.resample(area_def, resampler='nearest', cache_dir='cache/')

area = new_scn[recipe].attrs['area']
        
plt_scn = scn.resample(area)

image = get_enhanced_image(plt_scn[recipe]).data
crs = plt_scn[recipe].attrs['area'].to_cartopy_crs()

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

ax.add_feature(cfeature.COASTLINE, edgecolor=data['border_color_val'], linewidth=float(data['border_width_val']))
ax.add_feature(cfeature.BORDERS, edgecolor=data['border_color_val'], linewidth=float(data['border_width_val']))
ax.add_feature(cfeature.STATES, edgecolor=data['border_color_val'], linewidth=float(data['border_width_val']))
ax.gridlines(draw_labels=True, linewidth=0.5, color='gray', alpha=0.5, linestyle='--')

if len(plt_scn[recipe].shape) > 2:
    print("ax: ", ax)
    image.plot.imshow(vmin=0, vmax=1, add_colorbar=False, rgb='bands', ax=ax)
else:  # Single-band grayscale case
    image = np.squeeze(image)
    image.plot.imshow(vmin=0, vmax=1, add_colorbar=False, ax=ax)

plt.savefig(filename, dpi=300, bbox_inches='tight')
plt.close()

In this particular image below, I've set the extent for Florida (info below) and the recipe/composite to cimss_true_color_sunz_rayleigh:

  {
    "location": "Florida",
    "extent": [-87.634, -80.031, 24.396, 31.0],
    "satellite": "GOES-EAST (16)",
    "domain": "C",
    "projection": "Transverse Mercator",
    "central_longitude": -83.83250000000001,
    "central_latitude": 27.698,
    "standard_parallel": [24.396, 31.0]
  }

When cross-comparing this to products I've created in goes2go and to NESDIS satellite output, this appears to not match the resolution for the green/blue abi-bands (3 and 1). I'm confident this is just a lack of my understanding on how to sample properly or missing an attribute, but where am I going wrong?

My output image:

[djhoese]

Do you have an example image from those other tools? Or can you provide more details about what's signaling a lower resolution when you look at this image? I assume the size of your area is similar to the size of the images you're comparing against?

[guidocioni]

I'm replying because you use an approach that is quite similar to mine and I had this issue in the past.

1. Did you try to provide interpolation=none when calling imshow?
2. Also I think that you don't need to specify the area when resampling. You can just do

scn = satpy.Scene(**scene_inputs)
scn.load(['ash'], upper_right_corner='NE')
scn = scn.crop(ll_bbox=tuple(cust_area))
scn = scn.resample(resampler='native')

which should hopefully choose the highest possible resolution.

[ustropics]

Apologies on the late reply (university life before Thanksgiving break is always hectic 😛), but as @djhoese and @guidocioni essentially indicated there wasn't really an issue with the code perse. The actual issue was the dataset I was importing to create the image. The TL;DR if anyone has similar problems in the future, don't get the data from the ABI-L2-MCMIP product.

The ABI L2 Cloud and Moisture Imagery has two slightly different products:
ABI-L2-MCMIP – multichannel format where all bands are in one file. The resolution of all channels is reduced to 2000m, which is the root of the resolution issue.
ABI-L2-CMIP – full resolution data (red = 500 m, VIS/NIR = 1000 m, TIR = 2000 m); every band has its own file because the resolutions differs.

My new method is to check the composite the user is selecting first, then pull in the individual bands from the ABI instrument radiance product (ABI-L1B-Rad) and using glob to pull in the filenames. Satpy will automatically detect what composites that are available based on the bands you pull in (can use scn.available_composite_names() to check).