**Memorandum for:** GOES Aerosol and Smoke Product (GASP) users

**From:** Shobha Kondragunta (<Shobha.Kondragunta@noaa.gov>) and Matthew
Seybold (<Matthew.Seybold@noaa.gov>), co-leads for GASP, NOAA/NESDIS

**Subject:** Revisions to GOES Aerosol and Smoke Product (GASP)
Algorithm

**Date:** January 7, 2009

**Summary**

This document describes the changes made to the GOES Aerosol and Smoke
Product (GASP) algorithm and output file. Changes include four parts:
(1) bug fix related to the definition of relative azimuth angle, (2)
procedure to calculate Aerosol Optical Depth (AOD) standard deviation,
(3) addition of scattering angle information to the output file, and (4)
new procedure to calculate surface reflectance. Details of the changes
and impact on the product are described below:

(1) **Relative azimuth angle definition**

An error related to the definition of relative azimuth angle in the GASP
code was identified. In the original code, the relative azimuth angle is
subtracted from 180^o^ before it is used as a search index in the lookup
table. The look-up table is, however, written in the reverse order of
relative azimuth angle. This resulted in wrong indexing in the look-up
table and led to errors in surface reflectance retrieval and AOD
retrieval. In the new code, this bug has been fixed. Figure 1 shows an
example of expected changes in the output. [For users of GASP binary
files or imagery, the change will be transparent.]{.mark} For users who
run the GASP code, the errors in the following subroutines need to be
fixed:

**GOES-WEST:**

> In the program named as "goes10.pro", both on line 2127 and 3551:
>
> Originally it is written as:
>
> **pazi = 180. -- pazi**
>
> Change made:
>
> **delete this line**

**GOES-EAST:**

> In program named as "goes12.pro", both on line 2158 and 3608:
>
> Originally it is written as:
>
> **pazi = 180. -- pazi**
>
> Change made:
>
> **delete this line**

**(2) AOD standard deviation**

AOD standard deviation is one of the parameters used to screen cloud
contaminated AOD values. In the IDL reader (Appendix A) used to read
GASP binary files, pixels with AOD standard deviations \> 0.3 are
flagged as bad data. The procedure to calculate AOD standard deviation
has been modified from the original method. In the original method, AOD
standard deviation was calculated in a 3 X 3 box for all pixels, without
checking if any of the 9 pixels in this 3 X 3 box are cloudy. Not
checking if a pixel is cloud-free or not resulted in high AOD standard
deviation, when some of the pixels in the 3 X 3 box are contaminated by
clouds. In the new method, AOD standard deviation is calculated using a
5 X 5 box but only for pixels which are identified as clear by the cloud
mask, and the number of valid (cloud-free) pixels in the box is greater
than 10. Otherwise AOD in this central pixel is considered as invalid.
[This change will also be transparent to the users.]{.mark}

Figure 2 shows an example of original operational GOES-12 AOD retrievals
compared to retrievals with new AOD standard deviation calculation.

**(3) Addition of scattering angle to the output file**

Research work has also shown that AOD retrievals are very sensitive to
scattering angle and for certain scattering angle range (backward
scattering), sensitivity to AOD changes is poor and thus retrievals are
not very accurate. Based on various tests including AERONET comparisons,
we noticed that AOD retrievals are accurate for scattering angles
between 50^o^ and 170^o^. Figure 3 shows an example of GASP output with
and without scattering angle based screening. However, this screening is
very tight and reduces the data coverage for certain locations and times
of the day. Instead of screening the data for users, we decided to
provide all data but also provide scattering angle information in the
output files so users can use their discretion in using the data.
Scattering angle is written out as a parameter in the GASP binary output
files that have an extension of '.all' in the filename. For historic
data files that do not have scattering angle information, it can be
easily computed using the code provided below (Appendix B). We recommend
not using AOD data for scattering angles less than 50^o^ and greater
than 170^o^.

**(4) Surface reflectance calculation**

In the original approach, clear sky composite is created using counts in
channel 1. For any given day, 2nd darkest pixel (in units of counts) in
the previous 28 days including current day, is used to compute surface
reflectance. In converting counts to surface reflectance, solar zenith
angle corresponding to the day for which AOD is being calculated is
used. This approach creates a large artifact in the early morning and
late evening since changes of solar zenith angle during the 28 days can
be large, especially in the Spring and the Fall. In the new approach,
surface reflectance is retrieved for every day assuming a background
aerosol optical depth of 0.02 and a search is made for the 2^nd^ darkest
pixel (in surface reflectance) in the previous 28 days including current
day. This improves the surface reflectance retrieval because the correct
solar zenith angle is used.

