Terrain rgb

[kylebarron]

https://github.com/tmnnrs/terrain-rgb

gdalbuildvrt -vrtnodata -9999 -input_file_list terrain50_filelist.txt terrain50.vrt
gdalwarp -r cubicspline -s_srs EPSG:4269 -t_srs EPSG:3857 -dstnodata 0 -co COMPRESS=DEFLATE USGS_NED_13_n34w117_IMG.img terrain_cubicspline.tif
rio rgbify -b -10000 -i 0.1 --max-z 9 --min-z 5 --format png terrain_cubicspline.tif terrain.mbtiles

[kylebarron]

By the way, it looks like it makes 512x512 tiles by default.

[kylebarron]

For manually getting the data out of the tiles, you can probably use https://github.com/mapbox/tilebelt to find the locations of each tile:

> lon = -116.6794204
-116.6794204
> lat = 33.8143198
33.8143198
> tilebelt.pointToTile(lon, lat, 12)
[ 720, 1638, 12 ]
> tilebelt.pointToTileFraction(lon, lat, 12)
[ 720.447483448889, 1638.7684395282195, 12 ]

I think the tile fraction should let you get the pixel within the file that you want:

> x = .447483448889 * 512
229.111525831168
> y = .7684395282195 * 512
393.441038448384
> var getPixels = require("get-pixels")
undefined

var path = '/Users/kyle/github/mapping/nst-guide/hillshade/data/tmp-rgbify/terrain_cubicspline/12/720/1638.png'

let data
getPixels(path, function(err, pixels) {
  if(err) {
    console.log("Bad image path")
    return
  }
  data = pixels
  console.log("got pixels", pixels.shape.slice())
})

x = 229.111525831168
y = 393.441038448384
x = 229
y = 393

R = data.get(x, y, 0)
G = data.get(x, y, 1)
B = data.get(x, y, 2)

height = -10000 + ((R * 256 * 256 + G * 256 + B) * 0.1)
// > 3277.9 meters -> 10754 feet
-116.6794204
33.8143198
// https://nationalmap.gov/epqs/pqs.php?x=-116.6794204&y=33.8143198&units=Meters&output=json
// {"USGS_Elevation_Point_Query_Service":{"Elevation_Query":{"x":-116.6794204,"y":33.8143198,"Data_Source":"3DEP 1\/3 arc-second","Elevation":3281.75,"Units":"Meters"}}}

[kylebarron]

And with bilinear interpolation:

var getPixels = require("get-pixels")
var path = '/Users/kyle/github/mapping/nst-guide/hillshade/data/tmp-rgbify/terrain_cubicspline/12/720/1638.png'

let data
getPixels(path, function(err, pixels) {
  if(err) {
    console.log("Bad image path")
    return
  }
  data = pixels
  console.log("got pixels", pixels.shape.slice())
})

x = 229.111525831168
y = 393.441038448384
x_arr = [Math.floor(x), Math.ceil(x)]
y_arr = [Math.floor(y), Math.ceil(y)]

// Arrays must be in order of:
// val00, val01, val10, val11
let r_arr = []
let g_arr = []
let b_arr = []
for (const i of x_arr) {
  for (const j of y_arr) {
    r_arr.push(data.get(i, j, 0))
    g_arr.push(data.get(i, j, 1))
    b_arr.push(data.get(i, j, 2))
  }
}

const allEqual = arr => arr.every( v => v === arr[0] )

const bi_interp = (data, x, y) => {
  if (allEqual(data)) {
    return data[0]
  }

  // https://en.wikipedia.org/wiki/Bilinear_interpolation#Unit_square
  var linear_x = data[0]*(1-x)*(1-y) + data[1]*(1-x)*y + data[2]*x*(1-y) + data[3]*x*y
  return linear_x;
}

R = bi_interp(r_arr, x % 1, y % 1)
G = bi_interp(g_arr, x % 1, y % 1)
B = bi_interp(b_arr, x % 1, y % 1)

height = -10000 + ((R * 256 * 256 + G * 256 + B) * 0.1)
// > 3274.7 meters -> 10743 feet
// Correct height:
// https://nationalmap.gov/epqs/pqs.php?x=-116.6794204&y=33.8143198&units=Meters&output=json
// {"USGS_Elevation_Point_Query_Service":{"Elevation_Query":{"x":-116.6794204,"y":33.8143198,"Data_Source":"3DEP 1\/3 arc-second","Elevation":3281.75,"Units":"Meters"}}}

[kylebarron]

code

[CrazyFluffyPony]

Seems like manually editing this code has fixed it for me (no guarantees whatsoever guys):

/usr/local/lib/python3.9/dist-packages/rio_rgbify/mbtiler.py LINE 195

    w, s, e, n = transform_bounds(*[src_crs, "epsg:4326"] + bbox, densify_pts=0)

change it to:

    w, s, e, n = transform_bounds(*[src_crs, "epsg:4326"] + bbox, densify_pts=2)

the file depends on your python version obviously