Skip to content

Commit

Permalink
Merge pull request #14 from SDFIdk/update_gtiff_elevation_methods
Browse files Browse the repository at this point in the history
Update gtiff elevation methods
  • Loading branch information
iamfrank authored Sep 11, 2024
2 parents 609d861 + c764b74 commit 582b38b
Show file tree
Hide file tree
Showing 12 changed files with 566 additions and 116 deletions.
Binary file added assets/dk-terrain.tiff
Binary file not shown.
Binary file added assets/tiny-dk-terrain.tiff
Binary file not shown.
4 changes: 3 additions & 1 deletion index.js
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@ export {
getSTAC,
getDHM,
getTerrainGeoTIFF,
getDenmarkGeoTiff,
get,
post
} from './modules/api.js'
Expand All @@ -24,5 +25,6 @@ export {
} from './modules/saul-projection.js'

export {
getElevation
getElevation,
visualizeGeotiff
} from './modules/saul-elevation.js'
37 changes: 32 additions & 5 deletions modules/api.js
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
/**
* SAUL API utilities
* SAUL API utilities
*/

import { fromArrayBuffer } from 'geotiff'
Expand Down Expand Up @@ -202,18 +202,24 @@ async function consumeGeoTIFF(raw_data) {
return image
}

function calcSizeRatio(sizeX, bbox) {
const ratio = ( bbox[0] - bbox[2] ) / ( bbox[1] - bbox[3] )
return Math.round(sizeX * ratio)
}


/** 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_TOKENA: string, API_DHM_TOKENB: 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) {
function getTerrainGeoTIFF(stac_item, auth, resolution = 0.05, sizeX = 300) {

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 )
const sizeY = calcSizeRatio(sizeX, stac_item.bbox)
const width = stac_item.properties ? Math.round( stac_item.properties['proj:shape'][0] * resolution ) : sizeX
const height = stac_item.properties ? Math.round( stac_item.properties['proj:shape'][1] * resolution ) : sizeY

// GET request for DHM WCS data
let url = auth.API_DHM_WCS_BASEURL
Expand All @@ -232,12 +238,33 @@ function getTerrainGeoTIFF(stac_item, auth, resolution = 0.05) {
})
}

/**
* Fetches a geoTIFF with elevation data covering all of Denmark
* @param {Object} options
* @param {Object} options.src - URL to download a pre-generated GeoTiff.
* @param {Object} options.auth - API autentication data. See ../config.js.example for reference.
* @param {Number} options.size - Width of the geoTiff image to return
* @returns GeoTIFF raster with elevation data
*/
function getDenmarkGeoTiff(options) {
if (options.src) {
return fetch(options.src)
.then(response => response.arrayBuffer())
.then(arrayBuffer => consumeGeoTIFF(arrayBuffer))
} else {
const auth = options.auth
const size = options.size ? options.size : 1000
return getTerrainGeoTIFF({bbox: [430000,6040000,900000,6413000]}, auth, null, size)
}
}

