1

I am trying to process a Sentinel 1 image to add layover and shadow masking. However when I do that and try and export the image to google drive I get an error: Error: Image.select: Pattern 'VH' did not match any bands.

I am using Colab python to do this and my code for exporting is here:

image = ee.Image('COPERNICUS/S1_GRD/S1A_IW_GRDH_1SSV_20160126T015925_20160126T015950_009659_00E14D_108C')
roi = {'long' : -120.318572, 'lat' : 39.310503000000004}
coords = ee.Geometry.Point(roi['long'], roi['lat'])
clip_area = coords.buffer(2000)
image_post = preprocess(image)
image_clipped = image.clip(clip_area)
image_clipped_post = preprocess(image_clipped)
image_post_clipped = image_post.clip(clip_area)
ID = 'testing_processing_v2'
folder = 'test'
imgGeom = image_clipped_post.geometry()
geometry = ee.Geometry.Rectangle([roi['long'] - 0.05, roi['lat'] - 0.05, roi['long'] + 0.05, roi['lat'] + 0.05])
export_image(image_post_clipped.select('VV'), ID, folder, imgGeom)

The function I am using to preprocess the image is:

 # Implementation by Andreas Vollrath (ESA), inspired by Johannes Reiche (Wageningen)
def preprocess(image):
  return mask_border(
      mask_overlay(
        terrain_correction(image)
      )  
  )

# Implementation by Andreas Vollrath (ESA), inspired by Johannes Reiche (Wageningen)
def mask_border(image):
    totalSlices = ee.Number(image.get('totalSlices'))
    sliceNumber = ee.Number(image.get('sliceNumber'))
    middleSlice = ee.Image(sliceNumber.gt(1).And(sliceNumber.lt(totalSlices)))
    mask = image.select(['VV', 'VH']).mask().reduce(ee.Reducer.min()).floor()
    pixelsToMask = mask.Not()\
      .fastDistanceTransform(128, 'pixels').sqrt()
    metersToMask = pixelsToMask\
      .multiply(ee.Image.pixelArea().sqrt())\
      .rename('metersToMask')
    notBorder = metersToMask.gte(500).And(pixelsToMask.gt(2))
    angle = image.select('angle')
    return image\
        .updateMask(
            angle.gt(31).And(angle.lt(45))
            .And(middleSlice.Or(notBorder))
        )

# Implementation by Andreas Vollrath (ESA), inspired by Johannes Reiche (Wageningen)
def mask_overlay(image):
  return image.updateMask(
      image.select('shadow').And(image.select('layover'))
  )

 # Implementation by Andreas Vollrath (ESA), inspired by Johannes Reiche (Wageningen)
def terrain_correction(image):
  imgGeom = image.geometry()
  srtm = ee.Image('JAXA/ALOS/AW3D30_V1_1').clip(imgGeom) # 30m srtm 
  sigma0Pow = ee.Image.constant(10).pow(image.divide(10.0))

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

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

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

  # convert all to radians
  phi_rRad = phi_r.multiply(math.pi/180)
  alpha_sRad = alpha_s.multiply(math.pi/180)
  theta_iRad = theta_i.multiply(math.pi/180)
  ninetyRad = ee.Image.constant(90).multiply(math.pi/180)

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

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

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

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

  # apply model
  gamma0_Volume = gamma0.divide(volModel)
  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 
  alpha_rDeg = alpha_r.multiply(180/math.pi)
  layover = alpha_rDeg.lt(theta_i);

  # shadow where LIA > 90
  shadow = theta_liaDeg.lt(85);

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

  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']), 
    None, 
    False
  )

The things I have tried so far:

  1. check the bands on image are only 'VV' and 'angle'

  2. remove all references to 'VH' in preprocessing -> I get error - Error: Image.constant: Parameter 'value' is required. I would prefer not to remove the 'VH' as well since I may want to use this function with images that have 'VH' as well later.

  3. change order of clipping and processing -> didn't help

  4. When I remove the preprocessing step then this exports fine so I think it is in that code I am getting this error

I can include the code with all 'VH's removed but I was still getting that error.

I don't know if it matters but here is the code to the export function as well.

def export_image(image, ID, folder, geometry, scale = 10):
  image = image.toFloat()
  task = ee.batch.Export.image.toDrive(**{
  'image': image,
  'description': ID,
  'folder': folder,
  'scale': scale,
  'region': geometry
  })
  task.start()
  while task.active():
    print('Polling for task (id: {}).'.format(task.id))
    print('Status: ',task.status().get('state'))
    time.sleep(3)
  else:
    if task.status().get('state') == 'FAILED':
      print('FAILED')
    print("Done: " + ID)
  return()

