ImgMask.py

import ee

import Image as Img
import eoAuxData as eoAD



CLEAR_MASK = 1
WATER_MASK = 2
SNOW_MASK  = 3
SATU_MASK  = 4   #Radiometric satuation




#############################################################################################################
# Description: This function creates a clear-sky pixel mask for a given image (Image) based on the comparison
#              with a MODIS mosaic image (MODIS_mosaic).
#
# Revision history:  2023-Apr-14  Lixin Sun  Initial creation
#                    2023-May-25  Lixin Sun  Treat high spatial resolution vegetated and non-vegetated 
#                                            pixels differently.
#############################################################################################################
def mask_from_MODIS(Image, SsrData, MODIS_mosaic): 
  '''This function creates a clear-sky pixel mask for a given image (Image) based on the comparison with a
     MODIS mosaic image (MODIS_mosaic)

     Args:
       Image(ee.Image): a given Sentinel-2 or Landsat image;
       SsrData(Dictionary): an sensor data dictionary;
       MODIS_mosaic(ee.Image): an compariable MODIS mosaic image.'''    
  #==========================================================================================================
  # Apply gain and offset to both MODIS and S2/LS images
  #==========================================================================================================
  modis_ssrData = Img.SSR_META_DICT['MOD_SR']
  modis = Img.apply_gain_offset(MODIS_mosaic, modis_ssrData, 100.0, False)
  image = Img.apply_gain_offset(Image, SsrData, 100.0, False)

  #==========================================================================================================
  # Obtain related bands in high-resolution image and then calculate NDVI map
  #==========================================================================================================
  s2ls_blu  = image.select(SsrData['BLU'])
  s2ls_red  = image.select(SsrData['BLU'])
  s2ls_nir  = image.select(SsrData['NIR'])
  ndvi      = s2ls_nir.subtract(s2ls_red).divide(s2ls_nir.add(s2ls_red))

  #==========================================================================================================
  # Create two adjustment maps for BLUE and NIR bands, respectively 
  # For a high-resolution vegetated pixel, its blu <= MODIS_blue
  # For a high-resolution non-vegetated pixel, its blu >= MODIS_blue
  # For a high-resolution vegetated pixel, its NIR >= MODIS_NIR
  # For a high-resolution non-vegetated pixel, its NIR <= MODIS_NIR
  #==========================================================================================================
  blu_adjust = s2ls_blu.multiply(0).add(2)
  blu_adjust = blu_adjust.where(ndvi.lt(0.4), blu_adjust.add(3)) 
  
  nir_adjust = s2ls_nir.multiply(0).add(5)
  nir_adjust = nir_adjust.where(ndvi.lt(0.4), blu_adjust.subtract(3))

  #==========================================================================================================
  # Apply blue and NIR adjustment images to MODIS blue and NIR bands, respectively
  #==========================================================================================================
  modis_blu = modis.select(modis_ssrData['BLU']).add(blu_adjust)
  modis_nir = modis.select(modis_ssrData['NIR']).subtract(nir_adjust)

  #==========================================================================================================
  # Create a pixel mask
  #==========================================================================================================
  mask = s2ls_blu.multiply(0)
  
  cond = s2ls_blu.gt(modis_blu).Or(s2ls_nir.lt(modis_nir))
  mask = mask.where(cond, ee.Image(1))

  return cond





