Wring mathematical formula for Sentinel-1 in Google Earth Engine

Question

I am trying to find permittivity using Sentinel-1 using following expression

It is inverse of Modified Dubois model for soil moisture. I am new to GEE so I tried to break the expression in to small components in following code

var admin1 = ee.FeatureCollection("FAO/GAUL_SIMPLIFIED_500m/2015/level1");
var Bihar = admin1.filter(ee.Filter.eq('ADM1_NAME', 'Bihar'))
var roi = Bihar.geometry()


//Define the time interval
var start_date = ee.Date('2021-11-15');
var end_date = start_date.advance(5,'days');
var date_filter = ee.Filter.date(start_date, end_date);

// Filter the collection for the VV product from the descending track
var collection1 = ee.ImageCollection('COPERNICUS/S1_GRD')
    .filter(ee.Filter.eq('instrumentMode', 'IW'))
    //.filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
    .filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING'))
    .filterMetadata ('resolution_meters', 'equals', 10)
    .filter(date_filter)
    .filterBounds(roi)
    // .select(['VV','angle']);
print(collection1);
var vis_paraVV = {bands:'VV', min:-25, max:5};
var vis_paraVH = {bands:'VH', min:-25, max:5};
var vis_paraAN = {bands:'angle', min:30.3, max:46.8};
Map.centerObject(roi, 6);
Map.addLayer(collection1,vis_paraVV,'VV');
Map.addLayer(collection1,vis_paraVH,'VH');
Map.addLayer(collection1,vis_paraAN,'Angle');

////Mosaic
var Mosaic1 = collection1.mosaic().clip(roi);
//var Mo_Angl1 = collection1.select(['angle']).mosaic().clip(roi);
// var VH_mosaic = collection1.select(['VH']).mosaic();
// var AN_mosaic = collection1.select(['angle']).mosaic();

//////Speckel Filter boxcar
var boxcar = ee.Kernel.square({ radius: 1.5, units: 'pixels', normalize: true});

function addfilter(image){
  var VH_fil = ((image.select('VH').convolve(boxcar))).rename('VH_fil');
  var VV_fil = ((image.select('VV').convolve(boxcar))).rename('VV_fil');
  return image.addBands(VH_fil).addBands(VV_fil);
}
////apply filter to VV and VH
var Mosaic1 = addfilter(Mosaic1);
//print(Mosaic1);
Map.addLayer(Mosaic1,{bands:'VV_fil', min:-25, max:5},'VVMo');
Map.addLayer(Mosaic1,{bands:'VH_fil', min:-25, max:5},'VHMo');
//Map.addLayer(Mo_VVVH1,vis_paraAN,'AngleMo');

///Mathemitical opretions
function addmath(image) {
  var VHminusVV = ((image.select('VH_fil')).subtract(image.select('VV_fil'))).rename('Diff');
  var VHbyVV = ((image.select('VH_fil')).divide(image.select('VV_fil'))).rename('Q');
  var DgtoRd = ((image.select('angle')).divide(180).multiply(Math.PI)).rename('Thita');
  return image.addBands(VHminusVV).addBands(VHbyVV).addBands(DgtoRd);
}

///Apply math
var Mosaic1 = addmath(Mosaic1);

//print(Mosaic1);
//Map.addLayer(Mosaic1) 

///Components for equations
function addcomp(image) {
  var A = ((image.select('Thita').tan())).rename('tan0');
  var B = ((image.select('Thita').sin())).rename('sin0');
  var C = ((image.select('Thita').cos())).rename('cos0');
  return image.addBands(A).addBands(B).addBands(C);
}

///Apply Components for equations
var Mosaic1 = addcomp(Mosaic1);
//print(Mosaic1);

///Components for complex 
var addsin150 = function(image){
  var Th = image.select('Thita');
  return image.addBands(image.expression('sin(1.5 * Th)',{Th:Th}).rename('sin150'));
};
////Apply
var Mosaic1 = addsin150(Mosaic1);
//print(Mosaic1);

////Surface rougness
var addSR = function(image) {
  var Diff = image.select(['Diff']);
  return image.addBands(image.expression ('4.27 + (0.22 * Diff)', {Diff:Diff}).rename('SR'));
};
var Mosaic1 = addSR(Mosaic1);
print(Mosaic1);
Map.addLayer(Mosaic1,{bands:'SR', min:2.0458423897394926, max:3.9287323550649575}, 'SurfaceRougness');

// ////Modified Dubois Model
// ////E= (1/0.046*tan0)*log10[((VV*10^2.35*wavelength^-0.7)/(k*s*sin0)^1.1)*(sin^3(0)/cos^3(0)]////

