Merging GOES17, EUMETSAT, GK-2A Meteorological satellite image

Question

I want to generate global weather satellite image using GOES17, EUMETSAT, and GK-2A. I want make it Plate carree coordinate. (GOES 17 netcdf file convert to Plate Carree)

First, using Satpy, I made plate carree image.

from satpy import Scene
from glob import glob
from pyresample import create_area_def

area_def = create_area_def("my_area_def", "+proj=eqc +datum=WGS84", resolution=2000)
goes17 = glob('./samplefile/*')
goes17_scene = Scene(reader="abi_l1b", filenames=goes17)
goes17_scene.load('[C13]')
new_scn = goes17_scene.resample(area_def)

# save to geotiffs
new_scn.save_datasets()

like this method, I want make other satellite image and merging to 1 image file. but is there any simple or easiest way to generate global weather image? My final goal is generate numpy array of global satellite image.

-- My entire code --

from satpy import Scene, MultiScene
from glob import glob
from pyresample import create_area_def


area_def = create_area_def("my_area_def", "+proj=eqc +datum=WGS84", resolution=2000,)


goes17 = glob('E:/Global/GOES_17/OR_ABI-L1b-RadF-M6C13_G17_s20212130000319_e20212130009396_c20212130009445.nc')
goes17_scene = Scene(reader="abi_l1b", filenames=goes17)
goes17_scene.load(['C13'])

gk2a = glob('E:/Global/GK-2A/gk2a_ami_le1b_ir105_fd020ge_202108010000.nc')
gk2a_scene = Scene(reader="ami_l1b", filenames=gk2a)
gk2a_scene.load(['IR105'])

eumetsat = glob('E:/Global/EUMETSAT/MSG4-SEVI-MSG15-0100-NA-20210801000010.306000000Z-20210801001259-4774254.nat')
eumetsat_scene = Scene(reader='seviri_l1b_native', filenames=eumetsat)
eumetsat_scene.load(['IR_108'])


from satpy import MultiScene, DataQuery
mscn = MultiScene([goes17_scene, gk2a_scene, eumetsat_scene])


groups = {DataQuery(name='IR_group', wavelength=(10.35, 10.35, 10.8)): ['C13', 'IR105', 'IR_108']}
mscn.group(groups)


from pyresample.geometry import AreaDefinition
resampled = mscn.resample(area_def, reduce_data=False)

resampled.load(['IR_group'])
blended = resampled.blend() 
blended.show(['IR_group'])

Answer 1

There are some ways to do this inside Satpy, but typically people have specific ways they want the data joined together. That is a question you'll have to answer before you choose the code you want. First though you need to make a Scene for each separate satellite image you want in the final image and resample them to the same grid. A DynamicAreaDefinition (as you're using now) is not good for this overall process as each resampled Scene would be on a different final area (based on the satellite data being resampled that "froze" the DynamicAreaDefinition).

Your options for merging:

1. Satpy has a BackgroundCompositor where you can put one image on top of another. There is some documentation for creating a custom composite where you could make a composite like this. A series of these composites could be chained together to get the overall global composite you are looking for. You can put all the datasets in the same Scene to make things easier:

scn = Scene()
scn["C13"] = resampled_goes17_scene["C13"]
... and so on for the other sensors ...
scn.load(["custom_composite"])

2. Use the Satpy MultiScene, give it all of your resampled Scenes and run the "blend" method to join the images together. https://satpy.readthedocs.io/en/stable/multiscene.html#blending-scenes-in-multiscene
3. Use xarray and dask .where functions along with a custom mask array to say where each image should appear in the overall image. Some people do this type of thing with solar zenith angle have a nice blend between the images rather than just overlaying one on top of the other.
4. Create individual geotiffs for each resampled sensor and use GDAL's gdal_merge.py utility to join them into one geotiff.