//Exporting TS and attribute table, Reservoirs and Lakes with existing bathymetry
//now exporting existingIDs for new time periods, no more missingIDs  
// ****************************************************************************************************************** //
// USER INPUTS
// ****************************************************************************************************************** //
//geometries of water bodies
var Zarafshan_WaterBodies = ee.FeatureCollection("projects/ee-adriankreinerak/assets/Lakes_combination_v5");
var thisname='173003';//ID of the reservoir that will be used as an example for map display

//define the study period
var start_date = ee.Date.fromYMD(2017,1,1);
var end_date = ee.Date.fromYMD(2023,1,1);

//define asset paths and name
var SedZonesMapName= 'SedZone_00to22_';//sediment balance map export base name
var ts_name = 'TS_00to22_ID_';//base name of time series assets
var freqMapName='SurfWaterFreq_S2_17to22_';//frequency map export base name

//sediment balance maps asset path:
var path= 'projects/ee-adriankreinerak/assets/SedZones_Zarafshan/SedZones_1st_Iteration/';
//asset path for time-series (water levels and/or areas):
var pathTS='projects/ee-adriankreinerak/assets/TimeSeries';
//asset path for DEMs:
var path_dem = 'projects/ee-adriankreinerak/assets/DEMs_v2';
//asset path for frequency maps:
var path_Freq = 'projects/ee-adriankreinerak/assets/SurfFreqMaps_Lakes_combination_v2';

// ****************************************************************************************************************** //
// GET LIST OF IDs FOR SPECIFIC DATA TYPES
// ****************************************************************************************************************** //

var Table = ee.data.listAssets(path_dem);
var l =Table.assets.length;
var allids=[];
for (var i=0; i<l; i++) {//l
  allids.push(Table.assets[i].id.split('DEMfrom_'+freqMapName)[1]);
}
print('IDs with DEMs',allids);

var Table3 = ee.data.listAssets(pathTS + '/Landsat');
var l3 =Table3.assets.length;
var allids3=[];
for (var i=0; i<l3; i++) {//l
  allids3.push(Table3.assets[i].id.split(ts_name)[1]);
}
print('IDs with Landsat TS',allids3);

var Table4 = ee.data.listAssets(pathTS + '/Sentinel2');
var l4 =Table4.assets.length;
var allids4=[];
for (var i=0; i<l4; i++) {//l
  allids4.push(Table4.assets[i].id.split(ts_name)[1]);
}
print('IDs with Sentinel TS',allids4);

var Table = ee.data.listAssets(path);
var l =Table.assets.length;
var tableids=[];
for (var i=0; i<l; i++) {//l
  tableids.push(Table.assets[i].id.split(path + SedZonesMapName)[1].split('_1st_iteration')[0]);
}
print('IDs with Exported Sediment Balances',tableids);

var reservoirs_with_DEM=ee.FeatureCollection(ee.List(ee.List(tableids).iterate(function(x,previous){
  return ee.List(previous).add(Zarafshan_WaterBodies.filter(ee.Filter.eq('ID_new',ee.Number.parse(x))));
},ee.List([])))).flatten();
print('Reservoirs with DEMs',reservoirs_with_DEM);

// ****************************************************************************************************************** //
// IMPORT LAYERS AND CODES
// ****************************************************************************************************************** //
// Import external dependencies
var gee_codes = require('users/hydrosolutions/public_functions:SedimentBalance_T1L2');

//Import DEM
var dem_thislake = ee.Image(path_dem + "/DEMfrom_"  + freqMapName  + thisname);
//Import a shapefile of the area of interest
var aoi = Zarafshan_WaterBodies.filter(ee.Filter.eq('ID_new', ee.Number.parse(thisname)));
var thislake=aoi.first();
Map.centerObject(aoi);
//keep only the geometry, add a buffer
aoi=aoi.first().geometry();//.buffer(500) 500m buffer around polygon
Map.addLayer(ee.Image().paint(aoi, 0, 2),null, "area_of_interest");
Map.setOptions("HYBRID");

