Skip to content

Commit

Permalink
Merge pull request #9 from SDFIdk/fix_elevation_methods
Browse files Browse the repository at this point in the history 
fix some bugs in getWorldXYZ methods
  • Loading branch information
@iamfrank
iamfrank authored Oct 28, 2022
2 parents 085137a + 624f00f commit 62577a9
Show file tree
Hide file tree
Showing 3 changed files with 65 additions and 95 deletions.
81 changes: 3 additions & 78 deletions modules/saul-core.js
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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 {
Expand Down
22 changes: 20 additions & 2 deletions modules/saul-elevation.js
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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 ) )

Expand Down
57 changes: 42 additions & 15 deletions test/saul-elevation.test.js
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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'
Expand All @@ -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
Expand All @@ -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 => {
Expand All @@ -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 )

Expand All @@ -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)
}

})

Expand Down

0 comments on commit 62577a9

Please sign in to comment.