66 lines
2.1 KiB
Ruby
66 lines
2.1 KiB
Ruby
# frozen_string_literal: true
|
|
|
|
require "zlib"
|
|
|
|
module MoonModel
|
|
class DataSet
|
|
MAGIC = "MOONDATA1\n".b
|
|
DATUM_RADIUS_M = 1_737_400.0
|
|
|
|
attr_reader :width, :height, :metadata
|
|
|
|
def self.default_path
|
|
File.expand_path("../../data/moon_surface.bin.gz", __dir__)
|
|
end
|
|
|
|
def initialize(path = self.class.default_path)
|
|
raise ArgumentError, "lunar data asset not found: #{path}" unless File.file?(path)
|
|
|
|
raw = Zlib::GzipReader.open(path, &:read)
|
|
raise ArgumentError, "unrecognized lunar data asset" unless raw.start_with?(MAGIC)
|
|
|
|
offset = MAGIC.bytesize
|
|
header_length = raw.byteslice(offset, 4).unpack1("L<")
|
|
offset += 4
|
|
@metadata = JSON.parse(raw.byteslice(offset, header_length))
|
|
offset += header_length
|
|
@width = metadata.fetch("width")
|
|
@height = metadata.fetch("height")
|
|
elevation_bytes = width * height * 2
|
|
@elevation = raw.byteslice(offset, elevation_bytes).unpack("s<*")
|
|
@reflectance = raw.byteslice(offset + elevation_bytes, width * height).bytes
|
|
end
|
|
|
|
def elevation_m(latitude, longitude)
|
|
DATUM_RADIUS_M + sample(@elevation, latitude, longitude)
|
|
end
|
|
|
|
def reflectance(latitude, longitude)
|
|
value = sample(@reflectance, latitude, longitude)
|
|
low, high = metadata.fetch("reflectance_percentiles", [0.0, 255.0])
|
|
[[(value - low) / (high - low), 0.0].max, 1.0].min
|
|
end
|
|
|
|
def relief_range_m
|
|
minmax = metadata["elevation_offset_range_m"]
|
|
minmax ? minmax[1] - minmax[0] : @elevation.minmax.then { |a, b| b - a }
|
|
end
|
|
|
|
private
|
|
|
|
def sample(values, latitude, longitude)
|
|
x = (longitude % 360.0) / 360.0 * width
|
|
y = (90.0 - [[latitude, 90.0].min, -90.0].max) / 180.0 * (height - 1)
|
|
x0 = x.floor % width
|
|
x1 = (x0 + 1) % width
|
|
y0 = y.floor
|
|
y1 = [y0 + 1, height - 1].min
|
|
tx = x - x.floor
|
|
ty = y - y.floor
|
|
a = values[y0 * width + x0] * (1.0 - tx) + values[y0 * width + x1] * tx
|
|
b = values[y1 * width + x0] * (1.0 - tx) + values[y1 * width + x1] * tx
|
|
a * (1.0 - ty) + b * ty
|
|
end
|
|
end
|
|
end
|