var FreqMapCollection=ee.ImageCollection(path_Freq).map(function(img){
    return ee.Image(img).set('ID',ee.String(ee.String(img.get('system:index')).split(freqMapName).get(1)));
  });
  
var FreqMap=FreqMapCollection.filter(ee.Filter.eq('ID',thisname)).first();
FreqMap=FreqMap.multiply(100).round().divide(100);
var maxFreq=ee.Number(FreqMap.rename('b1').reduceRegion(ee.Reducer.max(),aoi,10).get('b1'));

var minDEM=ee.Number(dem_thislake.get('min_elev_perc_fit')).getInfo();
var maxDEM=ee.Number(dem_thislake.get('max_elev_perc_fit')).getInfo();
Map.addLayer(FreqMap,{min:0,max:1,palette:'white,blue'}, 'water surface probability S2',false);

var aoi_img=ee.Image(0).clip(aoi.buffer(500)).blend(ee.Image(1).clip(aoi)).float();
var dem_orig=dem_thislake;
var dem_defined=dem_thislake.updateMask(dem_thislake.lte(maxDEM)).updateMask(dem_thislake.gte(minDEM)).clip(aoi);
//var dem_filled = ee.Image(-99).rename('b1').where(dem_defined.gt(0),dem_defined).clip(aoi);
var slope_occ=(ee.Number(dem_thislake.get('max_elev_perc_fit')).subtract(ee.Number(dem_thislake.get('min_elev_perc_fit')))).divide(
  ee.Number(100).multiply(ee.Number(dem_thislake.get('max_occ')).subtract(ee.Number(dem_thislake.get('min_occ')))));
// print('slope_occ',slope_occ);
var dem_delta = FreqMap.subtract(ee.Image(ee.Number(dem_thislake.get('max_occ')))).multiply(100).multiply(slope_occ).clip(aoi);
var dem_filled=dem_orig.clip(aoi).where(dem_delta.gt(0),ee.Image(ee.Number(dem_thislake.get('min_elev_perc_fit'))).subtract(dem_delta))
  .updateMask(FreqMap.gte(ee.Number(dem_thislake.get('min_occ'))));
var dem_reservoir=dem_orig.where(dem_delta.gt(0),ee.Image(ee.Number(dem_thislake.get('min_elev_perc_fit'))).subtract(dem_delta));
var dem_reservoir2=dem_orig.where(dem_delta.gt(0),ee.Image(ee.Number(dem_thislake.get('min_elev_perc_fit'))).subtract(dem_delta))
  .where(FreqMap.lt(0.05),-99)//lets say where the occ probab is less than 5% this is outside the reservoir
  .where(aoi_img.eq(0),-99);//clip away those areas outside aoi with occprob >=5%? (it could be rice fields, ponds, etc..)

Map.addLayer(dem_thislake.updateMask(dem_thislake.lte(maxDEM)).updateMask(dem_thislake.gte(minDEM)), {min:minDEM,max:maxDEM,palette:'white,black'},'DEM',false);


// ****************************************************************************************************************** //
//ENSEMBLE OF LANDSAT AND SENTINEL-2 SCENES
// ****************************************************************************************************************** //
// Import external dependencies
var gee_codes = require('users/hydrosolutions/public_functions:SedimentBalance_T1L2');
// number of years
var years = end_date.difference(start_date,'year').round();
print('number of years considered for analysis',years);

var merged_landsat = ee.ImageCollection(gee_codes.get_landsat_images(aoi.buffer(500),start_date,end_date,
  {nodata_thresh: 30,
  cc_thresh: 30,
  nodata_thresh_l7: 50,
  })).sort('system:time_start');
print('Number of scenes in collection:',merged_landsat.size());

