// ------This code is used to extract burned areas and calculate RSEI in GEE
// ------Author: Shiqi Zhang, E-mail: zhangshiqiiiii@163.com
// ------NBR
function NBR(img){
 var nbr = img.expression(
   '(NIR-SWIR)/(NIR+SWIR)',
    {
      'NIR': img.select('B8A'),
      'SWIR': img.select('B12')
    }
  );
  return nbr;
}

// -----OTSU
var otsu = function(histogram) {
  var counts = ee.Array(ee.Dictionary(histogram).get('histogram')); 
  var means = ee.Array(ee.Dictionary(histogram).get('bucketMeans')); 
  var size = means.length().get([0]); 
  var total = counts.reduce(ee.Reducer.sum(), [0]).get([0]); 
  var sum = means.multiply(counts).reduce(ee.Reducer.sum(), [0]).get([0]); 
  var mean = sum.divide(total); 
  var indices = ee.List.sequence(1, size); 
  var bss = indices.map(function(i) { 
    var aCounts = counts.slice(0, 0, i); 
    var aCount = aCounts.reduce(ee.Reducer.sum(), [0]).get([0]); 
    var aMeans = means.slice(0, 0, i); 
    var aMean = aMeans.multiply(aCounts) 
      .reduce(ee.Reducer.sum(), [0]).get([0]) 
      .divide(aCount); 
    var bCount = total.subtract(aCount); 
    var bMean = sum.subtract(aCount.multiply(aMean)).divide(bCount); 
    return aCount.multiply(aMean.subtract(mean).pow(2)).add( 
          bCount.multiply(bMean.subtract(mean).pow(2)));
  }); 
  return means.sort(bss).get([-1]);
};

// ------RSEI
function Wet(img){
  return img.addBands(img.expression('B*(0.1509) + G*(0.1973) + R*(0.3279) + NIR*(0.3406) + SWIR1*(-0.7112) + SWIR2*(-0.4572)',{
       'B': img.select(['SR_B2']).multiply(0.0001),
       'G': img.select(['SR_B3']).multiply(0.0001),
       'R': img.select(['SR_B4']).multiply(0.0001),
       'NIR': img.select(['SR_B5']).multiply(0.0001),
       'SWIR1': img.select(['SR_B6']).multiply(0.0001),
       'SWIR2': img.select(['SR_B7']).multiply(0.0001)
     }).rename('Wet'))
}

function TOA_LST(img){
  var L8_B10 = img.select('B10').multiply(0.1)
  var ndvi =  img.normalizedDifference(['B5', 'B4'])
  var pv = (ndvi.subtract(ee.Image(0.05)).divide((ee.Image(0.95)).subtract(ee.Image(0.05)))).pow(ee.Number(2))
  var c = pv.expression("0.004*pv+0.986", {pv: pv})
  var LST = L8_B10.expression(
'(L8_B10/(1 + (0.00109* (L8_B10 / 1.438))*log(c)))-273.15', {
 L8_B10: L8_B10,
    c: c,
})
  return img.addBands(LST.rename('LST_TOA'))
}

function si(img){
  var si = img.expression('((SWIR1 + RED) - (NIR + BLUE)) / ((SWIR1 + RED) + (NIR + BLUE))', {
      'SWIR1': img.select('SR_B6').multiply(0.0001),
      'NIR': img.select('SR_B5').multiply(0.0001),
      'RED': img.select('SR_B4').multiply(0.0001),
      'BLUE': img.select('SR_B2').multiply(0.0001)
    }) 
  return img.addBands(si.rename('si')) 
}

function L8_RSEI(img){
  return L8_Wet(img).select('Wet').addBands(L8_ndvi(img).select('ndvi'))
              .addBands(L8_si(img).select('si'))
              .copyProperties(img,img.propertyNames())
}

