Open kylebarron opened 4 years ago
By the way, it looks like it makes 512x512 tiles by default.
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"}}}
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"}}}
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
https://github.com/tmnnrs/terrain-rgb