Skip to content

Commit

Permalink
Merge pull request #10 from SDFIdk/elevation_station
Browse files Browse the repository at this point in the history
Add elevation tests
  • Loading branch information
iamfrank authored Nov 28, 2022
2 parents 35d3b33 + dcf120c commit e7ba3b4
Show file tree
Hide file tree
Showing 10 changed files with 228 additions and 57 deletions.
2 changes: 1 addition & 1 deletion config.js.example
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@ export default {

// STAC TOKEN can be aquired from https://dataforsyningen.dk/
API_STAC_TOKEN: "",
API_STAC_BASEURL: "https://api.dataforsyningen.dk/skraafotoapi",
API_STAC_BASEURL: "https://api.dataforsyningen.dk/skraafoto_api_test/v1.0",

// DHM URL and auth can be aquired from https://datafordeler.dk/dataoversigt/danmarks-hoejdemodel-dhm/koter/
API_DHM_WCS_BASEURL: "https://services.datafordeler.dk/DHMNedboer/dhm_wcs/1.0.0/WCS",
Expand Down
5 changes: 3 additions & 2 deletions index.js
Original file line number Diff line number Diff line change
Expand Up @@ -2,8 +2,9 @@

export {
image2world,
world2image,
world2image, // Deprecated. Use getImageXY
getWorldXYZ,
getImageXY,
getZ, // Deprecated. Use getElevation
iterate // Deprecated. Use getWorldXYZ
} from './modules/saul-core.js'
Expand All @@ -12,6 +13,7 @@ export {
postSTAC,
getSTAC,
getDHM,
getTerrainGeoTIFF,
get,
post
} from './modules/api.js'
Expand All @@ -22,6 +24,5 @@ export {
} from './modules/saul-projection.js'

export {
getTerrainGeoTIFF,
getElevation
} from './modules/saul-elevation.js'
44 changes: 43 additions & 1 deletion modules/api.js
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,8 @@
* SAUL API utilities
*/

import { fromArrayBuffer } from 'geotiff'

let error_msg
let load_stack = []

Expand Down Expand Up @@ -184,10 +186,50 @@ function postSTAC(endpoint, data, auth) {
.then((data) => data)
}


/** Converts raw GeoTIFF arrayBuffer to image */
async function consumeGeoTIFF(raw_data) {
const tiff = await fromArrayBuffer(raw_data)
const image = await tiff.getImage()
return image
}


/** Fetches a GeoTIFF with elevation data matching the bounding box of a STAC API item (image)
* @param {object} stac_item - STAC API item from a featureCollection request
* @param {{API_DHM_WCS_BASEURL: string, API_DHM_USERNAME: string, API_DHM_PASSWORD: string}} auth - API autentication data. See ../config.js.example for reference.
* @param {number} [resolution] - Resolution (1 - 0.01). Higher number means more pixels and better precision.
* @returns {object} GeoTiff data
*/
function getTerrainGeoTIFF(stac_item, auth, resolution = 0.05) {

const bbox = stac_item.bbox
const width = Math.round( stac_item.properties['proj:shape'][0] * resolution )
const height = Math.round( stac_item.properties['proj:shape'][1] * resolution )

// GET request for DHM WCS data
let url = auth.API_DHM_WCS_BASEURL
url += '?SERVICE=WCS&COVERAGE=dhm_terraen&RESPONSE_CRS=epsg:25832&CRS=epsg:25832&FORMAT=GTiff&REQUEST=GetCoverage&VERSION=1.0.0'
url += `&username=${ auth.API_DHM_USERNAME }&password=${ auth.API_DHM_PASSWORD }` // TODO: Should be auth stuff
url += `&height=${ height }`
url += `&width=${ width }`
url += `&bbox=${ bbox[0]},${ bbox[1]},${ bbox[2]},${ bbox[3]}`

return get(url, {}, false)
.then((response) => {
return response.arrayBuffer()
})
.then((arrayBuffer) => {
return consumeGeoTIFF(arrayBuffer)
})
}


export {
get,
post,
getSTAC,
postSTAC,
getDHM
getDHM,
getTerrainGeoTIFF
}
30 changes: 25 additions & 5 deletions modules/saul-core.js
Original file line number Diff line number Diff line change
Expand Up @@ -57,18 +57,22 @@ function image2world(image_data, col, row, Z = 0) {
var X = (Z-Z0)*kx + X0
var Y = (Z-Z0)*ky + Y0

return[X,Y,Z]
return[
Math.round(X * 100) / 100,
Math.round(Y * 100) / 100,
Z
]
}

/**
* Converts lat,lon coordinates to x,y coordinates within a specific image
* Converts world lat,lon coordinates to x,y coordinates within a specific image
* @param {Object} image_data - skraafoto-stac-api image data
* @param {Number} Y - northing
* @param {Number} X - easting
* @param {Number} [Z] - elevation (geoide)
* @returns {array} [x,y] Column/row image coordinate
*/
function world2image(image_data, X, Y, Z = 0) {
function getImageXY(image_data, X, Y, Z = 0) {

// constants pulled from image_data
const xx0 = image_data.properties['pers:interior_orientation'].principal_point_offset[0]
Expand Down Expand Up @@ -109,7 +113,22 @@ function world2image(image_data, X, Y, Z = 0) {
var col = ((x_dot-xx0)+(dimX))*(-1)/pix
var row = ((y_dot-yy0)+(dimY))*(-1)/pix

return [col, row]
return [
Math.round(col),
Math.round(row)
]
}

/**
* (DEPRECATED - Use getImageXY) Converts lat,lon coordinates to x,y coordinates within a specific image
* @param {Object} image_data - skraafoto-stac-api image data
* @param {Number} Y - northing
* @param {Number} X - easting
* @param {Number} [Z] - elevation (geoide)
* @returns {array} [x,y] Column/row image coordinate
*/
function world2image(image_data, X, Y, Z = 0) {
return getImageXY(image_data, X, Y, Z)
}

/** Converts degress to radians */
Expand Down Expand Up @@ -216,8 +235,9 @@ async function getWorldXYZ(options, precision = 0.3) {

export {
image2world,
getWorldXYZ,
world2image,
getWorldXYZ,
getImageXY,
getZ,
iterate
}
Expand Down
38 changes: 1 addition & 37 deletions modules/saul-elevation.js
Original file line number Diff line number Diff line change
@@ -1,41 +1,5 @@
import { fromArrayBuffer } from 'geotiff'
import { get } from './api.js'
import { getTerrainGeoTIFF } from './api.js'

/** Fetches a GeoTIFF with elevation data matching the bounding box of a STAC API item (image)
* @param {object} stac_item - STAC API item from a featureCollection request
* @param {{API_DHM_WCS_BASEURL: string, API_DHM_USERNAME: string, API_DHM_PASSWORD: string}} auth - API autentication data. See ../config.js.example for reference.
* @param {number} [fidelity] - Resolution fidelity. Higher number means more pixels and better precision. Between 1 and 0.01.
* @returns {object} GeoTiff data
*/
function getTerrainGeoTIFF(stac_item, auth, fidelity = 0.05) {

const bbox = stac_item.bbox
const width = Math.round( stac_item.properties['proj:shape'][0] * fidelity )
const height = Math.round( stac_item.properties['proj:shape'][1] * fidelity )

// GET request for DHM WCS data
let url = auth.API_DHM_WCS_BASEURL
url += '?SERVICE=WCS&COVERAGE=dhm_terraen&RESPONSE_CRS=epsg:25832&CRS=epsg:25832&FORMAT=GTiff&REQUEST=GetCoverage&VERSION=1.0.0'
url += `&username=${ auth.API_DHM_USERNAME }&password=${ auth.API_DHM_PASSWORD }` // TODO: Should be auth stuff
url += `&height=${ height }`
url += `&width=${ width }`
url += `&bbox=${ bbox[0]},${ bbox[1]},${ bbox[2]},${ bbox[3]}`

return get(url, {}, false)
.then((response) => {
return response.arrayBuffer()
})
.then((arrayBuffer) => {
return consumeGeoTIFF(arrayBuffer)
})
}

/** Converts raw GeoTIFF arrayBuffer to image */
async function consumeGeoTIFF(raw_data) {
const tiff = await fromArrayBuffer(raw_data)
const image = await tiff.getImage()
return image
}

/** Constrain coordinates to be within a bounding box */
function constrainToBbox(bbox, x, y) {
Expand Down
2 changes: 1 addition & 1 deletion package.json
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@
"main": "index.js",
"type": "module",
"scripts": {
"test": "echo 'TESTS WILL ONLY RUN UNDER NODE v18' && node ./test/saul-projection.test.js && node ./test/saul-core.test.js && node ./test/saul-elevation.test.js",
"test": "echo 'TESTS WILL ONLY RUN UNDER NODE v18' && node ./test/saul-projection.test.js && node ./test/saul-core.test.js && node ./test/saul-combined.test.js && node ./test/saul-elevation.test.js",
"docs": "jsdoc ./modules/*.js --destination docs"
},
"repository": {
Expand Down
148 changes: 148 additions & 0 deletions test/saul-combined.test.js
Original file line number Diff line number Diff line change
@@ -0,0 +1,148 @@
// You should add your own config.js with proper tokens, passwords, etc.
// Check config.js.example for info on how to set it up
import auth from '../config.js'
import assert from 'assert'
import { getSTAC, getTerrainGeoTIFF,getZ, getWorldXYZ, getImageXY } from '../index.js'


// Vars
const resolution = 0.005 // Higher number means more points and better precision
const max_deviation = 0.5

/*
// North
let item_1 = {
id: '2021_83_29_2_0019_00003995',
xy: [1000,1000]
}
// East
let item_2 = {
id: '2021_83_29_4_0016_00002877'
}
// West
let item_3 = {
id: '2021_83_29_5_0022_00003276'
}
// South
let item_4 = {
id: '2021_83_29_3_0019_00003975'
}
// Nadir
let item_5 = {
id: '2021_83_29_1_0019_00003985'
}
*/

/*
// North
let item_1 = {
id: '2021_83_33_2_0003_00002821',
xy: [8000,5000]
}
// East
let item_2 = {
id: '2021_83_33_4_0002_00000845'
}
// West
let item_3 = {
id: '2021_83_33_5_0005_00003192'
}
// South
let item_4 = {
id: '2021_83_33_3_0004_00002970'
}
// Nadir
let item_5 = {
id: '2021_83_33_1_0004_00002965'
}
*/

// North
let item_1 = {
id: '2021_83_29_2_0019_00003995',
xy: [1000,1000]
}
// East
let item_2 = {
id: '2021_83_29_4_0016_00002877'
}
// West
let item_3 = {
id: '2021_83_29_5_0022_00003276'
}
// South
let item_4 = {
id: '2021_83_29_3_0019_00003975'
}
// Nadir
let item_5 = {
id: '2021_83_29_1_0019_00003985'
}

// Helper functions
function is_equalIsh(num1, num2, deviation = max_deviation) {
if (Math.abs(num1 - num2) > deviation) {
return false
} else {
return true
}
}
function generateSTACurl(stac_item_id) {
return `/search?limit=1&crs=http://www.opengis.net/def/crs/EPSG/0/25832&ids=${ stac_item_id }`
}
async function enrichData(item) {
let new_item = Object.assign({}, item)
// add STAC item
new_item.item = await getSTAC(generateSTACurl(new_item.id), auth)
new_item.item = new_item.item.features[0]
// add terrain GeoTIFF
new_item.terrain = await getTerrainGeoTIFF(new_item.item, auth, resolution)
// Add imageXY
new_item.xy = getImageXY(new_item.item, item_1.xyz[0], item_1.xyz[1], item_1.xyz[2])
// Add world XYZ (from imageYX)
new_item.xyz = await getWorldXYZ({
xy: new_item.xy,
image: new_item.item,
terrain: new_item.terrain
})
// add kote
new_item.kote = await getZ(new_item.xyz[0], new_item.xyz[1], auth)
// add resolution info
new_item.ext_x = Math.round((new_item.item.bbox[2] - new_item.item.bbox[0]) * 100) / 100
new_item.ext_y = Math.round((new_item.item.bbox[3] - new_item.item.bbox[1]) * 100) / 100
new_item.z_pr_m_x = Math.round((new_item.item.properties['proj:shape'][0] / new_item.ext_x) * 100 * resolution) / 100
new_item.z_pr_m_y = Math.round((new_item.item.properties['proj:shape'][1] / new_item.ext_y) * 100 * resolution) / 100
const width = Math.round( new_item.item.properties['proj:shape'][0] * resolution )
const height = Math.round( new_item.item.properties['proj:shape'][1] * resolution )
new_item.z_points = width * height
return new_item
}

// Fetch image data and terrain
item_1.item = await getSTAC(generateSTACurl(item_1.id), auth)
item_1.item = item_1.item.features[0]
item_1.terrain = await getTerrainGeoTIFF(item_1.item, auth, resolution)

// Given image XY for one image, get world XYZ
item_1.xyz = await getWorldXYZ({
xy: item_1.xy,
image: item_1.item,
terrain: item_1.terrain
})

item_1.kote = await getZ(item_1.xyz[0], item_1.xyz[1], auth)

// Get image and corresponding world coordinates for other images (amongst other things)
item_2 = await enrichData(item_2)
item_3 = await enrichData(item_3)
item_4 = await enrichData(item_4)
item_5 = await enrichData(item_5)

// Compare items
console.table({
north: item_1,
east: item_2,
west: item_3,
south: item_4,
nadir: item_5
})
8 changes: 3 additions & 5 deletions test/saul-core.test.js
Original file line number Diff line number Diff line change
Expand Up @@ -5,8 +5,7 @@
import auth from '../config.js'

import assert from 'assert'
import {getSTAC} from '../modules/api.js'
import {world2image, image2world, getZ, iterate} from '../modules/saul-core.js'
import {getSTAC, world2image, image2world, getZ, iterate} from '../index.js'

/*
const item = await getSTAC('/collections/skraafotos2019/items/2019_83_37_2_0046_00001113', auth)
Expand All @@ -24,8 +23,7 @@ const image_x = 41839.698939
const image_y = 3204.326930
const world_elevation = 38.874336

function is_equalIsh(num1, num2) {
const deviation = 0.05
function is_equalIsh(num1, num2, deviation = 0.05) {
if (Math.abs(num1 - num2) > deviation) {
return false
} else {
Expand Down Expand Up @@ -65,7 +63,7 @@ try {

let xy = world2image(item, world_x, world_y, world_elevation)

assert(is_equalIsh(xy[0], image_x), "world2image fail: The x coordinate values are not equal")
assert(is_equalIsh(xy[0], image_x, 1), `world2image fail: The x coordinate values ${xy[0]} ${image_x} are not equal`)
assert(is_equalIsh(xy[1], image_y), "world2image fail: The y coordinate values are not equal")

console.log("Test world2image OK")
Expand Down
Loading

0 comments on commit e7ba3b4

Please sign in to comment.