#############################################################################################################
# Description: This function extracts a type of mask from the specific bands in a given image. The masks were
#              created by the image provider.
#
# Note:        the returned mask is a 0/1 image with 1 represening specified targets (e.g., cloud, cirrus,
#              shadow, water or snow/ice)  
#
# Revision history:  2022-Jun-22  Lixin Sun  Initial creation
#                    2023-Jan-11  Lixin Sun  Added MODIS sensor code option and MODIS mosaic image.
#  
#############################################################################################################
def Img_VenderMask(Image, SsrData, MaskType, MODIS_mosaic = None):
  '''This function extracts a specified mask from the intrinsic QA band of a given Landsat image.

     Args:
       Image(ee.Image): a given image with the bands of a specified sensor;
       SsrData(Dictionary): a dictionary containing some info on a sensor  (LS 5/7/8/9 and S2);       
       MaskType(int): the mask type code (CLEAR_MASK, WATER_MASK, SNOW_MASK and SATU_MASK);
       MODIS_mosaic(ee.Image): An optional MODIS mosaic image.'''         
 
  ssr_code  = SsrData['SSR_CODE']
  data_unit = SsrData['DATA_UNIT']
  mask_type = int(MaskType)  

  if ssr_code == Img.MOD_sensor:  # For MODIS image
    cloudShadow = ee.Image.constant(1 << 2)
    cirrus1     = ee.Image.constant(1 << 8)
    cirrus2     = ee.Image.constant(1 << 9)
    cloud       = ee.Image.constant(1 << 10)
    snow        = ee.Image.constant(1 << 12)
    cloud2      = ee.Image.constant(1 << 13)
    snow2       = ee.Image.constant(1 << 15)

    qa = Image.select('StateQA')
    # Both flags should be set to zero, indicating clear conditions.
    mask = qa.bitwiseAnd(cloudShadow).Or(qa.bitwiseAnd(cirrus1)).Or(qa.bitwiseAnd(cirrus2)) \
          .Or(qa.bitwiseAnd(cloud)).Or(qa.bitwiseAnd(snow)).Or(qa.bitwiseAnd(cloud2)).Or(qa.bitwiseAnd(snow2))
    
    return mask
  elif ssr_code > Img.MAX_LS_CODE:  # For Sentinel-2 image
    # For Sentinel-2, only two bands, 'QA60' and 'SCL', include mask information  
    qa  = Image.select(['QA60']).uint16()
    scl = Image.select(['SCL']) if data_unit == 2 else qa.multiply(0)
  
    if mask_type == CLEAR_MASK:
      cloud  = ee.Image.constant(1 << 10)   # Opaque clouds
      cirrus = ee.Image.constant(1 << 11)   # Cirrus clouds

      mask = qa.bitwiseAnd(cloud).Or(qa.bitwiseAnd(cirrus))
      mask = mask.Or(scl.eq(3)).Or(scl.eq(8)).Or(scl.eq(9)).Or(scl.eq(10))

      if MODIS_mosaic != None:        
        modis_mask = mask_from_MODIS(Image, SsrData, MODIS_mosaic)
        mask = mask.And(modis_mask)
      
      return mask
    elif mask_type == WATER_MASK:
      return scl.eq(6)
    elif mask_type == SNOW_MASK:
      return scl.eq(11)
    elif mask_type == SATU_MASK:
      return scl.eq(1)
    else:
      return qa.multiply(ee.Image(0))
  else:   # For both BOA and TOA reflectance data of Landsat 5/7/8/9
    # For Landsat, only one band, 'QA_PIXEL', includes mask information
    '''
      =================================================
       TOA reflectance of LS 5/7/8/9
      =================================================
      Bit 1: Dilated Cloud
        0: Cloud is not dilated or no cloud
        1: cloud dilation
      Bit 2: Unused
      Bit 3: Cloud
        0: Cloud confidence is not high
        1: High confidence cloud
      Bit 4: Cloud Shadow
        0: Cloud Shadow Confidence is not high
        1: High confidence cloud shadow

      =================================================
       BOA reflectance of LS 5/7/8/9
      =================================================
      Bit 1: Dilated Cloud
      Bit 2: Cirrus (high confidence)
      Bit 3: Cloud
      Bit 4: Cloud Shadow
    '''
    qa = Image.select(['QA_PIXEL']).uint16()
    dilated = ee.Image.constant(1 << 1)
    cirrus  = ee.Image.constant(1 << 2)
    cloud   = ee.Image.constant(1 << 3)
    shadow  = ee.Image.constant(1 << 4)
    snow    = ee.Image.constant(1 << 5)
    water   = ee.Image.constant(1 << 7)

    if mask_type == CLEAR_MASK:
      mask = qa.bitwiseAnd(dilated).Or(qa.bitwiseAnd(cirrus)).Or(qa.bitwiseAnd(cloud)).Or(qa.bitwiseAnd(shadow))
      
      if MODIS_mosaic != None:        
        modis_mask = mask_from_MODIS(Image, SsrData, MODIS_mosaic)
        mask = mask.And(modis_mask)

      return mask
    elif mask_type == WATER_MASK:
      return qa.bitwiseAnd(water)   # Bit 7: Water
    elif mask_type == SNOW_MASK:
      return qa.bitwiseAnd(snow)    # Bit 5: Snow
    #elif mask_type == SATU_MASK:
      #sa = Image.select(['QA_RADSAT']).uint8()
      #return sa.bitwiseAnd(127)
      #return sa.multiply(0)
    else:
      return qa.multiply(0)





#############################################################################################################
# Description: This function creates a value-invalid pixel mask (1 => invalid value pixels) for an image.
# 
# Revision history:  2022-Nov-11  Lixin Sun  Initial creation
#                    2022-Nov-15  Lixin Sun  Removed the limit that maximum reflectance value must be 100.
#
#############################################################################################################
def Img_ValueMask(Image, SsrData, MaxRef):
  '''Creates a value-invalid pixel mask (1 => value_invalid pixel) for an image.

     Args:
       Image(ee.Image): a given ee.Image object;
       SsrData(Dictionary): A Dictionary with metadata for a sensor and data unit;
       MaxRef(int): a maximum reflectance value (1 or 100).'''
  #===========================================================================================================
  # Extract optical bands 
  #===========================================================================================================
  band_names = SsrData['OUT_BANDS']
  used_img   = Image.select(band_names) 

  #===========================================================================================================
  # Create valid value mask
  #===========================================================================================================
  max_val = ee.Number(1.05).multiply(MaxRef)
  min_val = ee.Number(-0.005).multiply(MaxRef)
  
  mask = used_img.lt(min_val).Or(used_img.gt(max_val))

  return mask.reduce(ee.Reducer.max()).rename(['ValMask'])




