Angular-based Radiometric Slope Correction of Sentinel-1

Angular-based Radiometric Slope Correction of Sentinel-1

📆 Project Period	December, 2023
📍 GitHub	gitlab.esa.int

Angular-based radiometric slope correction of Sentinel-1

Radiometric slope correction of Sentinel-1 data on Google Earth Engine.

This notebook contains the actual slope correction routine and is used to generate and export the resultant layers to Google Drive for further plotting. It also provides the supplementary material and code implemented in the related reference paper : Angular-Based Radiometric Slope Correction for Sentinel-1 on Google Earth Engine

Note:

This step will take up to a couple of hours and is fundametal to be finished before running the subsequent notebooks of the article.
You will need a valid Google account.
You will need to be registered for Earth Engine.
You will need about 6 GB of free space on your Google Drive

1 - Get connected to Google Earth Engine

import ee
# Trigger the authentication flow.
ee.Authenticate()

2 - Import necessary python libs

import numpy as np
ee.Initialize()

3 - The Model function

This is a nested function, so that we can actually map over a collection, by simultaneously giving the DEM, model of choice and buffer parameter.

4 - Create corrected data for volume and surface model on Earth Engine and export both to Google Drive

Within this code, we apply the slope correction function to two consecutive frames of Sentinel-1 using the volume and surface model seaprately. All layers are exported to Google Drive. The slooe correction function does not only export the corrected layers, but also all the angles created during routine for later plotting.

5 - Create buffered layover and shadow masks

This cell creates the buffered layover and shadow masks for 30, 50 and 100 meters used in Figure 6 of the article.

def slope_correction(collection, elevation, model, buffer=0):
    '''This function applies the slope correction on a collection of Sentinel-1 data

       :param collection: ee.Collection of Sentinel-1
       :param elevation: ee.Image of DEM
       :param model: model to be applied (volume/surface)
       :param buffer: buffer in meters for layover/shadow amsk

        :returns: ee.Image
    '''

    def _volumetric_model_SCF(theta_iRad, alpha_rRad):
        '''Code for calculation of volumetric model SCF

        :param theta_iRad: ee.Image of incidence angle in radians
        :param alpha_rRad: ee.Image of slope steepness in range

        :returns: ee.Image
        '''

        # create a 90 degree image in radians
        ninetyRad = ee.Image.constant(90).multiply(np.pi/180)

        # model
        nominator = (ninetyRad.subtract(theta_iRad).add(alpha_rRad)).tan()
        denominator = (ninetyRad.subtract(theta_iRad)).tan()
        return nominator.divide(denominator)


    def _surface_model_SCF(theta_iRad, alpha_rRad, alpha_azRad):
        '''Code for calculation of direct model SCF

        :param theta_iRad: ee.Image of incidence angle in radians
        :param alpha_rRad: ee.Image of slope steepness in range
        :param alpha_azRad: ee.Image of slope steepness in azimuth

        :returns: ee.Image
        '''

        # create a 90 degree image in radians
        ninetyRad = ee.Image.constant(90).multiply(np.pi/180)

        # model
        nominator = (ninetyRad.subtract(theta_iRad)).cos()
        denominator = (alpha_azRad.cos()
          .multiply((ninetyRad.subtract(theta_iRad).add(alpha_rRad)).cos()))

        return nominator.divide(denominator)


    def _erode(image, distance):
      '''Buffer function for raster

      :param image: ee.Image that shoudl be buffered
      :param distance: distance of buffer in meters

      :returns: ee.Image
      '''

      d = (image.Not().unmask(1)
          .fastDistanceTransform(30).sqrt()
          .multiply(ee.Image.pixelArea().sqrt()))

      return image.updateMask(d.gt(distance))


    def _masking(alpha_rRad, theta_iRad, buffer):
        '''Masking of layover and shadow


        :param alpha_rRad: ee.Image of slope steepness in range
        :param theta_iRad: ee.Image of incidence angle in radians
        :param buffer: buffer in meters

        :returns: ee.Image
        '''
        # layover, where slope > radar viewing angle
        layover = alpha_rRad.lt(theta_iRad).rename('layover')

        # shadow
        ninetyRad = ee.Image.constant(90).multiply(np.pi/180)
        shadow = alpha_rRad.gt(ee.Image.constant(-1).multiply(ninetyRad.subtract(theta_iRad))).rename('shadow')

        # add buffer to layover and shadow
        if buffer > 0:
            layover = _erode(layover, buffer)
            shadow = _erode(shadow, buffer)

        # combine layover and shadow
        no_data_mask = layover.And(shadow).rename('no_data_mask')

        return layover.addBands(shadow).addBands(no_data_mask)


    def _correct(image):
        '''This function applies the slope correction and adds layover and shadow masks

        '''

        # get the image geometry and projection
        geom = image.geometry()
        proj = image.select(1).projection()

        # calculate the look direction
        heading = (ee.Terrain.aspect(image.select('angle'))
                                     .reduceRegion(ee.Reducer.mean(), geom, 1000)
                                     .get('aspect'))


        # Sigma0 to Power of input image
        sigma0Pow = ee.Image.constant(10).pow(image.divide(10.0))

        # the numbering follows the article chapters
        # 2.1.1 Radar geometry
        theta_iRad = image.select('angle').multiply(np.pi/180)
        phi_iRad = ee.Image.constant(heading).multiply(np.pi/180)

        # 2.1.2 Terrain geometry
        alpha_sRad = ee.Terrain.slope(elevation).select('slope').multiply(np.pi/180).setDefaultProjection(proj).clip(geom)
        phi_sRad = ee.Terrain.aspect(elevation).select('aspect').multiply(np.pi/180).setDefaultProjection(proj).clip(geom)

        # we get the height, for export
        height = elevation.setDefaultProjection(proj).clip(geom)

        # 2.1.3 Model geometry
        #reduce to 3 angle
        phi_rRad = phi_iRad.subtract(phi_sRad)

        # slope steepness in range (eq. 2)
        alpha_rRad = (alpha_sRad.tan().multiply(phi_rRad.cos())).atan()

        # slope steepness in azimuth (eq 3)
        alpha_azRad = (alpha_sRad.tan().multiply(phi_rRad.sin())).atan()

        # local incidence angle (eq. 4)
        theta_liaRad = (alpha_azRad.cos().multiply((theta_iRad.subtract(alpha_rRad)).cos())).acos()
        theta_liaDeg = theta_liaRad.multiply(180/np.pi)

        # 2.2
        # Gamma_nought
        gamma0 = sigma0Pow.divide(theta_iRad.cos())
        gamma0dB = ee.Image.constant(10).multiply(gamma0.log10()).select(['VV', 'VH'], ['VV_gamma0', 'VH_gamma0'])
        ratio_gamma = (gamma0dB.select('VV_gamma0')
                        .subtract(gamma0dB.select('VH_gamma0'))
                        .rename('ratio_gamma0'))

        if model == 'volume':
            scf = _volumetric_model_SCF(theta_iRad, alpha_rRad)

        if model == 'surface':
            scf = _surface_model_SCF(theta_iRad, alpha_rRad, alpha_azRad)

        # apply model for Gamm0_f
        gamma0_flat = gamma0.divide(scf)
        gamma0_flatDB = (ee.Image.constant(10)
                         .multiply(gamma0_flat.log10())
                         .select(['VV', 'VH'],['VV_gamma0flat', 'VH_gamma0flat'])
                        )

        masks = _masking(alpha_rRad, theta_iRad, buffer)

        # calculate the ratio for RGB vis
        ratio_flat = (gamma0_flatDB.select('VV_gamma0flat')
                        .subtract(gamma0_flatDB.select('VH_gamma0flat'))
                        .rename('ratio_gamma0flat')
                     )

        return (image.rename(['VV_sigma0', 'VH_sigma0', 'incAngle'])
                      .addBands(gamma0dB)
                      .addBands(ratio_gamma)
                      .addBands(gamma0_flatDB)
                      .addBands(ratio_flat)
                      .addBands(alpha_rRad.rename('alpha_rRad'))
                      .addBands(alpha_azRad.rename('alpha_azRad'))
                      .addBands(phi_sRad.rename('aspect'))
                      .addBands(alpha_sRad.rename('slope'))
                      .addBands(theta_iRad.rename('theta_iRad'))
                      .addBands(theta_liaRad.rename('theta_liaRad'))
                      .addBands(masks)
                      .addBands(height.rename('elevation'))
                 )

    # run and return correction
    return collection.map(_correct)

