I need to calculate distances between two different pairs of pixels in a neighborhood.
There are two solutions. The first is to reproject the image to a plane squared grid (Albers Equal Area or UTM, for example) in its original scale and get a proportional pixel size in the whole area. However, I've noticed that reprojecting is time-consuming and is not suitable for large areas with small pixel size.
The second solution which I'm developing is to calculate the distances between to points in the original Mercator projection (which is a Spheroidal Equal Angular Grid), using haversine formulae and ee.Image.pixelLonLat().
When I calculate the distance between two positions with different Latitude, the algorithm works well. When calculating the distances for points with different Longitude, the result is zero. I checked the Longitude position within the neighborhood and curiously they are the same.
How could I overcome this? The code is below or from here
var retenv = ee.Geometry.Rectangle({coords: [-48,-23,-47,-22],
geodesic: false});
// Defining DEM
var demSRTM = ee.Image('USGS/SRTMGL1_003').clip(retenv);
// Defining positions
var demPositions = demSRTM.addBands(ee.Image.pixelLonLat()).clip(retenv);
print(demPositions, 'DEM with positions');
//Map.addLayer(demPositions, {bands: 'longitude', min: -48, max: -46.5}, "Longitude");
// Weights for a 3x3 kernel
var w00 = [0, 0, 0];
var w11 = [1, 0, 0];
var w12 = [0, 1, 0];
var w13 = [0, 0, 1];
// Neighborhood indices
var n1 = [w11, w00, w00];
var n2 = [w12, w00, w00];
var n3 = [w13, w00, w00];
var n4 = [w00, w11, w00];
var n5 = [w00, w12, w00];
var n6 = [w00, w13, w00];
var n7 = [w00, w00, w11];
var n8 = [w00, w00, w12];
var n9 = [w00, w00, w13];
// Kernel for each neighborhood index
var kerneln1 = ee.Kernel.fixed(3, 3, n1, 1, 1, false);
var kerneln2 = ee.Kernel.fixed(3, 3, n2, 1, 1, false);
var kerneln3 = ee.Kernel.fixed(3, 3, n3, 1, 1, false);
var kerneln4 = ee.Kernel.fixed(3, 3, n4, 1, 1, false);
var kerneln5 = ee.Kernel.fixed(3, 3, n5, 1, 1, false);
var kerneln6 = ee.Kernel.fixed(3, 3, n6, 1, 1, false);
var kerneln7 = ee.Kernel.fixed(3, 3, n7, 1, 1, false);
var kerneln8 = ee.Kernel.fixed(3, 3, n8, 1, 1, false);
var kerneln9 = ee.Kernel.fixed(3, 3, n9, 1, 1, false);
// Function to compute single neighborhood values
var addNParameters = function(image) {
var N1 = image.convolve(kerneln1);
var N2 = image.convolve(kerneln2);
var N3 = image.convolve(kerneln3);
var N4 = image.convolve(kerneln4);
var N5 = image.convolve(kerneln5);
var N6 = image.convolve(kerneln6);
var N7 = image.convolve(kerneln7);
var N8 = image.convolve(kerneln8);
var N9 = image.convolve(kerneln9);
return image.addBands([N1,N2,N3,N4,N5,N6,N7,N8,N9]);
};
// Adding the neighborhood values
// Longitude
var long = demPositions.select('longitude');
var longNParameters = addNParameters(long);
var longNParameters = longNParameters.rename(['longitude','longN1','longN2','longN3','longN4','longN5','longN6','longN7','longN8','longN9']);
print(longNParameters, 'Longitude with neighborhood parameters');
// Latitude
var lat = demPositions.select('latitude');
var latNParameters = addNParameters(lat);
var latNParameters = latNParameters.rename(['latitude','latN1','latN2','latN3','latN4','latN5','latN6','latN7','latN8','latN9']);
print(latNParameters, 'Latitude with neighborhood parameters');
// Function with haversine formula to retrieve distances between two points
var haversineFunction = function(imageLat, imageLong, latNeigh1, latNeigh2, longNeigh1, longNeigh2) {
var φ1 = imageLat.select(ee.String(latNeigh1)).divide(180).multiply(Math.PI); // to Radians
var φ2 = imageLat.select(ee.String(latNeigh2)).divide(180).multiply(Math.PI); // to Radians
var λ1 = imageLong.select(ee.String(longNeigh1)).divide(180).multiply(Math.PI); // to Radians
var λ2 = imageLong.select(ee.String(longNeigh2)).divide(180).multiply(Math.PI); // to Radians
var Δφ = φ2.subtract(φ1); // (φ2 - φ1) // here is where I find zero values for long distances
var Δλ = λ2.subtract(λ1); // (λ2 - λ1) // here is where I find zero values for long distances
var p1 = Δφ.divide(2).sin().multiply(Δφ.divide(2).sin()); // sin(Δφ/2) * sin(Δφ/2)
var p2 = φ1.cos().multiply(φ2.cos()); // cos(φ1) * cos(φ2)
var p3 = Δλ.divide(2).sin().multiply(Δλ.divide(2).sin()); // sin(Δλ/2) * sin(Δλ/2)
var a = p2.add(p3).multiply(p1); // a = p1 + p2 * p3
var p4 = ee.Image(ee.Number(1)).clip(retenv); // p4 = image with constant 1
var p5 = ee.Image(ee.Number(2)).clip(retenv); // p5 = image with constant 2
var p6 = p4.subtract(a).sqrt(); // sqrt(1-a)
var p7 = a.sqrt(); // sqrt(a)
var p8 = p6.atan2(p7); // atan2(p6,p7)
var c = p5.multiply(p8); // c = 2 * p8
var R = ee.Image(ee.Number(6371000)).clip(retenv); // radius of Earth
var d = R.multiply(c); // d = R * c which is the distance between two points
return d;
};
var lenghtOfE = haversineFunction(latNParameters, longNParameters, 'latN1', 'latN4', 'longN1', 'longN4');
var lenghtOfD = haversineFunction(latNParameters, longNParameters, 'latN4', 'latN7', 'longN4', 'longN7');
var lenghtOfC = haversineFunction(latNParameters, longNParameters, 'latN1', 'latN2', 'longN1', 'longN2');
var lenghtOfB = haversineFunction(latNParameters, longNParameters, 'latN4', 'latN5', 'longN4', 'longN5');
var lenghtOfA = haversineFunction(latNParameters, longNParameters, 'latN7', 'latN8', 'longN7', 'longN8');
Map.addLayer(lenghtOfE, {}, 'Lenght of e');
Map.addLayer(lenghtOfD, {}, 'Lenght of d');
Map.addLayer(lenghtOfC, {}, 'Lenght of c');
Map.addLayer(lenghtOfB, {}, 'Lenght of b');
Map.addLayer(lenghtOfA, {}, 'Lenght of a');
Map.centerObject(retenv, 9);