// ------PCA
var PCA = function(centered, scale, region) {
  var bandNames = centered.bandNames();
  var getNewBandNames = function(prefix) {
  var seq = ee.List.sequence(1, bandNames.length());
    return seq.map(function(b) {
      return ee.String(prefix).cat(ee.Number(b).int());
    });
  };
  var arrays = centered.toArray();
  var covar = arrays.reduceRegion({
    reducer: ee.Reducer.centeredCovariance(),
    geometry: region,
    scale: scale,
    maxPixels: 1e9
  });
  var covarArray = ee.Array(covar.get('array'));
  var eigens = covarArray.eigen();
  var eigenValues = eigens.slice(1, 0, 1);
  var eigenValuesList = eigenValues.toList().flatten()
  var total = eigenValuesList.reduce(ee.Reducer.sum())
  var percentageVariance = eigenValuesList.map(function(item) {
    return (ee.Number(item).divide(total)).multiply(100).format('%.2f')
  })
  var eigenVectors = eigens.slice(1, 1);
  var arrayImage = arrays.toArray(1);
  var principalComponents = ee.Image(eigenVectors).matrixMultiply(arrayImage);
  var sdImage = ee.Image(eigenValues.sqrt())
    .arrayProject([0]).arrayFlatten([getNewBandNames('sd')]);
  principalComponents=principalComponents
    .arrayProject([0])
    .arrayFlatten([getNewBandNames('pc')])
    .divide(sdImage);
  return principalComponents
};

// ------Get Image
function cloudfree_landsat (image){
  var qa = image.select('QA_PIXEL');
  var cloudsBitMask = 1 << 3;
  var cloudShadowBitMask = 1 << 4;
  var snowBitMask = 1 << 5;
  var mask = qa.bitwiseAnd(cloudsBitMask).eq(0)
              .and(qa.bitwiseAnd(cloudShadowBitMask).eq(0))
              .and(qa.bitwiseAnd(snowBitMask).eq(0));
  return image.updateMask(mask)  
}

var image_pre = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")
              .filterDate('2016-04-20', '2016-05-20')
              .filterBounds(roi)
              .map(cloudfree_landsat)
              .median()
              .clip(roi);
              
var image_post = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")
              .filterDate('2020-04-20', '2020-05-20')
              .filterBounds(roi)
              .map(cloudfree_landsat)
              .median()
              .clip(roi);

//TOA images
var image_toa_pre = ee.ImageCollection("LANDSAT/LC08/C02/T1_TOA")
              .filterDate('2016-04-20', '2016-05-20')
              .filterBounds(roi)
              .map(cloudfree_landsat)
              .median()
              .clip(roi);
              
var image_toa_post = ee.ImageCollection("LANDSAT/LC08/C02/T1_TOA")
              .filterDate('2020-04-20', '2020-05-20')
              .filterBounds(roi)
              .map(cloudfree_landsat)
              .median()
              .clip(roi);

function maskS2clouds(image) {
  var qa = image.select('QA60');

  // Bits 10 and 11 are clouds and cirrus, respectively.
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;

  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
      .and(qa.bitwiseAnd(cirrusBitMask).eq(0));

  return image.updateMask(mask).divide(10000);
}

var image_pre = ee.ImageCollection("COPERNICUS/S2")
              .filterDate('2019-04-20', '2019-05-20')
              .filterBounds(roi)
              .map(maskS2clouds)
              .median()
              .clip(roi);
              
var image_post = ee.ImageCollection("COPERNICUS/S2")
              .filterDate('2020-04-20', '2020-05-20')
              .filterBounds(roi)
              .map(maskS2clouds)
              .median()
              .clip(roi);

// -----Extract burned areas
//dNBR
var image_pre_burn = image_pre.addBands(NBR(image_pre).rename('NBR'));
var image_post_burn = image_post.addBands(NBR(image_post).rename('NBR'));

//Layer1
var dNBR = image_pre_burn.select("NBR").subtract(image_post_burn.select("NBR")).rename("dNBR");
Map.addLayer(dNBR, {}, 'dNBR', false);
var histogram_dnbr = dNBR.select('dNBR').reduceRegion({
  reducer: ee.Reducer.histogram(), 
  geometry: roi, 
  scale: 20, 
  maxPixels: 1e13
});