# geometry for AOI in Autria
geometry = ee.Geometry.Polygon([[[11.25, 47.35],
                    [11.25, 47.05],
                    [11.7, 47.05],
                    [11.7, 47.35]]])

# land cover of Austria
lc = ee.Image('users/andreasvollrath/DEM_Austria/FP_P1_HR_LandCoverMap2016_L2_clipped_wgs84_nn').rename('landcover')

# filter Sentinel-1 collection for AOI and selected dates
s1Collection = ee.ImageCollection('COPERNICUS/S1_GRD') \\
        .filterBounds(geometry) \\
        .filterDate('2016-08-15', '2016-08-16')


# paths to dem
dem = 'USGS/SRTMGL1_003'

# list of models
models = ['volume', 'surface']

# this is the scale we want the data to be sampled
scale = 10

# loop through all combinations and export to drive
for model in models:

    # get the respective collection and bands and mosaic to a single image
    corrected_image = slope_correction(
        s1Collection,
        ee.Image(dem),
        model
    ).mosaic()

    # we get geometry and projection of the image
    proj = corrected_image.select(1).projection()
    geom = corrected_image.clip(geometry).select(1).geometry()

    # add lc and bring everything to same projection/geometry
    added_LC = corrected_image.addBands(lc)
    image_reprojected = added_LC.reproject(proj, scale=scale).clip(geom)

    # get the bandlist
    bandlist = image_reprojected.bandNames().getInfo()

    # create an export job for each band
    for band in bandlist:
        task = ee.batch.Export.image.toDrive(
            image=image_reprojected.select(band).clip(geom),
            description=band,
            folder='slope_correction/{}_{}_buf_0'.format(dem.split('/')[-1], model),
            fileNamePrefix=band,
            region=geom.coordinates().getInfo(),
            scale=scale,
            maxPixels=1e12
            )
        task.start()
        print(task.status())

# list of buffers
buffers = [30, 50, 100]
model = 'volume'

for buffer in buffers:

    # get the respective collection and bands and mosaic to a single image
    corrected_image = slope_correction(
        s1Collection,
        ee.Image(dem),
        model,
        buffer=buffer
    ).mosaic()

    # create an export job for each band
    for band in ['layover', 'shadow']:
        task = ee.batch.Export.image.toDrive(
            image=corrected_image.select(band).clip(geom),
            description=band,
            folder='slope_correction/{}_volume_buf_{}'.format(dem.split('/')[-1], buffer),
            fileNamePrefix=band,
            region=geom.coordinates().getInfo(),
            scale=scale,
            maxPixels=1e12
            )
        task.start()
        print(task.status())