Map.addLayer(ee.Image(merged_landsat.toList(9999).get(0)), {bands: (['R', 'G', 'B']), min: 0, max: 0.3}, 'RGB Landsat',false);

 var s2_all = ee.ImageCollection(gee_codes.get_s2_images(aoi.buffer(500),start_date,end_date,
    {nodata_thresh: 30,
    cc_thresh: 30,
    cc_pix: 50,
    }));
  var edge_imgs_s2 = ee.ImageCollection(s2_all
      //returns the DEM clipped to the shorelines for each scene
       .map(function(img){return gee_codes.edgeOtsu_optical(img,dem_reservoir,
        {cannyThreshold: 0.7,
          cannySigma: 0.7,
          minimumValue: -0.3,
          maximumValue: 0.3,
          scale: 30,
          th_veg: 0.5,
          withDEM: true,
          keepMNDWI:true
        })})).sort('system:time_start'); 
Map.addLayer(ee.Image(s2_all.toList(9999).get(0)), {bands: (['R', 'G', 'B']), min: 0, max: 3000}, 'RGB Sentinel-2',false);
       
// ****************************************************************************************************************** //
//SHORELINES OF ALL LANDSAT SCENES
// ****************************************************************************************************************** //

var edge_imgs_landsat = ee.ImageCollection(merged_landsat
    //returns the DEM clipped to the shorelines for each scene
      .map(function(img){return gee_codes.edgeOtsu_optical(img,dem_reservoir,
        {cannyThreshold: 0.7,
          cannySigma: 0.7,
          minimumValue: -0.3,
          maximumValue: 0.3,
          scale: 30,
          th_veg: 0.5,
          withDEM: true,
          keepMNDWI:true
        })})).sort('system:time_start');
// print('edge_imgs_landsat',edge_imgs_landsat);

Map.addLayer(ee.Image(edge_imgs_landsat.toList(999).get(0)).updateMask(ee.Image(edge_imgs_landsat.toList(999).get(0)).select('veg_edge').neq(1)),{bands: (['dem_edge']),palette:[ 'white','blue'],min: -1, max: 1}, "first dem_edge Landsat",false);
Map.addLayer(ee.Image(edge_imgs_s2.toList(999).get(0)).updateMask(ee.Image(edge_imgs_s2.toList(999).get(0)).select('veg_edge').neq(1)),{bands: (['dem_edge']),palette:[ 'white','blue'],min: -1, max: 1}, "first dem_edge Sentinel",false);
//Map.addLayer(ee.Image(edge_imgs_landsat.toList(999).reverse().get(6)).updateMask(ee.Image(edge_imgs_landsat.toList(999).reverse().get(6)).select('veg_edge').neq(1)),{bands: (['dem_edge']),palette:[ 'white','blue'],min: -1, max: 1}, "sixth dem_edge",false);


var func4edgeOtsu = function(img){return gee_codes.edgeOtsu_optical(img,dem_reservoir,
          {cannyThreshold: 0.7,
            cannySigma: 0.7,
            minimumValue: -0.3,
            maximumValue: 0.3,
            scale: 30,
            th_veg: 0.5,
            withDEM: true,
            keepMNDWI:true
          })};

var existingIDs=[];
for (var i=0; i<tableids.length; i++) {//l
  if( allids.indexOf(tableids[i])>=0){
    existingIDs.push(tableids[i]);
  } 
}
allids=existingIDs;

var missingIDs=[];
for (var i=0; i<allids.length; i++) {//l
  if( tableids.indexOf(allids[i])<0){
    missingIDs.push(allids[i]);
  }
}
print('missing IDs SedZones', missingIDs);

var missingIDs=[];
for (var i=0; i<allids.length; i++) {//l
  if( allids3.indexOf(allids[i])<0){
    missingIDs.push(allids[i]);
  }
}

print('missing IDs TS Landsat', missingIDs);

var missingIDs=[];
var existingIDs=[];
for (var i=0; i<allids.length; i++) {//l
  if( allids4.indexOf(allids[i])<0){
    missingIDs.push(allids[i]);
    //print('Missing Lake ID :'+allids[i]);
  } else {
    existingIDs.push(allids[i]);
  }
}

print('missing IDs TS Sentinel', missingIDs);
// print('DEMs with Sediment Balances and Time Series', existingIDs);

