AVHRR Band 2 image shows black where bright clouds should be

[kathys]

Liam informed us that he was running into cases where the AVHRR Band 2 images were missing the clouds (bright regions). It seemed to only happen in the band 2 data. Nigel Atkinson informed us that the slope after applying the coefficients was sometimes negative in band 2, and therefore an offset had to be applied to ensure that it was always positive. Dave investigated and found that the negative slope check and offset was included in the code already, but found that the section of code where it was applied "polar2grid/avhrr/readers.py" line 299, did not recognize it as band 2, but thought it was band 1.

Here is the image that Liam provided, and I was able to reproduce.

[djhoese]

I've asked @mraspaud and @adybbroe of the Pytroll group about this as I think they were the last authors on the Satpy reader. The Satpy reader also seems to be effected by this based on code review, haven't looked at actual imagery.

[mraspaud]

Looking at the code, the negative slope parameter is currently corrected for band 3a, not 2 at the moment. So it might be a bug if 3a doesn't need to be corrected at all.

[djhoese]

From Nigel Atkinson:

In the level 1b file, the visible coefficients are stored as 4-byte integers. Scaling factors then convert them to real numbers which are applied to the measured counts. The coefficient is different depending on whether the counts are less than or greater than the high-gain/low-gain transition value (nominally 500). The slope for visible channels should always be positive (reflectance increases with count). With the pre-launch coefficients the channel 2 slope is always positive but with the operational coefs the stored number in the high-reflectance regime overflows the maximum 2147483647, i.e. it is negative when interpreted as a signed integer. So you have to modify it. The AAPP utility “avhrrin”, which converts from l1b to l1c (reflectance/BT), does this test for the slope in the high-reflectance regime as follows (subroutine avh_lbc.F):

        do c = 1, avh_mxvischn
           if (avh_calvis(1,1,c) .NE. 0) then     !test the operational slope
             vcoefset=1
           else
             vcoefset=3
           endif
           slope     = avh_calvis(1,vcoefset,c)*1.e-10
           intercept = avh_calvis(2,vcoefset,c)*1.e-7
           do i4ncount = 0 , avh_calvis(5,vcoefset,c)
              tabalb(i4ncount,c) = (slope * i4ncount + intercept)
              if(tabalb(i4ncount,c).lt.0.) tabalb(i4ncount,c)=0.
              if(tabalb(i4ncount,c).gt.160.) tabalb(i4ncount,c) = 160.
           enddo

           slope     = avh_calvis(3,vcoefset,c)*1.e-10
           if (slope .lt. 0) slope = slope + 0.4294967296  !Ensure positive
           intercept = avh_calvis(4,vcoefset,c)*1.e-7
           i4nmax =  65535
           do i4ncount = avh_calvis(5,vcoefset,c), i4nmax
              tabalb(i4ncount,c) = (slope * i4ncount + intercept)
              if(tabalb(i4ncount,c).lt.0.) tabalb(i4ncount,c) = 0.
              if(tabalb(i4ncount,c).gt.160.) tabalb(i4ncount,c) = 160.
           enddo
        enddo

Based on that code though this is probably checked for every visible channel.

[mraspaud]

So would the solution be to use unsigned ints here instead?
https://github.com/pytroll/satpy/blob/4cef629e7093b70352f088e631dc3dfafbedde77/satpy/readers/aapp_l1b.py#L451

[djhoese]

@mraspaud Logically that makes sense (overflowing max signed int but want to use it as unsigned int) but will only work if that piece of the code only applies to visible channels. If it also applies to IR channels then it may break cases where slope is supposed to be negative.

The other workaround is to make it so that slope check in Satpy is always checked for visible channels.

[mraspaud]

@djhoese could you extract the relevant fields from the data header so that we can make a test case?

[djhoese]

You mean a test case for a unit test in Satpy? As a starting point @kathys can probably get us the files Liam used.

[mraspaud]

yes, sorry, I thought you had the file already.

[djhoese]

@mraspaud If you want to play with this, the file in question is on bumi at the SSEC here:

bumi:/data/users/kathys/liam_avhrr/data/hrpt_noaa18_20210327_1411_81698.l1b

If there's a better place to put it, let me know.

[mraspaud]

I don't seem to have access to bumi anymore... google drive or dropbox? Won't really have time now though until after my vacation.

[djhoese]

@mraspaud Did you go through "ash" to connect to it? If so, maybe your account expired. Anyway here is the file: ftp://ftp.ssec.wisc.edu/pub/ssec/davidh/hrpt_noaa18_20210327_1411_81698.l1b

[mraspaud]

Yes, I'm using ash as gateway. I can't seem to access the ftp either, sorry...

[mraspaud]

to be more exact, the ftp is fine, but the /pub/ssec/davidh/ doesn't seem to exist.

[djhoese]

That's because it gets cleared out after 7 days. I'll put it up again later today.

[djhoese]

@mraspaud The file should be there now.

[mraspaud]

got it, thanks

[djhoese]

@kathys I don't know if you did this already, but if you change that if statement we pointed out ("polar2grid/avhrr/readers.py" line 299,) to check for 1 instead of 2, do you get the right image?

[djhoese]

Closing this as it is now fixed in Satpy and tested by @kathys in current P2G.

[kathys]

This problem is still seen in Metopb and Metopc band 3a geotiff images.

[djhoese]

With my local P2G version this is what I get for the original use case of this PR:

But this issue was originally just for band 2. Luckily Nina Hakansson fixed it in Satpy for the other bands: pytroll/satpy#2123

But this PR was only merged on July 29th, after your latest beta tarball. I think we're good here, but let's put it on the TODO list to recheck with the new tarball.

[kathys]

The metopb and metopc input files:
bumi:/data/users/kathys/test_data/avhrr/2022_examples/hrpt_M01_20220923_1630_51969.l1b
bumi:/data/users/kathys/test_data/avhrr/2022_examples/hrpt_M03_20220918_1545_20055.l1b

Output files:
bumi:/data/users/kathys/test_data/avhrr/metopc/new/
bumi:/data/users/kathys/test_data/avhrr/metopb/new/

You are right, the band 2 data now looks fine. It is only the band 3a data that does not display correctly, and 3a is only available on the metop avhrr satellites.

[djhoese]

Processing with modern Satpy shows it is fixed:

Closing. Feel free to reopen if you'd rather wait until you have a chance to test it.