![](media/image1.jpeg){width="2.943700787401575in"
height="2.748425196850394in"}![](media/image2.jpeg){width="2.943700787401575in"
height="2.748425196850394in"}

Figure 1: GOES-12 original processing (left panel) and processing with
azimuth angle bug fix (right panel)

> ![](media/image3.jpeg){width="2.944488188976378in"
> height="2.748031496062992in"}![](media/image4.jpeg){width="2.944488188976378in"
> height="2.748031496062992in"}
>
> Figure 2: GOES-11 data with old AOD standard deviation (left panel)
> and new AOD standard deviation (right panel)
>
> ![](media/image5.jpeg){width="2.944488188976378in"
> height="2.748031496062992in"}![](media/image6.jpeg){width="2.944488188976378in"
> height="2.748031496062992in"}
>
> Figure 3: GOES-11 processing with no scattering angle based screening
> (left panel) and with scattering angle based screening (right panel)

**Appendix A: IDL Reader to read GASP binary file**

function rd_aod_all,bcld=bcld,bad=bad,bcls=bcls,glat=glat,glon=glon,\$

dsn=dsn,dir=dir,bsfc=bsfc,bmos=bmos,bmsk=bmsk,bch1=bch1,bsig=bsig,\$

baodstd=baodstd,bsca=bsca,nsat=nsat,nx=nx,ny=ny,good=good

nfiles=1

for i=0,nfiles-1 do begin

print,nx,ny

; f1=dir+dsn

f1=dsn

;read in data

baod = bytarr(nx,ny)

bcls = bytarr(nx,ny)

baodstd=baod

bsfc = baod

bch1 = bsfc

bmos = bsfc

bcld = bsfc

bsig = bsfc

bmsk = bsfc

bcls = bsfc

bsca = bsfc

print,f1

openr,ilun,f1,/comp,/get_lun

readu,ilun,baod,bmsk,bcls,baodstd,bsfc,bch1,bmos,bcld,bsig,bsca

;check,baod

free_lun,ilun

cc=where(baod eq 0)

dd=where(baod eq 255)

baod = float(baod) / 100. - 0.5

; printf,2,baod

; close,2

baodstd = float(baodstd) / 100.

bsfc = float(bsfc) / 500. - .1

bch1 = float(bch1) / 600.

bmos = float(bmos) / 600.

bsig = float(bsig) / 250.- 0.5

bsca = float(bsca)

; newmask = bytarr(nx,ny)

; Note: here selection criteria may change based on user preference

; Normal criteria

good = where(baodstd lt 0.3 and bsig gt 0.01 and bsfc lt 0.15 and bsfc
gt 0.005 \$

and bcls gt 15 and baod lt 10.and bch1 gt 0. and bcld eq 1 and bsca gt
70 and bsca lt 170, \$

ngood,compl=bad,ncomp=nbad)

; good = where(bsfc lt 0.15 and bsfc gt 0.005 \$

; and baod lt 10.and bch1 gt 0. and bcld eq 1, \$

; ngood,compl=bad,ncomp=nbad)

; A little strict criteria

; good = where(baodstd lt 0.15 and bsig gt 0.01 and bsfc lt 0.15 and
bsfc gt 0.005 \$

; and bcls ge 25 and baod lt 10.and bch1 gt 0. and bcld eq 1, \$

; ngood,compl=bad,ncomp=nbad)

baod(bad)=-999.9

bch1(bad)=-999.9

;bsfc(bad)=-999.9

;bcld(bad)=-999.9

;bcls(bad)=-999.9

endfor

return,baod

stop

end

**Appendix B: Subroutine to calculate the scattering angle**

An IDL subroutine to calculate the viewing geometry for both GOES-EAST
and WEST satellite is attached. Besides the scattering angle, this
subroutine also calculates solar zenith angle, solar azimuth angle,
satellite viewing zenith angle, satellite azimuth angle and relative
(solar - satellite) azimuth angle. The required inputs for this
subroutines are as follow:

