2

I need to calculate distances between two different pairs of pixels in a neighborhood.

enter image description here

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);
2
  • Can you post the example code of what you have tried so far?
    – Kersten
    Jun 30, 2018 at 6:43
  • I added the code and updated the problem. The rounding-to-zero problem is related only with longitude. Jun 30, 2018 at 18:04

1 Answer 1

2

The problem was solved by checking again the code. The haversine function had a minimal mistake, where the correct is

var a = p2.multiply(p3).add(p1); // a = p1 + p2 * p3

instead of

var a = p2.add(p3).multiply(p1); // a = p1 + p2 * p3

The full code, where it is possible to find the distance between pixels of a neighborhood for a Spherical Grid Projection is here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.