Add a reader for FY-4A LMI level 2 data

[zxdawn]

This PR adds a reader for the FY-4A LMI (Lightning Mapping Imager) L2 data format which are NetCDF files.

• Tests added
• Tests passed
• Passes flake8 satpy
• Fully documented

[coveralls]

Coverage increased (+0.09%) to 89.428% when pulling cdba35f on zxdawn:lmi into bbfb44f on pytroll:master.

[codecov]

Codecov Report

Merging #1103 into master will increase coverage by 0.10%.
The diff coverage is 23.28%.

@@            Coverage Diff             @@
##           master    #1103      +/-   ##
==========================================
+ Coverage   89.34%   89.44%   +0.10%     
==========================================
  Files         195      201       +6     
  Lines       28785    29560     +775     
==========================================
+ Hits        25717    26439     +722     
- Misses       3068     3121      +53     

Impacted Files	Coverage Δ
satpy/readers/lmi_l2.py	23.28% <23.28%> (ø)
satpy/readers/seviri_base.py	90.90% <0.00%> (-7.96%)	⬇️
satpy/readers/tropomi_l2.py	89.32% <0.00%> (-4.80%)	⬇️
satpy/composites/__init__.py	78.76% <0.00%> (-0.21%)	⬇️
satpy/tests/test_scene.py	99.71% <0.00%> (-0.14%)	⬇️
satpy/readers/fci_l1c_fdhsi.py	96.15% <0.00%> (-0.03%)	⬇️
satpy/readers/__init__.py	95.09% <0.00%> (ø)
satpy/tests/test_composites.py	100.00% <0.00%> (ø)
satpy/tests/writer_tests/test_cf.py	98.53% <0.00%> (ø)
satpy/tests/reader_tests/test_utils.py	100.00% <0.00%> (ø)
... and 29 more

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update bbfb44f...cdba35f. Read the comment docs.

[zxdawn]

@djhoese I have to assign the y coords, although the units isn't meters. Otherwise, MultiScene won't work. Since there're LON and LAT variables, is the wrong units not matter?

BTW, if I assign the time coords which have the same dimension of y coords, the dimensions of y corrds at different times are still different. That's why MultiScene doesn't work.

[djhoese]

Why do you need to assign these? I don't think you need y coordinates here. Maybe I'm wrong and naming it y is a bad idea, but I think that is consistent with other readers and the rest of Satpy. The dimension should be y, but no coords for it. My thought was instead that you add a data.coords['time'] = [... series of times ...] then the MultiScene timeseries function could be updated to look at that coordinate. So time would be the only coordinate you add in the reader, the base reader will then add lon and lat for you.

[zxdawn]

@djhoese

Got it! Here's the method:

        tseries = pd.DatetimeIndex(np.repeat(data.attrs['start_time'], data.shape[0]))
        data = data.expand_dims({'time': tseries})
        data = data.rename({'x': 'y'})

Example:

<xarray.DataArray 'ER' (time: 31, y: 31)>
dask.array<where, shape=(31, 31), dtype=float32, chunksize=(31, 31), chunktype=numpy.ndarray>
Coordinates:
  * time     (time) datetime64[ns] 2019-07-25T05:25:10 ... 2019-07-25T05:25:10
Dimensions without coordinates: y

But, I will get this error when using MultiScene by mscn = MultiScene.from_files(filenames, reader='lmi_l2', group_keys=('start_time', 'subtask_num')) and imgs = mscn.blend(blend_function=timeseries):

ValueError: arguments without labels along dimension 'y' cannot be aligned because they have different dimension sizes: {8, 74, 44, 45, 86, 55, 24, 60, 93, 31}

The modification of MultiScene is mentioned in #1173.

[zxdawn]

Here's an example of plotting the number of events on 7800m resolution.
It looks summed values are correct.
But, I'm not sure how to create the correct 7800m AreaDefinition. I can't find the explanation of this resolution ... So, is it suitable to create the correct AreaDefinition based on the sub-satellite point and same information of FY4A AGRI?

Code

mscn = MultiScene.from_files(filenames, reader='lmi_l2', group_keys=('start_time', 'subtask_num'))
mscn.load(['ER'])

# resample to 7800m (resolution of LMI)
mscn = mscn.resample('china_7.8km', resampler='bucket_count')

# blend by time dim
imgs = mscn.blend(blend_function=timeseries)

# sum by time dim
sum_num = imgs['ER'].sum('time')

# copy attrs
sum_num.attrs = imgs['ER'].attrs

crs = sum_num.attrs['area'].to_cartopy_crs()
ax = plt.axes(projection=crs)
ax.coastlines()
ax.gridlines()
ax.set_global()
plt.imshow(sum_num, transform=crs, extent=crs.bounds, origin='upper')
cbar = plt.colorbar()

Figures

2019-07-25 05:15:10 to 05:34:49

2019-07-25 05:15:10 to 05:24:49

2019-07-25 05:25:10 to 05:34:49

[djhoese]

But, I'm not sure how to create the correct 7800m AreaDefinition. I can't find the explanation of this resolution ... So, is it suitable to create the correct AreaDefinition based on the sub-satellite point and same information of FY4A AGRI?

That seems reasonable. I think it depends on what you are doing and what you hope to get out of it.