> **jday -- the day of the year**
>
> **time -- UTC time of the observation in the format of "HHMMSS". HH is
> the hour, MM is the minutes and SS is the second.**
>
> **npixel -- number of pixel for each scanline:**
>
> **GOES-EAST: 2000**
>
> **GOES-WEST: 2500**
>
> **nlines -- number of scanline**
>
> **GOES-EAST: 850**
>
> **GOES-WEST: 912**
>
> **Sat_lat -- latitude of subsatellite point**
>
> **GOES-EAST: 0**
>
> **GOES-WEST: 0**
>
> **Sat_lon -- longitude of subsatellite point**
>
> **GOES-EAST: -75**
>
> **GOES-WEST: -135**
>
> **lat_geo -- array of latitude**
>
> **lon_geo -- array of longitude**
>
> The corresponding outputs are listed as follow:
>
> **solzen -- solar zenith angle**
>
> **az_not -- solar azimuth angle**
>
> **theta -- satellite viewing zenith angle**
>
> **azimuth -- satellite azimuth angle**
>
> **sca_ang -- scattering angle**
>
> **imgazi -- relative azimuth angle**
>
> **the relative azimuth angle is defined as solar azimuth angle minus
> satellite azimuth angle and in the range of 0 to 180.**
>
> **pro** **goes_gasp_geometry, \$
> jday,time,solzen,az_not,theta,azimuth,sca_ang,imgazi,lat_geo,lon_geo,npixel,nlines,
> \$**
>
> **sat_lat, sat_lon**
>
> **; ========input============**
>
> ; jday -- the day of the year
>
> ; time -- UTC time of the observation in the format of "HHMMSS". HH is
> the hour,
>
> ; MM is the minutes and SS is the second.
>
> ; npixel -- number of pixel for each scanline: 2000 (GOES-EAST), 2500
> (GOES-WEST)
>
> ; nlines -- number of scanline. 850 (GOES-EAST), 800 (GOES-WEST)
>
> ; lat_geo -- array of latitude
>
> ; lon_geo -- array of longitude
>
> **; ========output============**
>
> **;** solzen -- solar zenith angle
>
> ; az_not -- solar azimuth angle
>
> ; theta -- satellite viewing zenith angle
>
> **;** azimuth -- satellite azimuth angle
>
> **;** sca_ang -- scattering angle
>
> **;** imgazi -- relative azimuth angle
>
> imgnx=npixel
>
> imgny=nlines
>
> pi = !pi
>
> dtor = !pi / 180.
>
> lat = fltarr(imgnx,imgny)
>
> lon = fltarr(imgnx,imgny)
>
> lat=lat_geo
>
> lon=lon_geo
>
> rlat=0
>
> rlon = 0
>
> satlat = sat_lat
>
> satlon = sat_lon
>
> ;\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\--
>
> ; calculate solar angle :
>
> ;! solzen - solar zenith angle
>
> ;! az_not - solar azimuth angle
>
> HH = fix( FLOAT(TIME) / 10000.0)
>
> MM = fix((FLOAT(TIME)/100) - FLOAT(HH) \* 100)
>
> SS = TIME - HH \* 10000L - MM \* 100L
>
> tu = FLOAT(HH) + FLOAT(MM)/60.0 + FLOAT(SS)/3600.0
>
> pi = !pi
>
> rad = !pi / 180.
>
> glat = satlat
>
> glon = satlon
>
> HH = fix( FLOAT(TIME) / 10000.0)
>
> MM = fix((FLOAT(TIME)/100) - FLOAT(HH) \* 100)
>
> SS = TIME - HH \* 10000l - MM \* 100l
>
> DTIME = FLOAT(HH) + FLOAT(MM)/60.0 + FLOAT(SS)/3600.0
>
> H_ANG = -1.0 \* PI \* (DTIME-12.0) / 12.0
>
> nc = 0
>
> nd = 0
>
> DECLIN = -23.45 \* COS(FLOAT(jday)\*2.\*PI/365. + 10.\*2.\*PI/365.)
>
> na = where(lat ge -90. and lat le 90. and lat gt declin,na,\$
>
> complement=b,ncomplement=nb)
>
> gd = where(lat ge -90. and lat le 90.) ;,complement=ng)
>
> g=gd
>
> betaa=lat
>
> azimuth=lat
>
> betaa(\*)=0.
>
> mu_not = lat
>
> mu_not(\*,\*)=0.
>
> gd = where(lat ge -90. and lat le 90.) ;,complement=ng)
>
> ng = where(lat lt -90. or lat gt 90.) ;,complement=ng)
>
> MU_NOT(gd) = (SIN(LAT(gd)\*RAD)\*SIN(DECLIN\*RAD) + \$
>
> COS(LAT(gd)\*RAD)\*COS(DECLIN\*RAD)\*COS(H_ANG-LON(gd)\*RAD))
>
> solzen = acos(mu_not \< 1.) / rad
>
> mu_not = 0
>
> nc = 0
>
> nd = 0
>
> a = 0
>
> ; CALCULAT SATELLITE ZENITH ANGLE
>
> betaa=lat
>
> theta=lat
>
> betaa(\*)=0.
>
> betaa(g) = ACOS(COS(LAT(g)\*RAD) \* COS((GLON - LON(g))\*RAD))
>
> theta(g) = asin((42164.\*sin(betaa(g)) / \$
>
> sqrt(1.8084e9 - 5.3725e8\*cos(betaa(g)))) \< 1) / rad
>
> ;! CALCULATE SOLAR AZIMUTH (gasp)
>
> a = where(lat ge -90. and lat le 90. and lat gt declin,na,\$
>
> complement=b,ncomplement=nb)
>
> g=gd
>
> gd = 0
>
> az_not = lat
>
> az_not(\*,\*) = 0.
>
> if (na gt 0) then az_not(a) = ACOS( COS((LAT(a)-DECLIN)\*RAD) \*
> COS(H_ANG-LON(a)\*RAD))
>
> if (nb gt 0) then az_not(b) = ACOS( COS((LAT(b)-DECLIN)\*RAD) \* \$
>
> COS(LON(b)\*RAD-H_ANG))
>
> if (na gt 0) then AZ_NOT(a) = SIN(H_ANG-LON(a)\*RAD) / SIN(az_not(a))
>
> if (nb gt 0) then AZ_NOT(b) = SIN(LON(b)\*RAD-H_ANG) / SIN(az_not(b))
>
> az_not(g) = (az_not(g) \> (-1.)) \< 1.
>
> AZ_NOT(g) = ASIN( AZ_NOT(g) )
>
> if (na gt 0) then AZ_NOT(a) = PI - AZ_NOT(a)
>
> if (nb gt 0) then AZ_NOT(b) = PI + AZ_NOT(b)
>
> AZ_NOT(g) = AZ_NOT(g) /RAD
>
> nc = 0
>
> nd = 0
>
> a = 0
>
> if (nb gt 0) then c = where(az_not(b) lt
> 180.,nc,complement=d,ncomplement=nd)
>
> if (nc gt 0) then az_not(b(c)) = 180. - az_not(b(c))
>
> if (nd gt 0) then az_not(b(d)) = 540. - az_not(b(d))
>
> ;! CALCULATE THE SATELLITE AZIMUTH (GASP)
>
> betaa=lat
>
> betaa(\*)=0.
>
> betaa(g) = ACOS(COS(LAT(g)\*RAD) \* COS((GLON - LON(g))\*RAD))
>
> a = where(lat gt (-90) and lat le 90. and lat gt glat,na,\$
>
> complement=b,ncomplement=nb)
>
> azimuth=lat
>
> azimuth(\*,\*)=0.
>
> if (na gt 0) then AZIMUTH(a) = SIN((GLON-LON(a))\*RAD) / SIN(BETAa(a))
>
> if (nb gt 0) then AZIMUTH(b) = SIN((LON(b)-GLON)\*RAD) / SIN(BETAa(b))
>
> azimuth(g) = (azimuth(g) \> (-1.)) \< 1.
>
> azimuth(g) = asin(azimuth(g))
>
> if (na gt 0) then azimuth(a) = pi - azimuth(a)
>
> a = 0
>
> if (b(0) ge 0) then azimuth(b) = pi + azimuth(b)
>
> azimuth(g) = azimuth(g) / rad
>
> if (nb gt 0) then c = where(azimuth(b) lt 180.,nc)
>
> if (nb gt 0) then d = where(azimuth(b) ge 180.,nd)
>
> if (nb gt 0 and nc gt 0) then azimuth(b(c)) = 180. - azimuth(b(c))
>
> if (nb gt 0 and nd gt 0) then azimuth(b(d)) = 540. - azimuth(b(d))
>
> imgsza=solzen
>
> imgvza=theta
>
> imgazi=(az_not-azimuth)
>
> idx1=where (imgazi lt 0.,n1)
>
> if (n1 gt 0 ) then imgazi(idx1)=imgazi(idx1)+2.\*pi
>
> idx2=where (imgazi gt (2.\*pi),n2)
>
> if (n2 gt 0 ) then imgazi(idx2)=imgazi(idx2)-2.\*pi
>
> rad = dtor
>
> cosga =-cos(imgsza\*rad)\*cos(imgvza\*rad)
> -sin(imgsza\*rad)\*sin(imgvza\*rad)\*cos(imgazi\*rad)
>
> cosga = (cosga \> (-1.)) \< 1.
>
> sca_ang=acos(cosga) / rad

end