// ****************************************************************************************************************** //
//CALCULATE SEDIMENT BALANCE PER PIXEL WITHIN AOI
//sediment balance is required for TS generation, because pixels with unstable terrain will be excluded from the shoreline considered to measure water levels
// ****************************************************************************************************************** //
//if already exported, load the annual shoreline anomalies from the asset
var annual_deltah=ee.Image(path +  SedZonesMapName + thisname + '_1st_iteration');
//Sen's slope of shoreline anomalies per pixel --> sediment balance
var sed_thislake=ee.Image(gee_codes.stacktoimage(annual_deltah));

Map.addLayer(sed_thislake,{bands: (['slope']),palette:['red', 'white','blue'],min: -0.05, max: 0.05}, "slope initial");

function get_new_ts_toexport(aoi,thisname,sed_thislake,edge_imgs_landsat,aoi_img){
  sed_thislake=ee.Image(sed_thislake);
  aoi=ee.Geometry(aoi);
  thisname=ee.String(thisname);
  edge_imgs_landsat=ee.ImageCollection(edge_imgs_landsat).map(function(img){return ee.Image(img).addBands(ee.Image(img).select('dem_edge')/*.updateMask(river_img.neq(1))*/ ,null,true)});
  var edge_img4DEM = ee.ImageCollection(gee_codes.get_corrected_elevation_DEM(edge_imgs_landsat, sed_thislake, {
    filter_outliers: true, //optional: filter out the noisiest shorelines (dh_std >= 1.960 standard deviations of the entire ensemble)
    pvalues: true, //optional: consider only pixels with p-values <0.1
    additional_exports: true, //calculate also the water surface area and additional properties
    },aoi_img));

  var edge_imgs_landsat_merged=edge_img4DEM;
    //.filter(ee.Filter.gt('dem_elev', -100));
    //.filter(ee.Filter.eq('SATELLITE', 'LANDSAT_5'));
    
  var mndwi_thresh=-0.3;
  var elev_data=edge_imgs_landsat_merged.filter(ee.Filter.neq('otsu_th', mndwi_thresh));

  //return features without geometry
  var export_fc = ee.FeatureCollection(elev_data.map(function(feature){
  var export_feature = ee.Feature(aoi).centroid(1000).copyProperties(feature,null,dem_thislake.propertyNames())
    .set('lake_id',feature.get('lake_id'))
    .set('system:time_start',feature.get('system:time_start'));
    return export_feature }));  
  return export_fc.sort('system:time_start');
}

var new_table=get_new_ts_toexport(aoi,thisname,sed_thislake,edge_imgs_landsat,aoi_img);
// print('new_table',new_table);

var s2_all = ee.ImageCollection(gee_codes.get_s2_images(aoi,start_date,end_date,
  {nodata_thresh: 30,
  cc_thresh: 30,
  cc_pix: 50,
  }));
var edge_imgs_s2 = ee.ImageCollection(s2_all
    //returns the DEM clipped to the shorelines for each scene
     .map(func4edgeOtsu)); 
 var new_table_s2=get_new_ts_toexport(aoi.buffer(500),thisname,sed_thislake,edge_imgs_s2,aoi_img);
// print('new_table_s2',new_table_s2);

new_table_s2=ee.FeatureCollection(new_table_s2.map(function(ft){
  return ee.Feature(null).copyProperties(ft,['dem_elev'])
    .set('s2_elev',ft.get('dem_elev')).set('system:time_start',ft.get('system:time_start'))}));

print(ui.Chart.feature.byFeature(ee.FeatureCollection(new_table), 'system:time_start',['area_water']).setOptions({
        lineWidth: 1,
        pointSize: 2,
        legend: {position:"top"},
        fontSize: 14,
        title: 'Water Areas based on Landsat',//(no outlier cleaning)
        series: {
          0: {targetAxisIndex: 0,color: '#ff0000'},
        },
        vAxes: {
          0: {title: 'Area [ha]'}, //00FF00
        },
      }))

print(ui.Chart.feature.byFeature(ee.FeatureCollection(new_table), 'system:time_start',['dem_elev']).setOptions({
        lineWidth: 1,
        pointSize: 2,
        legend: {position:"top"},
        fontSize: 14,
        title: 'Water Levels based on Landsat',//(no outlier cleaning)
        series: {
          0: {targetAxisIndex: 0,color: '#ff0000'},
        },
        vAxes: {
          0: {title: 'Elevation [m asl]'}, //00FF00
        },
      }));
      
