2

I am using python API in google earth engine (colab) to retrieve local incidence angle from sentinel-1 GRD data. However, the .select ('') function only allows the VV, VH, HH, HV, and 'angle' band selection. The angle band is viewing angle of the sensor, and I need a local incidence angle (retrieved as a by-product of the terrain correction step). Does anyone have the answer? Thanks

2 Answers 2

4

The terrain correction implemented by Andreas Vollrath calculates the local incidence angle. You can modify the implementation to return that angle.

// Implementation by Andreas Vollrath (ESA), inspired by Johannes Reiche (Wageningen)
function terrainCorrection(image) {
  var imgGeom = image.geometry()
  var srtm = ee.Image('USGS/SRTMGL1_003').clip(imgGeom) // 30m srtm 
  var sigma0Pow = ee.Image.constant(10).pow(image.divide(10.0))

  // Article ( numbers relate to chapters)
  // 2.1.1 Radar geometry 
  var theta_i = image.select('angle')
  var phi_i = ee.Terrain.aspect(theta_i)
    .reduceRegion(ee.Reducer.mean(), theta_i.get('system:footprint'), 1000)
    .get('aspect')

  // 2.1.2 Terrain geometry
  var alpha_s = ee.Terrain.slope(srtm).select('slope')
  var phi_s = ee.Terrain.aspect(srtm).select('aspect')

  // 2.1.3 Model geometry
  // reduce to 3 angle
  var phi_r = ee.Image.constant(phi_i).subtract(phi_s)

  // convert all to radians
  var phi_rRad = phi_r.multiply(Math.PI / 180)
  var alpha_sRad = alpha_s.multiply(Math.PI / 180)
  var theta_iRad = theta_i.multiply(Math.PI / 180)
  var ninetyRad = ee.Image.constant(90).multiply(Math.PI / 180)

  // slope steepness in range (eq. 2)
  var alpha_r = (alpha_sRad.tan().multiply(phi_rRad.cos())).atan()

  // slope steepness in azimuth (eq 3)
  var alpha_az = (alpha_sRad.tan().multiply(phi_rRad.sin())).atan()

  // local incidence angle (eq. 4)
  var theta_lia = (alpha_az.cos().multiply((theta_iRad.subtract(alpha_r)).cos())).acos()
  var theta_liaDeg = theta_lia.multiply(180 / Math.PI)
  // 2.2 
  // Gamma_nought_flat
  var gamma0 = sigma0Pow.divide(theta_iRad.cos())
  var gamma0dB = ee.Image.constant(10).multiply(gamma0.log10())
  var ratio_1 = gamma0dB.select('VV').subtract(gamma0dB.select('VH'))

  // Volumetric Model
  var nominator = (ninetyRad.subtract(theta_iRad).add(alpha_r)).tan()
  var denominator = (ninetyRad.subtract(theta_iRad)).tan()
  var volModel = (nominator.divide(denominator)).abs()

  // apply model
  var gamma0_Volume = gamma0.divide(volModel)
  var gamma0_VolumeDB = ee.Image.constant(10).multiply(gamma0_Volume.log10())

  // we add a layover/shadow maskto the original implmentation
  // layover, where slope > radar viewing angle 
  var alpha_rDeg = alpha_r.multiply(180 / Math.PI)
  var layover = alpha_rDeg.lt(theta_i);

  // shadow where LIA > 90
  var shadow = theta_liaDeg.lt(85)

  // calculate the ratio for RGB vis
  var ratio = gamma0_VolumeDB.select('VV').subtract(gamma0_VolumeDB.select('VH'))

  var output = gamma0_VolumeDB.addBands(ratio).addBands(alpha_r).addBands(phi_s).addBands(theta_iRad)
    .addBands(layover).addBands(shadow).addBands(gamma0dB).addBands(ratio_1)

  return image.addBands(
    output.select(['VV', 'VH', 'slope_1', 'slope_2'], ['VV', 'VH', 'layover', 'shadow']),
    null,
    true
  )
}
1
  • I tried to modify the code to get the local incidence angle but it is showing values of incidence angle between 0 and 1 (int). The code is provided below Mar 4, 2020 at 10:28
0

Code:

var s1 = ee.ImageCollection('COPERNICUS/S1_GRD')
    .filterDate('2019-01-01', '2019-01-31')
    .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
    .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH'))
    .filter(ee.Filter.eq('instrumentMode', 'IW'))
    .filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING'))
    .filter(ee.Filter.eq('resolution', 'H'))
    .filter(ee.Filter.eq('resolution_meters', 10))
    .filterBounds(geometry)
    .select('VV','VH','angle');
