diff --git a/modules/saul-core.js b/modules/saul-core.js index 7994bea..4b067c4 100644 --- a/modules/saul-core.js +++ b/modules/saul-core.js @@ -5,8 +5,6 @@ import { getDHM } from './api.js' import { getElevation } from './saul-elevation.js' -let iterations - /** * Converts x,y coordinates from an image to real world lat,lon coordinates * @param {object} image_data - skraafoto-stac-api image data @@ -184,84 +182,11 @@ function iterate(image_data, col, row, auth, limit = 1) { */ async function getWorldXYZ(options) { - iterations = 1 - const world_xyz = image2world(options.image, options.xy[0], options.xy[1], 0) - const better_world_xyz = await compareXYZ({ - xy: options.xy, - image: options.image, - terrain: options.terrain, - world_xyz: Array.from(world_xyz), - limit: 0.5 - }) - return better_world_xyz -} - -/** - * Adjusts elevation value if values are too far apart - * @param {number} z - elevation - * @param {number} delta_y - discrepancy value - * @returns {number} Adjusted elevation number - */ -function adjustZ(z, delta_y) { - let new_z - if (delta_y < 0) { - new_z = z + 0.1 - } else { - new_z = z - 0.1 - } - console.log('z', new_z) - return new_z -} - -/** - * Compares differences i x and y values - * @param {number} delta_x - x difference - * @param {number} delta_y - y difference - * @param {number} limit - difference may be within this limit - * @returns {boolean} Returns true if differences are within accepted limit - */ -function compareXY(delta_x, delta_y, limit) { - if (Math.abs(delta_x) > limit || Math.abs(delta_y) > limit) { - return false - } else { - return true - } -} - -/** - * Iterates world coordinate calculations in order to get a more precise result - * @param {object} options - * @param {array} options.xy - image coordinates - * @param {object} options.image - STAC API item - * @param {object} options.terrain - GeoTIFF object from getTerrainGeoTIFF() output - * @param {array} options.world_xy - Array with (currently guessed) world coordinates (indexes 0 and 1) and elevation (index 2). Coordinates are EPSG:25832 format. - * @param {number} options.limit - Return when two guesses are within this limit - * @returns {array} Array with world coordinates (indexes 0 and 1) and elevation (index 2). Coordinates are EPSG:25832 format. - */ -async function compareXYZ(options) { - - iterations++ - if (iterations > 50) { - console.error('too many iterations') - return false - } - - const z = await getElevation(options.xy[0], options.xy[1], options.terrain) - const world_xyz = image2world(options.image, options.xy[0], options.xy[1], z) - const delta_x = world_xyz[0] - options.world_xyz[0] - const delta_y = world_xyz[1] - options.world_xyz[1] + world_xyz[2] = await getElevation(world_xyz[0], world_xyz[1], options.terrain) + const better_world_xyz = image2world(options.image, options.xy[0], options.xy[1], world_xyz[2]) - if (!compareXY(delta_x, delta_y, options.limit)) { - - const new_options = Object.assign({},options) - new_options.world_xyz = Array.from(world_xyz) - new_options.world_xyz[2] = adjustZ(world_xyz[2], delta_x, delta_y) - return await compareXYZ(new_options) - - } else { - return world_xyz - } + return better_world_xyz } export { diff --git a/modules/saul-elevation.js b/modules/saul-elevation.js index e11d790..0e7d2d5 100644 --- a/modules/saul-elevation.js +++ b/modules/saul-elevation.js @@ -36,6 +36,22 @@ async function consumeGeoTIFF(raw_data) { return image } +function constrainToBbox(bbox, x, y) { + let new_x = x + let new_y = y + if (x < bbox[0]) { + new_x = bbox[0] + 1 + } else if (bbox[2] < x) { + new_x = bbox[2] - 1 + } + if (y < bbox[1]) { + new_y = bbox[1] + 1 + } else if (bbox[3] < y) { + new_y = bbox[3] - 1 + } + return [new_x, new_y] +} + /** Calculates nearest elevation for given coordinate using geoTIFF terrain data * Heavily inspired by https://towardsdatascience.com/geotiff-coordinate-querying-with-javascript-5e6caaaf88cf * @param {number} x - EPSG:25832 X coordinate @@ -51,8 +67,10 @@ async function getElevation(x, y, geoTiff) { const bboxWidth = bbox[2] - bbox[0] const bboxHeight = bbox[3] - bbox[1] - const widthPct = ( x - bbox[0] ) / bboxWidth - const heightPct = ( y - bbox[1] ) / bboxHeight + const xy = constrainToBbox(bbox, x, y) + + const widthPct = ( xy[0] - bbox[0] ) / bboxWidth + const heightPct = ( xy[1] - bbox[1] ) / bboxHeight const xPx = Math.floor( width * widthPct ) const yPx = Math.floor( height * ( 1 - heightPct ) ) diff --git a/test/saul-elevation.test.js b/test/saul-elevation.test.js index b464206..0d66a0c 100644 --- a/test/saul-elevation.test.js +++ b/test/saul-elevation.test.js @@ -5,7 +5,7 @@ import auth from '../config.js' import assert from 'assert' import { getTerrainGeoTIFF, getElevation } from '../modules/saul-elevation.js' import { get } from '../modules/api.js' -import { getZ, getWorldXYZ } from '../modules/saul-core.js' +import { getZ, getWorldXYZ, world2image } from '../modules/saul-core.js' // Vars const stac_item = '2021_83_29_2_0019_00003995' @@ -16,8 +16,8 @@ const max_deviation = 0.5 let url_stac = 'https://api.dataforsyningen.dk/skraafotoapi_test/search?limit=1&crs=http://www.opengis.net/def/crs/EPSG/0/25832&token=9b554b6c854184c3b0f377ffc7481585' url_stac += `&ids=${ stac_item }` -function is_equalIsh(num1, num2) { - if (Math.abs(num1 - num2) > max_deviation) { +function is_equalIsh(num1, num2, deviation = max_deviation) { + if (Math.abs(num1 - num2) > deviation) { return false } else { return true @@ -30,13 +30,38 @@ function getRandomCoordinate(bbox) { return [x,y] } -function testDataAnumberOfTimes(data, bbox, times) { +function testGetElevationAnumberOfTimes(data, bbox, times) { for (let i = 1; i <= times; i++) { const world1 = getRandomCoordinate(bbox) compareElevations(world1[0], world1[1], data) } } +function getRandomImageXY(image) { + const shapex = image.properties['proj:shape'][0] + const shapey = image.properties['proj:shape'][1] + const x = Math.floor(shapex - Math.random() * shapex) + const y = Math.floor(shapey - Math.random() * shapey) + return [x,y] +} + +function testGetWorldXYZAnumberOfTimes(item, terrain, times) { + for (let i = 1; i <= times; i++) { + const xy = getRandomImageXY(item) + getWorldXYZ({ + xy: xy, + image: item, + terrain: terrain + }).then(world_xy => { + compareElevations(world_xy[0], world_xy[1], terrain) // lower left corner of image + const image_coords = world2image(item, world_xy[0], world_xy[1], world_xy[2]) + assert(is_equalIsh(image_coords[0], xy[0], 0.75), `Image x coordinates ${image_coords[0]} / ${xy[0]} do not match`) + assert(is_equalIsh(image_coords[1], xy[1], 0.75), `Image y coordinates ${image_coords[1]} / ${xy[1]} do not match`) + console.log('getWorldXYZ => world2image OK') + }) + } +} + function compareElevations(x,y,geotiff) { getElevation(x, y, geotiff) .then(elevation => { @@ -52,8 +77,6 @@ function compareElevations(x,y,geotiff) { get(url_stac) .then((json) => { - console.log('Testing getTerrainGeoTIFF and getElevation') - const width = Math.round( json.features[0].properties['proj:shape'][0] * fidelity ) const height = Math.round( json.features[0].properties['proj:shape'][1] * fidelity ) @@ -62,15 +85,19 @@ get(url_stac) getTerrainGeoTIFF(json.features[0], auth, fidelity) .then(data => { - testDataAnumberOfTimes(data, json.features[0].bbox, 5) - - getWorldXYZ({ - xy: [1,1], - image: json.features[0], - terrain: data - }).then(world_xy => { - compareElevations(world_xy[0], world_xy[1], data) // lower left corner of image - }) + // Testing getElevation + try { + testGetElevationAnumberOfTimes(data, json.features[0].bbox, 5) + } catch(error) { + console.error(error) + } + + // Testing getWorldXYZ + try { + testGetWorldXYZAnumberOfTimes(json.features[0], data, 10) + } catch(error) { + console.error(error) + } })