2

I asked Accessing slope and aspect data using Leaflet? last week about getting slope and aspect as a function of lat and lng in leaflet, and I am in the process of trying to code out a good answer. This article contains some great ideas about how to do this. After playing with it, I came up with this function:

   function getSlope(point){
   // get the xy position of the point within its DOM container:
   let xy = leafletMap.latLngToContainerPoint(point);
   // get the 4 pixels around it in xy coordinates:
   let rightxy = { x: xy.x + 1,  y: xy.y },
      leftxy = { x: xy.x - 1, y: xy.y },
      topxy = { x: xy.x, y: xy.y - 1 },
      bottomxy = { x: xy.x, y: xy.y + 1 };

   // Translate those 4 points back to lat lngs
   let right = leafletMap.containerPointToLatLng(rightxy),
      left = leafletMap.containerPointToLatLng(leftxy),
      top = leafletMap.containerPointToLatLng(topxy),
      bottom = leafletMap.containerPointToLatLng(bottomxy);

   // get the elevation of those points and get the difference
   let dzdx = getElevation(right) - getElevation(left);
   let dzdy = getElevation(top) - getElevation(bottom);

   // Get the angle of slope based on those differences
   let slopeAngle = Math.atan( Math.sqrt(dzdx ** 2 + dzdy ** 2)) * (180 / Math.PI);

   return slopeAngle;

}

You'll see the getElevation function in there, which pulls data from a mapbox rgb-encoded elevation tileset. The function looks like this:

function getElevation(location){
   var color = colorPicker.getColor(location)
   let R = color[0];
   let G = color[1];
   let B = color[2];

   let height = -10000 + ((R * 256 * 256 + G * 256 + B) * 0.1)
   return height;
}

The elevation function has been working consistently just fine. The issue is that I get wildly varying results for slope and aspect at different zoom levels. I understand that the zoom level will determine the color of the pixel being put into the elevation function, but the variance between, say, zoom level 16 and zoom level 10, seems huge. I wrote a more fleshed out function to be able to pick things apart. A little convoluted, but a good tester function:

function getSlopeTest(point){

   let slopeTest = {};
   // get the xy position of the point within its DOM container:
   slopeTest.xy = leafletMap.latLngToLayerPoint(point);
   // get the 4 pixels around it in xy coordinates:
   slopeTest.rightxy = { x: slopeTest.xy.x + 1,  y: slopeTest.xy.y };
   slopeTest.leftxy = { x: slopeTest.xy.x - 1, y: slopeTest.xy.y };
   slopeTest.topxy = { x: slopeTest.xy.x, y: slopeTest.xy.y - 1 };
   slopeTest.bottomxy = { x: slopeTest.xy.x, y: slopeTest.xy.y + 1 };

   // Translate those 4 points back to lat lngs
   slopeTest.right = leafletMap.layerPointToLatLng(slopeTest.rightxy);
   slopeTest.left = leafletMap.layerPointToLatLng(slopeTest.leftxy);
   slopeTest.top = leafletMap.layerPointToLatLng(slopeTest.topxy);
   slopeTest.bottom = leafletMap.layerPointToLatLng(slopeTest.bottomxy);

   // Get 4 points by destination function rathan than mapping to and from xy
   // slopeTest.right = L.GeometryUtil.destination(point, -90, 1);
   // slopeTest.left = L.GeometryUtil.destination(point, 90, 1);
   // slopeTest.top = L.GeometryUtil.destination(point, 0, 1);
   // slopeTest.bottom = L.GeometryUtil.destination(point, 180, 1);

   // Get the elevations of those 4 points
   slopeTest.rightElevation = getElevation( slopeTest.right );
   slopeTest.leftElevation = getElevation( slopeTest.left );
   slopeTest.topElevation = getElevation( slopeTest.top );
   slopeTest.bottomElevation = getElevation( slopeTest.bottom );

   // Get the distance between the points
   slopeTest.dx = L.GeometryUtil.distance(leafletMap, slopeTest.right, slopeTest.left)
   slopeTest.dy = L.GeometryUtil.distance(leafletMap, slopeTest.top, slopeTest.bottom)

   // get the elevation of those points and get the difference
   slopeTest.dzdx = ( slopeTest.rightElevation - slopeTest.leftElevation) / slopeTest.dx;
   slopeTest.dzdy = ( slopeTest.topElevation - slopeTest.bottomElevation ) / slopeTest.dy;

   // Get the angle of slope based on those differences
   slopeTest.slopeAngle = Math.atan( Math.sqrt(slopeTest.dzdx ** 2 + slopeTest.dzdy ** 2)) * (180 / Math.PI);
   slopeTest.slopeAspect = Math.atan2( slopeTest.dzdy, slopeTest.dzdx ) * (180 / Math.PI);


   return slopeTest;

}

The strangest thing to me is that no matter what zoom level I'm at, the slopeTest.dx and .dy are always the same. This is supposed to be a comparison of distance between two pixels, whose number should be higher when the map is more zoomed out. But it seems to always give me a result of 2.

Here is the codepen to see whats happening. I included some buttons to drop a marker whose popup contains the slope and aspect values at that point. If you zoom in and out (leaving the map center alone) and drop a marker at various zoom levels, you'll see that the slope and aspect results are very different between zoom levels, especially as you zoom further out. You can try the "Create Marker with Info Grid" to see a layout of some of the values.

Realistically I don't need a reliable measure of slope and aspect at zoom levels under 10. I can tolerate some variance in my slope and aspect between zoom levels over 10, but I'm getting unexpectedly large differences between zoom levels to the point where I'm not sure if my calculation are even correct anymore.

If anyone can spot the issues here can you let me know?

  • Since you are using pixel dimension as offset when calculation slope, you should expect different result when zooming in/out. Pixel dimension on map changes with zoom change and so also point coordinates you are using to calculate slope. – TomazicM Oct 28 at 15:54
  • Right. I said that my .dx and .dy were always the same regardless of zoom. Apparently L.GeometryUtil.distance is not what I thought. I switched it to map.distance(), a native leaflet function, which gives the correct distance between points, in meters. This gives different values at different zoom levels, as expected. Now my values are a little less outrageous. For example, for zoom levels 18 to 10, my slope and aspect values still vary from 29 to 51 degrees, and 0 to 75, respectively. Now that I've correct for the dx and dy issue, what else might be causing this inconsistency? – Seth Lutske Oct 28 at 18:42
  • Now it's problem of accuracy. When you zoom out, each DEM pixel at given zoom is an average of four pixels at prevous higher zoom, so you are bound to get different result. Best example would be a very sharp image at 1024 x 1024 pixels, resized to 256 x 256 pixels. It will loose a great deal of sharpness/accuracy. – TomazicM Oct 28 at 18:59
  • Would it be possible to instead access a different zoom-level tile and pull the DEM values from that? For example, even if the map is being viewed at zoom level 8, can a function getSlope(point) somehow load the tile at that point for a given zoom level (say zoom 16), then calculate values based on that, without changing the user view? This is all I can think of to try to get more accuracy and answer my initial question. – Seth Lutske Oct 28 at 19:14
  • Yes, that would be the solution. It can be done, but you would have to write your own function to get the desired zoom tile for lat,lng you want to get accurate height. – TomazicM Oct 28 at 19:35

Your Answer

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

Browse other questions tagged or ask your own question.