export {
get,
post,
getSTAC,
postSTAC,
getDHM,
getTerrainGeoTIFF
getTerrainGeoTIFF,
getDenmarkGeoTiff,
consumeGeoTIFF
}
4 changes: 2 additions & 2 deletions modules/saul-core.js
Original file line number Diff line number Diff line change
Expand Up @@ -202,8 +202,8 @@ async function compareZ(options) {
const delta = Math.abs(world_xyz[2] - options.z)

// Check delta vs. limit. If the difference is too big, try building coordinates with the new z
// Make sure the recursive loop ends at some point (iteration = 50) even if no precise z is found
if (options.iteration < 50 && delta > options.limit) {
// Make sure the recursive loop ends after 20 iterations even if no precise z is found
if (options.iteration < 20 && delta > options.limit) {
options.iteration++
options.z = world_xyz[2]
return compareZ(options)
Expand Down
120 changes: 87 additions & 33 deletions modules/saul-elevation.js
Original file line number Diff line number Diff line change
@@ -1,52 +1,106 @@
import { getTerrainGeoTIFF } from './api.js'
/**
* Maps a coordinate from the one coordinate system to another based on GeoTiff image dimensions.
* @param {Array} coord - EPSG:25832 X,Y coordinate pair
* @param {Array} bbox - Bounding box within input coordinate systtem. Consists of 2 EPSG:25832 coordinate pairs. Example: [x1,y1,x2,y2]
* @param {Array} imageDimensions - Width and height of a GeoTiff image. This will set the bounding box of the output coordinate system. Note: This coordinate system will have 0,0 set in the top right position of the image.
*/
function convertCoordinate(coord, bbox1, imageDimensions) {

// Extract the coordinates from the bounding boxes
const [x1_min, y1_min, x1_max, y1_max] = bbox1

// Extract the coordinate to convert
const [x1, y1] = coord

/** Constrain coordinates to be within a bounding box */
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]
// Calculate the width and height of both bounding boxes
const width1 = x1_max - x1_min
const height1 = y1_max - y1_min
const width2 = imageDimensions[0]
const height2 = imageDimensions[1]

// Calculate the relative position in the first coordinate system
const x_relative = (x1 - x1_min) / width1
const y_relative = (y1 - y1_min) / height1

// Calculate the new position in the second coordinate system
const x2 = Math.abs(x_relative * width2)
const y2 = height2 - Math.abs(y_relative * height2)

return [x2, y2]
}

/** Calculates nearest elevation for given coordinate using geoTIFF terrain data
* Heavily inspired by https://towardsdatascience.com/geotiff-coordinate-querying-with-javascript-5e6caaaf88cf
/**
* Calculates nearest elevation for given coordinate using geoTIFF terrain data
* @param {number} x - EPSG:25832 X coordinate
* @param {number} y - EPSG:25832 Y coordinate
* @param {object} geoTiff - GeoTiff data output from getTerrainGeoTIFF() method
* @param {object} geoTiff - GeoTiff data output from getTerrainGeoTIFF() or getDenmarkGeoTiff() method
* @returns {number} Elevation in meters
*/
async function getElevation(x, y, geoTiff) {

const bbox = geoTiff.getBoundingBox()
const height = geoTiff.getHeight()
const width = geoTiff.getWidth()
const bboxWidth = bbox[2] - bbox[0]
const bboxHeight = bbox[3] - bbox[1]

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 ) )
// Convert lat/lon to geotiff x/y
const [xPx, yPx] = convertCoordinate([x,y], bbox, [width, height])

// Fetch window of raster values around x/y
const window = [
Math.floor(xPx),
Math.floor(yPx),
Math.ceil(xPx),
Math.ceil(yPx)
]

// Fetch elevation value from geoTiff rasters
const elevation_data = await geoTiff.readRasters({
window: window,
fillValue: 0
})

return elevation_data[0][0]
}

const window = [ xPx, yPx, xPx + 1, yPx + 1 ]
const elevation_data = await geoTiff.readRasters( {window} )
/** Creates an ASCII visualization of a terrain GeoTiff. Useful for display in a terminal.
* NOTE: You may want to use a downsmapled GeoTiff to fit within your terminal window. A GeoTiff pixel is represented by 3 characters each.
* @param {object} gTiff - GeoTiff output from getTerrainGeoTIFF() or getDenmarkGeoTiff() method
*/
async function visualizeGeotiff(gTiff) {
console.log('--- GTIFF visualization ---')
const float32Arr = await gTiff.readRasters()
const tiffWidth = gTiff.getWidth()
const rasters = float32Arr[0]
let lines = []
let line = ''
for (let i = 0; i < rasters.length; i++) {
if (i % tiffWidth === 0) {
lines.push(line)
line = ''
}
line += zeroPadNumber(rasters[i])
}
for (const l of lines) {
console.log(l)
}
}

return Math.round( elevation_data[0][0] * 10) / 10
/** Converts numbers to strings in a 3-letter format. Ex. `002`, `0-3` (for negative numbers), or `123` */
function zeroPadNumber(input) {
if (input === 0) {
return "..."
} else if (input < 0) {
return `0${input.toFixed(0)}`
} else if (input < 10) {
return `00${Math.floor(input)}`
} else if (input < 100) {
return `0${Math.floor(input)}`
} else {
return `${Math.floor(input)}`
}
}

export {
getTerrainGeoTIFF,
getElevation
getElevation,
visualizeGeotiff
}
Loading

0 comments on commit 582b38b

Please sign in to comment.