sentinle-1-2021a-4splits

// Function returns true if leap year, false otherwise.
var isLeapYear = function(year) {
  return (year % 4 === 0 && year !== 100) || year % 400 === 0;
};

// Gets the number of days in any month in a year.
var getDaysInMonth = function(year, month) {
  var leapYear = isLeapYear(parseInt(year, 10));
  var daysInMonth = 0;
  var monthNumber = parseInt(month, 10);

  if (monthNumber === 2) {
    daysInMonth = 28;
    //leapyear ? daysInMonth = 29 : daysInMonth = 28;
  }
  else if (
    monthNumber === 4 ||
    monthNumber === 6 ||
    monthNumber === 9 ||
    monthNumber === 11
  ){
    daysInMonth = 30;
  }
  else {
    daysInMonth = 31;
  }
  
  return (daysInMonth);
};


// Adds rectangle (the study area) and create a Feature object using the rectangle.
// var region = ee.Geometry.Rectangle(-5.0,47.2,-1.0,49.0); //-4.82,48.1,-3.4,48.9
var region = ee.Geometry.Rectangle(-4.4,48.3,-3.8,48.7)
var studyArea = ee.Feature(region, { name: 'france'});


/* Show object properies in the Console. */
print(studyArea);

/* Sets the map view to the center
 * of the study area. */
//Map.setCenter(-3.5,48.5,9);
//Map.centerObject(studyArea);


// Loads the Sentinel-1 image collection.
var sentinel1 = ee.ImageCollection('COPERNICUS/S1_GRD');

/* Coordinate pairs, set start and end date
 * for filtering the collection. */
var point = ee.Geometry.Point(-3.5, 48.5);

var currentYear = new Date().getFullYear();

// User inputs: year and month. Data is available from 12/2014!
var year = prompt('Which year you want to process (between 2014 and ' + currentYear + ')?');
// Use leading 0 if month number < than 10!
var month = prompt('Which month you want to process?', '07');

var daysInMonth = getDaysInMonth(year, month);

var start = ee.Date(year + '-' + month + '-01');
var finish = ee.Date(year + '-' + month + '-' + daysInMonth);
print(finish);
print(daysInMonth);

 // Filtering based on metadata properties.
var vh = sentinel1
  // Filter to get images with VV and VH dual polarization.
  .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
  .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH'))
  // Filter to get images collected in interferometric wide swath mode.
  .filter(ee.Filter.eq('instrumentMode', 'IW'))
  .filterDate(start, finish);
  
var TER = 0;
var tertest = '0';
if (TER===1) {
  vh = vh.map(terrainCorrection);
  print('terrain correction');
  tertest = '1';
}
if (TER===0) {
  vh = vh;
  print('No terrain correction');
  tertest = '0';
}

// Filter to get images from different look angles.
var vhAscending = vh.filter(ee.Filter.eq('orbitProperties_pass', 'ASCENDING'));
var vhDescending = vh.filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING'));


// Create a composite from means at different polarizations and look angles.
var composite = ee.Image.cat([
  vhAscending.select('VV').mean(),
  vhDescending.select('VH').mean(),
  vhDescending.select('VV').mean()
]).focal_median();

// 4 splits
var compositeAsVH = ee.Image.cat([
  vhAscending.select('VH').mean()
]).focal_median();
var compositeClippedAsVH = compositeAsVH.clip(studyArea).multiply(10000).toInt();

var compositeAsVV = ee.Image.cat([
  vhAscending.select('VV').mean()
]).focal_median();
var compositeClippedAsVV = compositeAsVV.clip(studyArea).multiply(10000).toInt();

var compositeDesVH = ee.Image.cat([
  vhDescending.select('VH').mean()
]).focal_median();
var compositeClippedDesVH = compositeDesVH.clip(studyArea).multiply(10000).toInt();

var compositeDesVV = ee.Image.cat([
  vhDescending.select('VV').mean()
]).focal_median();
var compositeClippedDesVV = compositeDesVV.clip(studyArea).multiply(10000).toInt();