This might be too much to get help on I have just been stuck for quite a few days and figured I would at least try here as well.

2

There are a couple of issues with the code:

  • You're invoking export_image() with the wrong argument order.
  • Like you pointed out, your image doesn't contain a VH band. You have to remove all references to it. Instead of doing that, I typically filter out images that doesn't contain both VV and VH.
  • There is problems with your geometry. I didn't try to figure out what, because you've clipped your image, so you actually don't need to provide the geometry at all.

This exports for me:

import ee
import math
import time
ee.Initialize()

def preprocess(image):
  return mask_border(
      mask_overlay(
        terrain_correction(image)
      )  
  )

def mask_border(image):
    totalSlices = ee.Number(image.get('totalSlices'))
    sliceNumber = ee.Number(image.get('sliceNumber'))
    middleSlice = ee.Image(sliceNumber.gt(1).And(sliceNumber.lt(totalSlices)))
    mask = image.select(['VV']).mask().reduce(ee.Reducer.min()).floor()
    pixelsToMask = mask.Not()\
      .fastDistanceTransform(128, 'pixels').sqrt()
    metersToMask = pixelsToMask\
      .multiply(ee.Image.pixelArea().sqrt())\
      .rename('metersToMask')
    notBorder = metersToMask.gte(500).And(pixelsToMask.gt(2))
    angle = image.select('angle')
    return image\
        .updateMask(
            angle.gt(31).And(angle.lt(45))
            .And(middleSlice.Or(notBorder))
        )

def mask_overlay(image):
  return image.updateMask(
      image.select('shadow').And(image.select('layover'))
  )

 # Implementation by Andreas Vollrath (ESA), inspired by Johannes Reiche (Wageningen)
def terrain_correction(image):
  imgGeom = image.geometry()
  srtm = ee.Image('JAXA/ALOS/AW3D30_V1_1').clip(imgGeom) # 30m srtm 
  sigma0Pow = ee.Image.constant(10).pow(image.divide(10.0))

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

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

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

  # convert all to radians
  phi_rRad = phi_r.multiply(math.pi/180)
  alpha_sRad = alpha_s.multiply(math.pi/180)
  theta_iRad = theta_i.multiply(math.pi/180)
  ninetyRad = ee.Image.constant(90).multiply(math.pi/180)

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

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

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

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

  # apply model
  gamma0_Volume = gamma0.divide(volModel)
  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 
  alpha_rDeg = alpha_r.multiply(180/math.pi)
  layover = alpha_rDeg.lt(theta_i);

  # shadow where LIA > 90
  shadow = theta_liaDeg.lt(85);

  # calculate the ratio for RGB vis
  ratio = gamma0_VolumeDB.select('VV')

  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', 'slope_1', 'slope_2'], ['VV', 'layover', 'shadow']), 
    None, 
    False
  )

def export_image(image, ID, folder, scale = 10):
  image = image.toFloat()
  task = ee.batch.Export.image.toDrive(**{
  'image': image,
  'description': ID,
  'folder': folder,
  'scale': scale
#   'region': geometry
  })
  task.start()
  while task.active():
    print('Polling for task (id: {}).'.format(task.id))
    print('Status: ',task.status().get('state'))
    time.sleep(3)
  else:
    if task.status().get('state') == 'FAILED':
      print('FAILED')
    print("Done: " + ID)
  return()


image = ee.Image('COPERNICUS/S1_GRD/S1A_IW_GRDH_1SSV_20160126T015925_20160126T015950_009659_00E14D_108C')
roi = {'long' : -120.318572, 'lat' : 39.310503000000004}
coords = ee.Geometry.Point(roi['long'], roi['lat'])
clip_area = coords.buffer(2000)
image_post = preprocess(image)
image_clipped = image.clip(clip_area)
image_clipped_post = preprocess(image_clipped)
image_post_clipped = image_post.clip(clip_area)
ID = 'testing_processing_v2'
folder = 'test'
geometry = ee.Geometry.Rectangle([roi['long'] - 0.05, roi['lat'] - 0.05, roi['long'] + 0.05, roi['lat'] + 0.05])
export_image(image=image_post_clipped.select('VV'), ID=ID, folder=folder)
  • Thanks so much! I didn't realize that I didn't need to provide a geometry with clipped image. I have relatively few images to work with so I think I will just create an alternative set of functions to preprocess with ['VH'] included as well. Thanks again! – ZachK Jan 21 at 20:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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