#############################################################################################################
# Description: This function creates a vegetation pixel mask (1 => vegetation) for an image.
#
# Revision history:  2022-Nov-11  Lixin Sun  Initial creation
#                    2022-Nov-15  Lixin Sun  Removed the limit that maximum reflectance value must be 100.
# 
#############################################################################################################  
def Img_VegMask(Image, SsrData):
  '''Creates a vegetation pixel mask (1 => vegetation) for an image.

     Args:
       Image(ee.Image): a given ee.Image object;
       SsrData(Dictionary): A Dictionary containing metadata associated with a sensor and data unit.'''
 
  blu  = Image.select(SsrData['BLU'])
  grn  = Image.select(SsrData['GRN'])
  red  = Image.select(SsrData['RED'])

  ndvi = Image.normalizedDifference([SsrData['NIR'], SsrData['RED']])

  return ee.Image(ndvi.gt(ee.Image(0.3)).And(grn.gt(blu)).And(grn.gt(red)))
  





#############################################################################################################
# Description: This function creates a non-vegetated pixel mask (1 => non-vegetated pixels) for an image.
#
# Revision history:  2021-Jul-07  Lixin Sun  Initial creation.
#                    2022-Nov-15  Lixin Sun  Removed the limit that maximum reflectance value must be 100.
#
#############################################################################################################
def Img_NonVegMask(Image, SsrData, MaxRef, indx_name):
  '''Creates a non-vegetated pixel mask for an image.

     Args:
       Image(ee.Image): a given ee.Image object;
       SsrData(Dictionary): A Dictionary containing metadata associated with a sensor and data unit;
       MaxRef(int): a maximum reflectance value (1 or 100);
       indx_name(string): the name string of build-up index.'''

  red = Image.select(SsrData['RED'])
  nir = Image.select(SsrData['NIR'])  
  sw1 = Image.select(SsrData['SW1'])
  sw2 = Image.select(SsrData['SW2'])
  
  img0 = ee.Image(0)

  if indx_name.find('lxi') > -1:  # the index developed by Lixin Sun
    ndvi = Image.normalizedDifference([SsrData['NIR'], SsrData['RED']])
    nir_thresh = ee.Number(0.08).multiply(MaxRef)
    condition = ndvi.lt(ee.Image(0.3)).And(nir.gt(nir_thresh))   # NDVI < 0.1 and NIR > 8
  
    mask = ndvi.multiply(ee.Image(0))  
    return mask.where(condition, ee.Image(1))
    
  elif indx_name.find('nbi') > -1:  # (SWIR*RED)/NIR
    sw = sw1.add(sw2).divide(ee.Image(2))
    return sw.multiply(red).divide(nir)

  elif indx_name.find('ndbi') > -1: # (SWIR1 - NIR)/(SWIR1 + NIR)
    mask = sw1.subtract(nir).divide(sw1.add(nir))
    return mask.where(mask.lt(img0), img0)

  elif indx_name.find('bui') > -1: # NDBI - NDVI
    ndvi = nir.subtract(red).divide(nir.add(red))
    ndbi = sw1.subtract(nir).divide(sw1.add(nir))
    mask = ndbi.subtract(ndvi)
    
    return mask.where(mask.lt(img0).Or(ndvi.lt(img0)), img0)
  else:
    return red.multiply(img0)

  
 