// Clip composite image with our study area.
var compositeClipped = composite.clip(studyArea);

// Display as a composite of polarization and backscattering characteristics.
Map.addLayer(
  compositeClipped,
  {
    min: [-25, -20, -25],
    max: [0, 10, 0]
  },
  'composite'
);

// Create an empty image into which to paint the features, cast to byte.
var empty = ee.Image().byte();

// Paint all the polygon edges with the same number and width, display.
var outline = empty.paint(
  {
    featureCollection: studyArea,
    color: "black",
    width: 2
  }
);

Map.addLayer(outline, { palette: 'f00' }, 'Study Area');
print(compositeClipped);

var composite_out = ee.Image.cat([
  compositeClippedAsVH,
  compositeClippedAsVV,
  compositeClippedDesVH,
  compositeClippedDesVV
])

// Save the composite image as GeoTIFF.
Export.image.toDrive(
  {
    image: composite_out,
    description: 'Sen1_terCor'+tertest+'_asVH_asVV_desVH_desVV'+year+month,
    scale: 10,
    fileFormat: 'GeoTIFF',
    maxPixels: 3784216672400,
    region: studyArea
  }
);

// 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'], ['VV', 'VH']),
    null,
    true
  );
}
 
 
function powerToDb(img){
  return ee.Image(10).multiply(img.log10());
}
 
function dbToPower(img){
  return ee.Image(10).pow(img.divide(10));
}
 
function gammaMap(img){
 
  var ksize = 7;
  var enl = 5;
  var bandNames = img.bandNames();
   
  // Convert image from dB to natural values
  var nat_img = dbToPower(img);
 
  // Square kernel, ksize should be odd (typically 3, 5 or 7)
  var weights = ee.List.repeat(ee.List.repeat(1,ksize),ksize);
   
  // ~~(ksize/2) does integer division in JavaScript
  var kernel = ee.Kernel.fixed(ksize,ksize, weights, ~~(ksize/2), ~~(ksize/2), false);
 
  // Get mean and variance
  var mean = nat_img.reduceNeighborhood(ee.Reducer.mean(), kernel);
  var variance = nat_img.reduceNeighborhood(ee.Reducer.variance(), kernel);
 
  // "Pure speckle" threshold
  var ci = variance.sqrt().divide(mean);  // square root of inverse of enl
 
  // If ci <= cu, the kernel lies in a "pure speckle" area -> return simple mean
  var cu = 1.0/Math.sqrt(enl);
   
  // If cu < ci < cmax the kernel lies in the low textured speckle area -> return the filtered value
  var cmax = Math.sqrt(2.0) * cu
 
  var alpha = ee.Image(1.0 + cu*cu).divide(ci.multiply(ci).subtract(cu*cu));
  var b = alpha.subtract(enl + 1.0)
  var d = mean.multiply(mean).multiply(b).multiply(b).add(alpha.multiply(mean).multiply(nat_img).multiply(4.0*enl));
  var f = b.multiply(mean).add(d.sqrt()).divide(alpha.multiply(2.0));
   
  var caster = ee.Dictionary.fromLists(bandNames,ee.List.repeat('float',3));
  var img1 = powerToDb(mean.updateMask(ci.lte(cu))).rename(bandNames).cast(caster);
  var img2 = powerToDb(f.updateMask(ci.gt(cu)).updateMask(ci.lt(cmax))).rename(bandNames).cast(caster);
  var img3 = img.updateMask(ci.gte(cmax)).rename(bandNames).cast(caster);
   
  // If ci > cmax do not filter at all (i.e. we don't do anything, other then masking)
  var result = ee.ImageCollection([img1,img2,img3])
    .reduce(ee.Reducer.firstNonNull()).rename(bandNames);
   
  // Compose a 3 band image with the mean filtered "pure speckle", the "low textured" filtered and the unfiltered portions
  return result;
}