FY3D MERSI-2 issues

[kathys]

Working with FY3D MERSI2 images using P2G V3.0 I encountered these issues:

• [x] Nearest neighbor resampling is being used when it should be EWA with D=40 and d=1.
• [x] There is an extra default image being created, natural_color.
• [x] There is a discolored border visible along the left side of the true color image that is not present on Version 2.3 images:

Data files: bumi:/data/users/kathys/test_data/fy3d_mersi/day/l1b

Images and log files: v3.0: bumi:/data/users/kathys/test_data/fy3d_mersi/day/geotiff_v3 v2.3: bumi:/data/users/kathys/test_data/fy3d_mersi/day/geotiff

[djhoese]

The natural color RGB seems to have been accidentally added by @joleenf in https://github.com/ssec/polar2grid/pull/359. I'm not sure if that was copied from ABI or something else. It also seems deliberate although I can't find the discussion since it is in the docstring at the top of the module.

[djhoese]

Nevermind, it was actually me who added it to the docs. Looks like the old "natural_color" was there for a long time (2019 when the reader was added), but I think we wanted it to be more like VIIRS so we switch to "false_color". Instead of completely replacing it we had both available but natural_color was never a default. In Joleen's #359 it was added as a default.

So I think it gets documented, but not a default. I'll fix that.

[djhoese]

@kathys Could you try forcing the EWA resampling with --method ewa --weight-delta-max 40 --weight-distance-max 1 and see if that changes that left-most column?

[kathys]

@djhoese Done. It did not make a difference to the left-most column.

[djhoese]

I was going to say if the resampling wasn't wrong in the first case that I would think this is something in the data that was masked by the bugs in EWA that were fixed in Polar2Grid 3.0. I guess I'll have to dive in.

[djhoese]

I know you had concerns about RSRs in pyspectral, it looks like something might need to get fixed in pyspectral because when running locally I get this:

WARNING  : rsr data may not be up to date: /home/davidh/.local/share/pyspectral/rsr_mersi-2_FY-3D.h5
4811182it [00:01, 3592830.47it/s]
WARNING  : Could not get the reflectance correction using band name: 1
WARNING  : Will try use the wavelength, however, this may be ambiguous!
WARNING  : A wavelength is provided instead of band name - disregard the relative spectral responses and assume it is the effective wavelength: 0.470000 (micro meter)
WARNING  : Could not get the reflectance correction using band name: 3
WARNING  : Will try use the wavelength, however, this may be ambiguous!
WARNING  : A wavelength is provided instead of band name - disregard the relative spectral responses and assume it is the effective wavelength: 0.650000 (micro meter)
WARNING  : Could not get the reflectance correction using band name: 2
WARNING  : Will try use the wavelength, however, this may be ambiguous!
WARNING  : A wavelength is provided instead of band name - disregard the relative spectral responses and assume it is the effective wavelength: 0.550000 (micro meter)

This usually means pyspectral doesn't know how to map satpy's name for the channel (ex. "1") to the name in the file so it uses the wavelength as an estimate. Probably the same results, but it should be fixed.

[djhoese]

Something is weird with this data. Here is edge of band 20:

And each of the bands seems to be a slightly different size/width. I'll check on the raw data to see if it is like this.

[djhoese]

The 8 first columns of band 2 seem to be bad:

But that doesn't explain why this wouldn't be completely masked off since that is what the ratio sharpening is supposed to do.

[djhoese]

Here's band 3 (top), 2, and 1 (bottom), on the left edge in the middle of row of the swath:

[djhoese]

Ah the bands are all the same resolution! So there is no ratio sharpening and therefore no "this band is NaN so all the bands should be NaN".

[djhoese]

Still debugging this, but discovered something that explains some of the confusion I was having. If I use pure Satpy I can reproduce this if I do:

1. Load
2. Resample lcc_fit using EWA
3. Save to geotiff

But if I instead:

1. Load
2. Resample to the finest resolution using the "native" resampler
3. Resample lcc_fit using EWA
4. Save to geotiff

Then I don't see the issue.

[djhoese]

Ok, so the main issue seems to be angle datasets. The angles (solar/satellite azimuth/zenith) are all at 1000m resolution while the band data is at 250m. Looking at the right-side of the lon/lat arrays (which is geographically the left-side) I see:

DataID(name='longitude', resolution=250, modifiers=())
[[-78.49014  -78.50028  -78.51044  -78.5206   -78.530785 -78.54098  -78.5512   -78.561424 -78.57166  -78.581924]
 [-78.4917   -78.50185  -78.512    -78.52217  -78.53236  -78.54255  -78.552765 -78.562996 -78.57324  -78.58349 ]
 [-78.49326  -78.50341  -78.513565 -78.523735 -78.53392  -78.54412  -78.55434  -78.56457  -78.574814 -78.58507 ]

DataID(name='longitude', resolution=1000, modifiers=())
[[-78.17625  -78.21511  -78.254196 -78.293495 -78.33301  -78.37275  -78.41272  -78.45291  -78.49333  -78.534   ]
 [-78.18235  -78.22123  -78.26034  -78.29965  -78.33919  -78.378944 -78.41894  -78.459145 -78.499596 -78.54027 ]
 [-78.18847  -78.22737  -78.26649  -78.305824 -78.34538  -78.38516  -78.42517  -78.4654   -78.50587  -78.54657 ]

Looking at the right-most pixels, I see that in the 250m longitude data we have ~5 pixels before we even get to similar longitude values in the 1000m longitude data. When I looked at the opposite side of the scene (the eastern most side) I only see ~1 pixel extra of the 1000m data, but this is expected (in my opinion) since the foot print of the 1000m pixels should be larger and "further out" than the 250m data.

The point is is that it seems that the 1000m data does not go out as far on the western side of this data as the 250m data by a decent number of pixels. This leads to resampling having out-of-swath fill pixels for the angle datasets where it has valid data pixels for the band datasets. The rayleigh correction ends up having a correction of 0 (since it is based on the angles) and that's why we get this strip of non-corrected swath.

Now to figure out why this didn't happen in the P2G v2.3.

[djhoese]

Figured it out. Here is the v3.0 version:

And v2.3:

The blue band is additional data that is not available in the v2.3 image. The reason it is blue is that it is not getting rayleigh corrected. The reason it isn't getting corrected is because the angle data does not exist there and results in a correction of 0. The reason it has a value of 0 instead of NaN (like in the v2.3 version) due to a change from @adybbroe here:

https://github.com/pytroll/pyspectral/blob/e2744cd6b54d28f7bc19c7637f9c7218d7b4247f/pyspectral/rayleigh.py#L270

This change follows the logic that "if we can't compute a correction then don't correct - do a correction of 0", rather than the old logic which is more like "if we can't compute a correction then this data is invalid in some way, mask it out".

[djhoese]

This change is being discussed on the pytroll slack.

[djhoese]

A fix for pyspectral's 0'ing out of the correction array is waiting approval and merge here: https://github.com/pytroll/pyspectral/pull/159

[djhoese]

The above mentioned pyspectral PR has been merged and released. The fix will be included in the next tarball.