Example of resolving GPS lat/long to elevation from a WGS84 TIFF

[Pomax]

The examples current show how to map tile corners pixel indices to GPS coordinates, but it does not show how to do the reverse, which is one of the most important things GDAL code should be able to do =)

Can the example be amended, or can a second example be written, that shows up how to start with a lat/long coordinate and getting the associated elevation data for that coordinate based on WGS84 (or tiff-stored) projection?

E.g. in the python version of GDAL, the following does the trick:

class ALOSTile():
    def __init__(self, tile_path):
        """
        Load a tile, and cache it. Individual tiles are relatively
        small, so for performance we preload the entire tile's
        elevation data into RAM.
        """
        self.tile_path = tile_path
        self.dataset = gdal.Open(tile_path, gdal.GA_ReadOnly)

        if self.dataset is None:
            raise Exception(f'Could not load GDAL file{tile_path}')

        src = osr.SpatialReference()
        src.SetWellKnownGeogCS("WGS84")
        dest = osr.SpatialReference(self.dataset.GetProjection())
        self.ct = osr.CoordinateTransformation(src, dest)
        self.grid = self.dataset.GetRasterBand(1).ReadAsArray()

    def lookup(self, lat, lon):
        """
        see https://gis.stackexchange.com/a/415337/219296
        """
        try:
            forward_transform = self.dataset.GetGeoTransform()
            reverse_transform = gdal.InvGeoTransform(forward_transform)
            x, y = [int(v) for v in gdal.ApplyGeoTransform(
                reverse_transform, lon, lat)]
            return int(self.grid[y][x])

        except:
            return None

But I don't see an InvGeoTransform in this codebase anywhere, so... how would one do this?

[jfoclpf]

why are you posting code in python if this lib is for Javascript?

[Pomax]

...did you read the text? Which explains exactly why I'm showing GDAL code that uses the official GDAL Python API?

It illustrates something that folks can do in "regular GDAL", and so is something people would expect to also be able to do in node-gdal but is seemingly unsupported. Which is odd given that node-gdal is supposedly Node bindings for GDAL itself.

[jfoclpf]

...did you read the text? Which explains exactly why I'm showing GDAL code that uses the official GDAL Python API?

It illustrates something that folks can do in "regular GDAL", and so is something people would expect to also be able to do in node-gdal but is seemingly unsupported. Which is odd given that node-gdal is supposedly Node bindings for GDAL itself.

Ok, got it, but this project seems abandoned

[Pomax]

It sure does.

I ended up just digging into the GeoTiIFF image format specification and wrote my own code for working with elevation data on top of the plain tiff package, although geotiff might be a better choice if you don't want to manually parse the geokeys byte dictionary.

I just submitted a PR to geotiff over on https://github.com/geotiffjs/geotiff.js/pull/353/files to show how to use it to resolve GPS coordinates to elevation values, but I essentially gave up trying to find a decent node implementation of the GDAL API.

[jfoclpf]

nice contribution @Pomax

exactly what I was looking for, thank you; if you may want to reply also here, you'll get the points and the prize :)

[jfoclpf]

@Pomax but I realized now you use another library geotiff, I will check it out

[Pomax]

yeah, I think the solution here is mostly "don't use node-gdal" although you might be able to construct the matrices using this project... see https://gis.stackexchange.com/questions/384221/calculating-inverse-polynomial-transforms-for-pixel-sampling-when-map-georeferen/452575#452575 for some more details on that.

[jfoclpf]

Btw @Pomax to install and use that geoTiff npm package in production do you need extra SW like imagemagick or it's ready to use after npm install geotiff? The instructions mentions some further SW installations and I got confused 😕

[Pomax]

That's if you want to build it from source with the full test suite it uses, which you don't need to. You just want to install it with npm.

[jfoclpf]

BTW @Pomax I found another solution using a fork of this repo, it works (I'm just stuck with coordinate transformation, but it works), here: https://github.com/mmomtchev/node-gdal-async/blob/master/examples/serve.js

[Pomax]

I used gdal-async for a little bit, but for what I need (elevation from WGS084 data) turns out that just parsing the tiff is actually enough, so that saves me having to use a larger library =)

Like...a lot larger. It's huge. This one's ~75MB, gdal-async is almost 200MB. The geotiff package is less than 2mb... and the simplest possible option, tiff, is 180kb! So... yeah I'm using tiff! =D