#############################################################################################################
# Description: The function creates a water mask (1 ==> water) for an image.
#
# Revision history:  2021-Jun-13  Lixin Sun  Initial creation
#                    2022-Nov-15  Lixin Sun  Removed the limit that maximum reflectance value must be 100.
#
#############################################################################################################
def Img_WaterMask(Image, SsrData, MaxRef):
  '''This function creates a water mask (1 => water) for an image.

     Args:
       Image(ee.Image): a given ee.Image object;
       SsrData(Dictionary): A Dictionary containing metadata associated with a sensor and data unit;
       MaxRef(int): a maximum reflectance value (1 or 100).'''
 
  grn = Image.select(SsrData['GRN'])
  nir = Image.select(SsrData['NIR'])
  sw1 = Image.select(SsrData['SW1'])
  sw2 = Image.select(SsrData['SW2'])

  sw_mean = sw1.add(sw2).divide(2.0)

  # calculate NDWI map
  ndwi = grn.subtract(sw_mean).divide(grn.add(sw_mean))

  # create three tests for determining water pixels 
  sw_mean_thresh = ee.Number(0.02).multiply(MaxRef)
  nir_thresh1    = ee.Number(0.15).multiply(MaxRef)
  nir_thresh2    = ee.Number(0.10).multiply(MaxRef)

  test0 = sw_mean.lt(sw_mean_thresh).And(ndwi.gt(0.3))
  test1 = nir.lt(nir_thresh1).And(ndwi.gt(0.3))
  test2 = nir.lt(nir_thresh2).And(ndwi.gt(0.2))

  intrin_mask = Img_VenderMask(Image, SsrData, WATER_MASK)

  return test0.Or(test1).Or(test2).Or(intrin_mask)





#############################################################################################################
# Description: This function creates a snow pixel mask (1 => snow) for an image.
#
# Revision history:  2022-Jun-14  Lixin Sun  Initial creation
#                    2022-Nov-15  Lixin Sun  Removed the limit that maximum reflectance value must be 100.
# 
#############################################################################################################
def Img_SnowMask(Image, SsrData, MaxRef):
  '''Creates a snow/ice pixel mask (1 => snow/ice) for an image.

     Args:
       Image(ee.Image): a given ee.Image object;
       SsrData(Dictionary): A Dictionary containing metadata associated with a sensor and data unit;
       MaxRef(int): a maximum reflectance value (1 or 100).'''

  grn  = Image.select(SsrData['GRN'])  
  ndsi = Image.normalizedDifference([SsrData['GRN'], SsrData['SW1']])

  grn_thresh = ee.Number(0.1).multiply(MaxRef)
  data_mask = ee.Image(ndsi.gt(0.2).And(grn.gt(grn_thresh)))  
  
  intrin_mask = Img_VenderMask(Image, SsrData, SNOW_MASK)

  return data_mask.Or(intrin_mask)


 



#############################################################################################################
# Description: This function creates a valid pixel mask (mask out cloud, shadow, invalid value and saturated
#              pixels)for an image.
#
# Note:        The given image could be acquired either by a Sentinel-2 or a Landsat sensor.
#
# Revision history:  2020-Dec-20  Lixin Sun  Initial creation.
#                    2022-Jun-24  Lixin Sun  Modified according to the changes of the called functions.
#                    2022_Jul-26  Lixin Sun  Replaced "Img_VenderMask" for CLEAR_MASK with "Img_ClearMask", 
#                                            which includes customized cloud and shadow detections.
#############################################################################################################
def Img_ValidMask(Image, SsrData, MaxRef):
  '''Creates a valid pixel mask (mask out cloud, shadow, invalid value and saturated pixels) for an image.

     Args:
       Image(ee.Image): a given ee.Image object;
       SsrData(Dictionary): A Dictionary containing metadata associated with a sensor and data unit;
       MaxRef(int): a maximum reflectance value (1 or 100).'''

  clear_mask = Img_VenderMask(Image, SsrData, CLEAR_MASK)
  satur_mask = Img_VenderMask(Image, SsrData, SATU_MASK)
  value_mask = Img_ValueMask (Image, SsrData, MaxRef)

  return clear_mask.Or(satur_mask).Or(value_mask).rename(['ValidMask'])
  #return clear_mask.Or(value_mask).rename(['ValidMask'])




  

#############################################################################################################
# Description: This function creates a mask that mask out the land outside Canada and optionally water based
#              on a land cover map.
#
# Revision history:  2023-Feb-16  Lixin Sun  Initial creation.
#
#############################################################################################################
def Can_land_mask(Year, mask_water):
  '''Creates a mask that mask out the land outside Canada and optionally water.

     Args:      
       Year(int or string): A target year;
       mask_water(Boolean): Flag indicating if water bodies are masked out as well.'''
  #==========================================================================================================
  # Choose a proper land cover image collection based on a given "Year"
  #==========================================================================================================
  Can_LC = eoAD.get_CanLC(int(Year)).uint8()
  
  if mask_water == True:
    Can_LC = Can_LC.where(Can_LC.gt(17), ee.Image(0))  #class ID = 18 and 19 represent water and snow/ice, respectively

  #mask = ccrs_LC.selfMask()  # Mask ouit pixels with value equal to zero
  mask = Can_LC.where(Can_LC.gt(0), ee.Image(1))

  # Perform an erosion followed by a dilation
  box = ee.Kernel.circle(radius = 1, units = 'pixels', normalize = True)
  mask = mask.focalMin(kernel = box, iterations = 1).focalMax(kernel = box, iterations = 1)

  return mask.unmask()