print(ui.Chart.feature.byFeature(ee.FeatureCollection([new_table_s2]).flatten().sort('system:time_start'), 'system:time_start',['s2_elev']).setOptions({
        lineWidth: 1,
        pointSize: 2,
        legend: {position:"top"},
        fontSize: 14,
        title: 'Water Levels based on Sentinel-2',//(no outlier cleaning)
        series: {
          0: {targetAxisIndex: 0,color: '#ff0000'},
          1: {targetAxisIndex: 0,color: 'black'}
        },
        vAxes: {
          0: {title: 'Elevation [m asl]'}, //00FF00
        },
      }));

print("missingIDs for loop", missingIDs);

missingIDs=reservoirs_with_DEM.aggregate_array('ID_new').getInfo();

for (var i=0; i<missingIDs.length; i++) {
  var thisname=missingIDs[i].toString(); 
  var dem_thislake = ee.Image(path_dem + "/DEMfrom_"  + freqMapName + thisname);
  //Import a shapefile of the area of interest
  var aoi = Zarafshan_WaterBodies.filter(ee.Filter.eq('ID_new',ee.Number.parse(thisname)));
  var thislake=ee.Feature(aoi.first());
  aoi=aoi.first().geometry();//.buffer(500) 500m buffer around polygon

  var aoi_img=ee.Image(0).clip(aoi.buffer(500)).blend(ee.Image(1).clip(aoi)).float();
  
  var FreqMap=FreqMapCollection.filter(ee.Filter.eq('ID',thisname)).first();

  var annual_deltah=ee.Image(path +  SedZonesMapName + thisname + '_1st_iteration');
  //Sen's slope of shoreline anomalies per pixel --> sediment balance
  var sed_thislake=ee.Image(gee_codes.stacktoimage(annual_deltah));
  // ****************************************************************************************************************** //
  //TIME SERIES
  // ****************************************************************************************************************** //    

  var merged_landsat = ee.ImageCollection(gee_codes.get_landsat_images(aoi.buffer(500),start_date,end_date,
    {nodata_thresh: 30,
    cc_thresh: 30,
    nodata_thresh_l7: 50,
    }));
  var edge_imgs_landsat = ee.ImageCollection(merged_landsat
      //returns the DEM clipped to the shorelines for each scene
       .map(func4edgeOtsu)); 
       
  var new_table=get_new_ts_toexport(aoi,thisname,sed_thislake,edge_imgs_landsat,aoi_img);
  //print('new_table',new_table);
    Export.table.toAsset({ 
      collection: new_table,
      description: 'Landsat_TS_17to22_ID_' + thisname,
      assetId: pathTS + '/Landsat/TS_17to22_ID_' + thisname,
    });
   /* Export.table.toDrive({
      collection: new_table,
      description: 'Landsat_TS_17to22_ID_' + thisname,
      folder: 'Landsat',
    });*/
    
  var s2_all = ee.ImageCollection(gee_codes.get_s2_images(aoi.buffer(500),start_date,end_date,
    {nodata_thresh: 30,
    cc_thresh: 30,
    cc_pix: 50,
    }));
  var edge_imgs_s2 = ee.ImageCollection(s2_all
      //returns the DEM clipped to the shorelines for each scene
       .map(func4edgeOtsu)); 
  var new_table_s2=get_new_ts_toexport(aoi,thisname,sed_thislake,edge_imgs_s2,aoi_img);
  // print('new_table_s2',new_table_s2);
  Export.table.toAsset({ 
    collection: new_table_s2,
    description: 'Sentinel2_TS_17to22_ID_' + thisname,
    assetId: pathTS+ '/Sentinel2/TS_17to22_ID_' + thisname,
  });
/*  Export.table.toDrive({
        collection: new_table_s2,
        description: 'Sentinel2_TS_17to22_ID_' + thisname,
        folder: 'Sentinel2',
      });*/
}