print('Collection: S1', s1);

var image = ee.Image(
  s1.first())
print('Test image', image);

// Implementation by Andreas Vollrath (ESA), inspired by Johannes Reiche (Wageningen)
 var terrainCorrection= function(image) {
  var imgGeom = image.geometry()
  var srtm = ee.Image('CGIAR/SRTM90_V4').clip(imgGeom) // 30m srtm 
  var sigma0Pow = ee.Image.constant(10).pow(image.divide(10.0))

  // Article ( numbers relate to chapters)
  // 2.1.1 Radar geometry 
  var theta_i = image.select('angle')
  var phi_i = ee.Terrain.aspect(theta_i)
    .reduceRegion(ee.Reducer.mean(), theta_i.get('system:footprint'), 1000)
    .get('aspect')

  // 2.1.2 Terrain geometry
  var alpha_s = ee.Terrain.slope(srtm).select('slope')
  var phi_s = ee.Terrain.aspect(srtm).select('aspect')

  // 2.1.3 Model geometry
  // reduce to 3 angle
  var phi_r = ee.Image.constant(phi_i).subtract(phi_s)

  // convert all to radians
  var phi_rRad = phi_r.multiply(Math.PI / 180)
  var alpha_sRad = alpha_s.multiply(Math.PI / 180)
  var theta_iRad = theta_i.multiply(Math.PI / 180)
  var ninetyRad = ee.Image.constant(90).multiply(Math.PI / 180)

  // slope steepness in range (eq. 2)
  var alpha_r = (alpha_sRad.tan().multiply(phi_rRad.cos())).atan()

  // slope steepness in azimuth (eq 3)
  var alpha_az = (alpha_sRad.tan().multiply(phi_rRad.sin())).atan()

  // local incidence angle (eq. 4)
  var theta_lia = (alpha_az.cos().multiply((theta_iRad.subtract(alpha_r)).cos())).acos()
  var theta_liaDeg = theta_lia.multiply(180 / Math.PI)
  var local=theta_lia.toFloat();
  // 2.2 
  // Gamma_nought_flat
  var gamma0 = sigma0Pow.divide(theta_iRad.cos())
  var gamma0dB = ee.Image.constant(10).multiply(gamma0.log10())
  var ratio_1 = gamma0dB.select('VV').subtract(gamma0dB.select('VH'))

  // Volumetric Model
  var nominator = (ninetyRad.subtract(theta_iRad).add(alpha_r)).tan()
  var denominator = (ninetyRad.subtract(theta_iRad)).tan()
  var volModel = (nominator.divide(denominator)).abs()

  // apply model
  var gamma0_Volume = gamma0.divide(volModel)
  var gamma0_VolumeDB = ee.Image.constant(10).multiply(gamma0_Volume.log10())

  // we add a layover/shadow maskto the original implmentation
  // layover, where slope > radar viewing angle 
  var alpha_rDeg = alpha_r.multiply(180 / Math.PI)
  var layover = alpha_rDeg.lt(theta_i);

  // shadow where LIA > 90
  var shadow = theta_liaDeg.lt(85)

  // calculate the ratio for RGB vis
  var ratio = gamma0_VolumeDB.select('VV').subtract(gamma0_VolumeDB.select('VH'))

  var output = gamma0_VolumeDB.addBands(ratio).addBands(alpha_r).addBands(phi_s).addBands(theta_iRad)
    .addBands(layover).addBands(shadow).addBands(gamma0dB).addBands(ratio_1).addBands(srtm).addBands(local)

  return image.addBands(
    output.select(['VV', 'VH', 'slope_1', 'slope_2'], ['VV', 'VH', 'layover', 'local']),
    null,
    true
  )
}

var test_output= terrainCorrection(image)
print("test_output",test_output)


var s1_all= s1.map(terrainCorrection)
print('s1_all', s1_all);
Map.addLayer(
  s1_all.mean(),
  {bands:'local', min:-20, max:2},
  'Incidence angle'
);

Output:

List (5 elements)
0:"VV", double, EPSG:32631, 28746x21405 px
1:"VH", double, EPSG:32631, 28746x21405 px
2:"angle", float, EPSG:32631, 21x10 px
3:"layover", int ∈ [0, 1], EPSG:4326, 5352x2670 px
4:"local", int ∈ [0, 1], EPSG:4326, 5352x2670 px
1
  • It is not showing the required range of local incidence angle, also srtm DEM is giving values between 0 and 1 Mar 4, 2020 at 10:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.