var dNBR_threshold = otsu(ee.Dictionary(histogram_dnbr).get('dNBR')); 
print('dNBR_threshold ', dNBR_threshold);

var mask2 = dNBR.gt(dNBR_threshold);
var Burned_area = mask2;
Burned_area = Burned_area.updateMask(mask2);
Map.addLayer(Burned_area, {palette: ['#FF0000']}, "Burned_area");

//Layer2
image_pre_burn = image_pre_burn.updateMask(Burned_area.mask().not());
image_post_burn = image_post_burn.updateMask(Burned_area.mask().not());
var dNBR_lighter = image_pre_burn.select("NBR").subtract(image_post_burn.select("NBR")).rename("dNBR_lighter");
Map.addLayer(dNBR_lighter, {}, 'dNBR_lighter', false);
var histogram_dnbr_lighter = dNBR_lighter.select('dNBR_lighter').reduceRegion({
  reducer: ee.Reducer.histogram(), 
  geometry: roi, 
  scale: 20, 
  maxPixels: 1e13
});
var dNBR_threshold_lighter = otsu(ee.Dictionary(histogram_dnbr_lighter).get('dNBR_lighter'));
print('dNBR_threshold_lighter', dNBR_threshold_lighter);

var mask3 = dNBR_lighter.gt(dNBR_threshold_lighter);
var Burned_area_lighter = mask3;
Burned_area_lighter = Burned_area_lighter.updateMask(mask3);
Map.addLayer(Burned_area_lighter, {palette: ['#FFA500']}, "Burned_area_lighter");

//download images
Export.image.toDrive({
  image: Burned_area,
  description: 'Heavy_fire_area_2020',
  scale: 30,
  region: roi,
  maxPixels: 1e13
});

Export.image.toDrive({
  image: Burned_area_lighter,
  description: 'Mild_fire_area_2020',
  scale: 30,
  region: roi,
  maxPixels: 1e13
});

// -----Calculate RSEI
var img_normalize = function(img){
      var minMax = img.reduceRegion({
            reducer:ee.Reducer.minMax(),
            geometry: roi,
            scale: 30,
            maxPixels: 1e13,
            tileScale: 16
        })
      var year = img.get('year')
      var normalize  = ee.ImageCollection.fromImages(
            img.bandNames().map(function(name){
                  name = ee.String(name);
                  var band = img.select(name);
                  return band.unitScale(ee.Number(minMax.get(name.cat('_min'))), ee.Number(minMax.get(name.cat('_max'))));
                    
              })
        ).toBands().rename(img.bandNames());
        return normalize;
}

var RSEI_pre = ee.Image(L8_RSEI(image_pre));
var RSEI_post = ee.Image(L8_RSEI(image_post));
var TOA_LST_pre = ee.Image(L8_TOA_LST(image_toa_pre));
var TOA_LST_post = ee.Image(L8_TOA_LST(image_toa_post));
RSEI_pre = RSEI_pre.addBands(TOA_LST_pre.select('LST_TOA'));
RSEI_post = RSEI_post.addBands(TOA_LST_post.select('LST_TOA'));
RSEI_pre = img_normalize(RSEI_pre).addBands(ee.Image.constant(2020).int8().rename('time')).set({'year': 2019});
RSEI_post = img_normalize(RSEI_post).addBands(ee.Image.constant(2021).int8().rename('time')).set({'year': 2020});
Map.addLayer(RSEI_pre,{},'RSEI_pre', false);
Map.addLayer(RSEI_post,{},'RSEI_post', false);