// ////functions for Components of Equation///
/// com1 = 1/0.046*tan0
/// com2 = VV*10^2.35*wavelenght^-0.7
/// com3 = (k*S*sin0)^1.1
/// com4 = sin0^3/cos0^3

var addcom1 = function(image){
   var T = image.select('tan0');
  return image.addBands(image.expression('(1 / (0.046 * T))',{T:T}).rename('com1'));
};

var addcom2 = function(image){
  var VV = image.select('VV_fil');
  return image.addBands(image.expression('(VV * (223.872113857) * (0.30142728421))',{VV:VV}).rename('com2')); 
};

var addcom3 = function(image){
  var SR = image.select('SR');
  var S = image.select('sin0');
  return image.addBands(image.expression('(1.13280427 * SR * S)**1.1',{SR:SR,S:S}).rename('com3'));
};

var addcom4 = function(image){
  var S = image.select('sin0');
  var C = image.select('cos0');
  return image.addBands(image.expression('(S**3) / (C**3)',{S:S,C:C}).rename('com4'));
};
var Mosaic1 = addcom1(addcom2(addcom3(addcom4(Mosaic1))));

////Calculations
////com5 = (VV*10^2.35*wavelength^-0.7)/(k*s*sin0)^1.1)
var addcom5 = function(image){
  var c2 = image.select('com2');
  var c3 = image.select('com3');
  return image.addBands(image.expression('c2 / c3',{c2:c2,c3:c3}).rename('com5'));
};
///Apply
var Mosaic1 = addcom5(Mosaic1);

////com6 = (VV*10^2.35*wavelength^-0.7)/(k*s*sin0)^1.1)*(sin^3(0)/cos^3(0)
function addcom6(image) {
  var com6 = ((image.select('com5').multiply(image.select('com4'))).rename('com6'));
  return image.addBands(com6);
}
///Apply 
var Mosaic1 = addcom6(Mosaic1);

////com7 = log10[((VV*10^2.35*wavelength^-0.7)/(k*s*sin0)^1.1)*(sin^3(0)/cos^3(0)]
var addcom7 = function (image){
  var c6 = image.select('com6');
  return image.addBands(image.expression('log(c6)',{c6:c6}).rename('com7'));
};
///Apply 
var Mosaic1 = addcom7(Mosaic1);
print (Mosaic1);
Map.addLayer(Mosaic1);

But I am getting zero value which is wrong. I want to know if my formula expression in GEE is written wrong and what can be done to improve it.

Answer 1

There is one bug in your expression for calculating the base-10 logarithm, that is log(c6) will give the natural logarithm, not the base-10 log. Rather, specify log10(c6). But this is not necessarily the issue to your output.

I ran the code setting the bands 'sin0','cos0','tan0','vv_fil', and 'SR' to be constant images of value 1. Using your code, the final output to this trial calculation gave the same ouput as my calculator did. This brings me to believe that the issue is not in the various expressions that you've defined to break-down the calculation, but rather the input data itself (so the processing of the Sentinel data).

To clarify, you are not getting zero values in your last step. Rather, the values are set to null values because you attempt to take the logarithm of a negative number (it seems like com6 < 0 for all pixels).

Answer 2

I believe the issue with obtaining a zero value is that you are inputting dB values, which are on the log scale and all negative, instead of linear values from Sentinel-1. Below is the code to convert from dB to linear scale.

///Convert dB to Linear scale
function toLinear(db) {
  return db.addBands(
    ee.Image().expression('pow(10, db / 10)', {
      db: db.select(['VV', 'VH'])
    }),
    null, true // Replace the bands to keep image properties
  )
}

var collection1 = collection1.map(toLinear)
Map.addLayer(collection1,{},'linear values')

https://code.earthengine.google.com/2810a7c24b822c71aa86fbf619c7b65f

See also: Convert Sentinel-1 images data from dB to linear

Edit: A couple other things for those that stumble on this thread. The 'log10' should not be in this equation. It should just be 'log'. This will result in the same Dielectric Constant as the formula found in: Srinivasa Rao, S., Dinesh kumar, S., Das, S.N., Nagaraju, M.S.S., Venugopal, M.V., Rajankar, P., Laghate, P., Reddy, M.S., Joshi, A.K., Sharma, J.R., 2013. Modified Dubois Model for Estimating Soil Moisture with Dual Polarized SAR Data. J Indian Soc Remote Sens 41, 865–872. https://doi.org/10.1007/s12524-013-0274-3.

In the script above, 'com1' must be multiplied by 'com7'.

Lastly, the Surface Roughness must be in dB, not linear values.