Atmospheric Correction Implementation

[Matthew McCarthy]

Hello Robert and all,

I would like to implement Py6S for automated atmospheric correction of high-resolution satellite imagery, and want to clarify a few things before scaling up.

1. When applying the xa, xb, xc coefficients to convert from TOA radiance to surface reflectance, are these accounting for Earth-Sun distance, irradiance, and cos(theta)? Theta being the solar zenith angle.

2. Because I want to generate these coefficients and correct every pixel in an image, should I use 0 (zero) for my radiance input to AtmosCorrLambertianFromRadiance(radiance)? I am creating the coefficients for each of the multispectral channels individually.

3. Is there any validation I can cite to indicate that the coefficient method is similarly robust as the LUT approach with Py6S?

4. I noticed that changing the aerosol model (e.g. AeroProfile.Continental vs Maritime) and Lambertian vs BRDF functions does not change the coefficients - is that correct?

5. Is there a way to account for off-nadir viewing angles in the atmospheric correction scheme?

Thanks,

Matt McCarthy

Remote Sensing Group

Oak Ridge National Lab

[Robin Wilson]

Hi Matt,

Responses below:

On 24 September 2020 at 15:59:55, Matthew McCarthy (mj...@mail.usf.edu) wrote:

Hello Robert and all,

I would like to implement Py6S for automated atmospheric correction of high-resolution satellite imagery, and want to clarify a few things before scaling up.

1. When applying the xa, xb, xc coefficients to convert from TOA radiance to surface reflectance, are these accounting for Earth-Sun distance, irradiance, and cos(theta)? Theta being the solar zenith angle.

Yes they should take into account all of those.

2. Because I want to generate these coefficients and correct every pixel in an image, should I use 0 (zero) for my radiance input to AtmosCorrLambertianFromRadiance(radiance)? I am creating the coefficients for each of the multispectral channels individually.

It won't matter what you put in as the radiance input - the xa, xb and xc coefficients will be the same regardless what you enter, as they aren't based on the input radiance at all.

3. Is there any validation I can cite to indicate that the coefficient method is similarly robust as the LUT approach with Py6S?

I'm not aware of any good references on this. It will be more accurate to do a full LUT-based approach, as that can take into account the variation in various parameters across the image - for example, view angle, AOT changes and so on.

4. I noticed that changing the aerosol model (e.g. AeroProfile.Continental vs Maritime) and Lambertian vs BRDF functions does not change the coefficients - is that correct?

A change to the aerosol model should change the coefficients - I've tested this on my machine and I get different numbers when switching from Maritime to Continental. Can you send the code you're using, in case there are any errors in how you're setting the aerosol profile?

The coefficients don't take into account any angular settings for BRDF (eg. angular variation in reflectance), so there will be no difference using Lambertian vs BRDF functions.

5. Is there a way to account for off-nadir viewing angles in the atmospheric correction scheme?

Only by creating a lookup-table for viewing angles. The simple way to do this would be to create a lookup table that 'looked-up' the values of xa, xb and xc, so that you basically had a list of 10 viewing angles and then for each of them you had an xa, xb and xc value. Then, to look up the coefficients for a given view angle you interpolate within that table to get an estimate of the xa, xb and xc values (or, just find the nearest one, if you're not so concerned about absolutely accuracy). This would be fairly easy to implement using a for loop and some arrays/lists.

Hope that helps,

Robin

Thanks,

Matt McCarthy

Remote Sensing Group

Oak Ridge National Lab

--
You received this message because you are subscribed to the Google Groups "Py6S" group.
To unsubscribe from this group and stop receiving emails from it, send an email to py6s+uns...@googlegroups.com.
To view this discussion on the web, visit https://groups.google.com/d/msgid/py6s/1e5b8459-656e-4839-85b7-dd7829bffbb9n%40googlegroups.com.

[Matthew McCarthy]

Hi Robin,

Thanks for the quick reply and explanations. Please find below my code. I've run it with different predefined aerosol profiles and find the same coefficients.

    s = SixS()

    # Derive scene elevation from DEM

    # USGS Elevation Point Query Service

    url = r'https://nationalmap.gov/epqs/pqs.php?'

    df = pd.DataFrame({

        'lat': [lat],

        'lon': [lon]

    })

    elevations = []

    for lat, lon in zip(df['lat'], df['lon']):

        params = {

            'output': 'json',

            'x': lon,

            'y': lat,

            'units': 'Meters'

        }

        # format query string and return query value

        result = requests.get((url + urllib.parse.urlencode(params)))

        elevations.append(result.json()['USGS_Elevation_Point_Query_Service']['Elevation_Query']['Elevation'])

        elevation = elevations[0]/1000 # Convert meters to kilometers for Py6S

    s.altitudes.set_sensor_satellite_level()

    s.altitudes.set_target_custom_altitude(elevation)

    # Derive atmospheric profile from metadata (Water Vapor estimate)

    datetime = str(day) + '/' + str(month) + '/' + str(year)

    latitude = round(lat, -1)

    s.atmos_profile = AtmosProfile.FromLatitudeAndDate(latitude,datetime)

    # Aerosol model

#    s.aeroprofile = AeroProfile.PredefinedType(AeroProfile.Continental)

    s.aeroprofile = AeroProfile.PredefinedType(AeroProfile.Maritime)

    # Atmospheric Correction model

    s.atmos_corr = AtmosCorr.AtmosCorrLambertianFromRadiance(0)

    # Geometry

    solar_z = 90-meansunEL

    view_z = 90-meansatEL

    longitude = round(lon, -1)

    s.geometry = Geometry.User()

    s.geometry.solar_z = solar_z

    s.geometry.solar_a = meansunAZ

    s.geometry.view_z = view_z

    s.geometry.view_a = meansatAZ

    s.geometry.day = day

    s.geometry.month = month

    # Wavelengths

    a = np.zeros(8)

    b = np.zeros(8)

    c = np.zeros(8)

    for i in range(8):

        s.wavelength = Wavelength(wv_waves[i])

        s.run()

        a[i] = s.outputs.coef_xa

        b[i] = s.outputs.coef_xb

        c[i] = s.outputs.coef_xc

    print(a,b,c)

The coefficients for this image and each of 8 wavelengths are returning as:

[0.00426 0.00301 0.00317 0.00337 0.00354 0.00395 0.00608 0.00884] [0.25073 0.16275 0.10519 0.07921 0.06583 0.05668 0.04808 0.04565] [0.25481 0.20967 0.17038 0.14924 0.13649 0.12636 0.11373 0.10842]   

[Robin Wilson]

Hi,

I think there's a typo in your code - it should be s.aero_profile rather than s.aeroprofile.

Let me know if that solves it,

Robin

To view this discussion on the web, visit https://groups.google.com/d/msgid/py6s/c52048cc-d97e-4647-906e-41074b29d9f3n%40googlegroups.com.

[Matthew McCarthy]

Hi Robin,

Thanks, that did it! FYI the documentation (https://py6s.readthedocs.io/en/latest/params.html#aerosol-profiles) uses the "s.aeroprofile" syntax rather than the correct syntax you indicated.

Best,

Matt

[Robin Wilson]

Wow, those errors must have been there for years and no-one else noticed (or at least, no-one ever told me!).

I've fixed them now, the docs are just rebuilding.

Thanks for letting me know!

Robin

To view this discussion on the web, visit https://groups.google.com/d/msgid/py6s/675ecfae-fd6c-44d7-af08-2b485f92b771n%40googlegroups.com.

[Matthew McCarthy]

My pleasure. Thanks again!

Matt