var pre_years = ee.List.sequence(2017, 2019);
var pre_list = pre_years.map(function(year) {
  var L8 = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2").filterBounds(roi)
                .filter(ee.Filter.calendarRange(year, year, 'year'))
                .filter(ee.Filter.calendarRange(4, 6, 'month'))
                .map(cloudfree_landsat).map(L8_RSEI).median().clip(roi);
  var L8_TOA = ee.ImageCollection("LANDSAT/LC08/C02/T1_TOA").filterBounds(roi)
                .filter(ee.Filter.calendarRange(year, year, 'year'))
                .filter(ee.Filter.calendarRange(4, 6, 'month'))
                .map(cloudfree_landsat).map(L8_TOA_LST).median().clip(roi);
  L8 = L8.addBands(L8_TOA.select('LST_TOA'));
  return img_normalize(L8).addBands(ee.Image.constant(year).int8().rename('time')).set({'year': year});
  })
var pre_images = ee.ImageCollection(pre_list).merge(RSEI_pre);

var pre_Interpol = require("users/cannanshiki/Function:RSEI_interpolation").Interpol;
var pre_interpolate = pre_Interpol(pre_images, 'time').map(function(image){return image.select(['Wet','ndvi','LST_TOA','si'])}).toList(4);
print(pre_interpolate);
Map.addLayer(ee.Image(pre_interpolate.get(3)), {}, 'RSEI_pre_interpol');
var pre_rsei_before = ee.Image(pre_interpolate.get(3));

var post_years = ee.List.sequence(2020, 2022);
var post_list = post_years.map(function(year) {
  var L8 = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2").filterBounds(roi)
                .filter(ee.Filter.calendarRange(year, year, 'year'))
                .filter(ee.Filter.calendarRange(4, 6, 'month'))
                .map(cloudfree_landsat).map(L8_RSEI).median().clip(roi);
  var L8_TOA = ee.ImageCollection("LANDSAT/LC08/C02/T1_TOA").filterBounds(roi)
                .filter(ee.Filter.calendarRange(year, year, 'year'))
                .filter(ee.Filter.calendarRange(4, 6, 'month'))
                .map(cloudfree_landsat).map(L8_TOA_LST).median().clip(roi);
  L8 = L8.addBands(L8_TOA.select('LST_TOA'));
  return img_normalize(L8).addBands(ee.Image.constant(year).int8().rename('time')).set({'year': year});
  })
var post_images = ee.ImageCollection(post_list).merge(RSEI_post);

var post_Interpol = require("users/cannanshiki/Function:RSEI_interpolation").Interpol;
var post_interpolate = post_Interpol(post_images, 'time').map(function(image){return image.select(['Wet','ndvi','LST_TOA','si'])}).toList(4);
print(post_interpolate);
Map.addLayer(ee.Image(post_interpolate.get(3)), {}, 'RSEI_post_interpol');
var post_rsei_before = ee.Image(post_interpolate.get(3));

var RSEI_pre_PCA = PCA(ee.Image(pre_rsei_before), 30, roi);
var rsei_pre = ee.Image(ee.Image(0).expression('constant - pc1' , {
            'constant': ee.Image(1),
            pc1: ee.Image(RSEI_pre_PCA).select('pc1')
        })).rename('rsei');
rsei_pre = img_normalize(rsei_pre);
Map.addLayer(ee.Image(rsei_pre), {}, 'PCA_RSEI_pre', false);

var RSEI_post_PCA = PCA(ee.Image(post_rsei_before), 30, roi);
var rsei_post = ee.Image(ee.Image(0).expression('constant - pc1' , {
            'constant': ee.Image(1),
            pc1: ee.Image(RSEI_post_PCA).select('pc1')
        })).rename('rsei');
rsei_post = img_normalize(rsei_post);
Map.addLayer(ee.Image(rsei_post), {}, 'PCA_RSEI_post', false);

//download images
Export.image.toDrive({
  image: ee.Image(rsei_pre),
  description: 'PCA_RSEI_pre_2020',
  scale: 30,
  region: roi,
  maxPixels: 1e13
});

Export.image.toDrive({
  image: ee.Image(rsei_post),
  description: 'PCA_RSEI_post_2020',
  scale: 30,
  region: roi,
  maxPixels: 1e13
});