import {
  existsSync,
  readFileSync,
  readdirSync,
  writeFileSync,
  copyFileSync,
} from "fs";
import { mkdir } from "fs/promises";
import { basename, join } from "path";
import url from "url";
import tiff from "tiff";

const __dirname = url.fileURLToPath(new URL(".", import.meta.url));
const { floor, ceil, max } = Math;

const SEA_LEVEL = 0;
const ALOS_VOID_VALUE = -9999;
const CACHE_DIR = join(__dirname, `cache`);

await mkdir(CACHE_DIR, { recursive: true });

class ALOSTile {
  constructor(tilePath, coarseLevel = 10) {
    this.tilePath = tilePath;
    this.coarseLevel = coarseLevel;
    // copy to cache dir for faster loading
    const filename = join(`.`, CACHE_DIR, basename(tilePath));
    if (!existsSync(filename)) copyFileSync(tilePath, filename);
    this.init(filename);
  }

  init(filename) {
    const file = readFileSync(filename);
    const image = tiff.decode(file.buffer);
    const block = (this.block = image[0]);
    const fields = block.fields;
    let [sx, sy, sz] = fields.get(33550);
    let [px, py, k, gx, gy, gz] = fields.get(33922);
    sy = -sy;
    this.forward = [gx, sx, 0, gy, 0, sy];
    this.reverse = [-gx / sx, 1 / sx, 0, -gy / sy, 0, 1 / sy];
    this.pixels = block.data;
    const params = [sx, sy, gx, gy];
    this.formCoarseTile(block.width, block.height, params);
  }

  formCoarseTile(width, height, [sx, sy, gx, gy]) {
    // form a much smaller, coarse lookup map
    const { coarseLevel, pixels: p } = this;
    this.coarsePixels = [];
    for (let i = 0; i < p.length; i += coarseLevel) {
      this.coarsePixels[i / coarseLevel] = max(...p.slice(i, i + coarseLevel));
    }
    const [w, h] = [width / coarseLevel, height / coarseLevel];
    for (let i = 0; i < w; i += coarseLevel) {
      let list = [];
      for (let j = 0; j < coarseLevel; j++) list.push(p[i + j * w]);
      this.coarsePixels[i / coarseLevel] = max(...list);
    }
    this.coarsePixels = new Uint16Array(this.coarsePixels);
    const [sxC, syC] = [sx * coarseLevel, sy * coarseLevel];
    this.coarseForward = [gx, sxC, 0, gy, 0, syC];
    this.coarseReverse = [-gx / sxC, 1 / sxC, 0, -gy / syC, 0, 1 / syC];
  }

  pixelToGeo(x, y, coarse = false) {
    // returns [lat, long] (it does NOT return [long, lat]!)
    const F = coarse ? this.coarseForward : this.forward;
    return [F[3] + F[4] * x + F[5] * y, F[0] + F[1] * x + F[2] * y];
  }

  geoToPixel(lat, long, coarse = false) {
    // returns [x, y]
    const R = coarse ? this.coarseReverse : this.reverse;
    return [R[0] + R[1] * long + R[2] * lat, R[3] + R[4] * long + R[5] * lat];
  }

  lookup(lat, long, coarse = false) {
    const [x, y] = this.geoToPixel(lat, long, coarse);
    const pos = (x | 0) + (y | 0) * this.block.width;
    let value = (coarse ? this.coarseForward : this.pixels)[pos];
    // the highest point on earth is 8848m
    if (value === undefined || value > 10000) value = ALOS_VOID_VALUE;
    return value;
  }
}

[jfoclpf]

my GeoTiffs have 2Gb, thus 200Mb makes no great difference ;)

[Pomax]

Of course it does. Why would you host your tiff files with the code that accesses them? This is what large file network locations were invented for =)

(Also holy crap why are your tiffs 2GB? That's insane, start using tiles =)

[jfoclpf]

I am using tiles, what I meant is that the total amount of tiffs makes 2Gb

Thanks for the tip regarding storage of large files

[Pomax]

Oh, yeah then that's nothing. I'm using the ALOS World 3D (30m) dataset, it's 450GB, but each tile's only 25MB. So the data lives on the network, and the code runs on my computer, only locally caching the tiles it actually needs. The problem isn't just "disk space", it's also the fact that you never just run npm install once, you always run it more than you want to, and having to pull down 200MB each time it needs to properly rerun is pretty crazy. I'll take the 100x smaller library every time =)