kp/l2decodeV5
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Initial commit

Browse source
This commit is contained in:
kp 2025-02-19 14:28:12 -06:00
commit d052cb8152
8 changed files with 1907 additions and 0 deletions

6
.editorconfig Normal file
View file

@ -0,0 +1,6 @@
[*]
insert_final_newline = true
[*.{rs}]
indent_style = tab
indent_size = 4

1
.gitignore vendored Normal file
View file

@ -0,0 +1 @@
/target

1039
Cargo.lock generated Normal file
View file

File diff suppressed because it is too large Load diff

9
Cargo.toml Normal file
View file

@ -0,0 +1,9 @@
[package]
name = "l2decodeV5"
version = "0.1.0"
edition = "2021"
[dependencies]
byteorder = "1.4"
bzip2 = "0.5.1"
image = "0.25.5"

11
README.md Normal file
View file

@ -0,0 +1,11 @@
# l2decodeV5
This is a prototype WSR-88D Level-II (NEXRAD weather radar) data decoder and image generator written in Rust. I've never really used Rust before, so I decided to rewrite one of my old C routines to get a feel for how things work. This code probably isn't the best, but it's a start. I've come to like Rust a lot in the process of writing this over the past couple of days!
Notes:
- Only processes reflectivity data for the first elevation scan. The code is structured such that it's possible to get other products out of it though.
- Expected input is an Archive-II file assuming RDA build 19 and later. It probably wouldn't take much to get it to support data from sites running older builds of the RDA software.
## Output Example
![KRAX data c. 2020](https://i.imgur.com/2hn5zWR.png)

BIN
output.png Normal file
View file

Binary file not shown.
Side by side

After

(image error) Size: 910 KiB

841
src/main.rs Normal file
View file

@ -0,0 +1,841 @@
//
// NEXRAD Level-II data decoder and image generator with WSR-88D Build
// 10 Message 31 support
//
use std::io::{self, Read, Seek};
use std::fs::File;
use std::mem;
use std::str;
use std::cmp;
use byteorder::{BigEndian, ReadBytesExt};
use bzip2::read::{BzDecoder};
use image::{ImageBuffer, Rgba};
const MSG_LENGTH: usize = 2432;
const MOMENT_LENGTH: usize = 2400;
#[repr(u8)]
enum RadialStatus {
StartElevationScan,
InterRadial,
EndElevationScan,
StartVolumeScan,
EndVolumeScan
}
#[repr(u8)]
enum DataBlockType {
Volume,
Elevation,
Radial,
Reflectivity,
Velocity,
SpectrumWidth,
Zdr, // Differential reflectivity
Phi, // Differential phase
Rho, // Correlation coefficient
}
// --------
// ARCHIVE II STRUCTS
// These are implemented based on what's written in the ICD from NOAA:
// https://www.roc.noaa.gov/public-documents/icds/2620002Y.pdf
// --------
// At the start of every Archive II file is a 24-byte record
// describing certain attributes of the volume scan.
#[repr(C, packed)]
struct VolumeHeader {
ver: [u8; 2],
period: [u8; 1],
ext: [u8; 3],
modified_julian_date: i32,
modified_time_millis: i32,
icao_id: [u8; 4]
}
// Level II data packet -- used for all message types
#[repr(C, packed)]
struct Packet {
ctm: [i16; 6],
size: u16,
id: u8,
message_type: u8,
seq: u16,
gen_julian_date: u16,
gen_time_millis: u32,
num_message_segments: u16,
message_segment: u16,
data: [u8; 2404]
}
// Header data for message 31 (digital radar data)
#[repr(C, packed)]
#[derive(Copy, Clone)]
struct Message31Header {
site_id: [u8; 4],
collection_time: u32,
modified_julian_date: u16,
radial_idx: u16,
radial_azimuth: f32,
compression_indicator: u8,
spare: u8,
radial_length: u16,
azimuth_res_spacing_code: u8,
radial_status: u8,
elevation_number: u8,
cut_sector_number: u8,
elevation_angle: f32,
radial_spot_blanking_status: u8,
azimuth_indexing_mode: u8,
data_block_count: u16,
data_block_pointers: [u32; 9]
}
// Data block (descriptor of generic data moment type)
#[repr(C, packed)]
#[derive(Copy, Clone)]
struct DataBlock {
block_type: u8,
moment_name: [u8; 3]
}
impl Default for DataBlock {
fn default() -> Self {
Self {
block_type: 0,
moment_name: [0; 3]
}
}
}
#[repr(C, packed)]
#[derive(Default, Copy, Clone)]
struct VolumeData {
data_block: DataBlock,
size: u16,
ver_major: u8,
ver_minor: u8,
lat: f32,
lon: f32,
site_height: i16,
feedhorn_height: u16,
calibration_constant: f32,
hor_tx_power: f32,
ver_tx_power: f32,
system_diff_refl: f32,
initial_system_diff_refl: f32,
vcp_num: u16,
processing_status: u16
}
#[repr(C, packed)]
#[derive(Default, Copy, Clone)]
struct ElevationData {
data_block: DataBlock,
size: u16,
atmos: [u8; 2],
calibration_constant: f32
}
#[repr(C, packed)]
#[derive(Default, Copy, Clone)]
struct RadialData {
data_block: DataBlock,
size: u16,
unambiguous_range: u16,
noise_level_hor: f32,
noise_level_ver: f32,
nyquist_velocity: u16,
spare: [u8; 2],
calibration_constant_channel_hor: f32,
calibration_constant_channel_ver: f32
}
#[repr(C, packed)]
#[derive(Default, Copy, Clone)]
struct GenericDataMoment {
data_block: DataBlock,
reserved: u32,
num_gates: u16,
range: u16,
range_sample_interval: u16,
threshold: u16,
snr_threshold: u16,
control_flags: i8,
word_size: u8,
scale: f32,
offset: f32
}
#[repr(C, packed)]
#[derive(Copy, Clone)]
struct DataMoment {
generic_data_moment: GenericDataMoment,
data: [u16; MOMENT_LENGTH]
}
impl Default for DataMoment {
fn default() -> Self {
Self {
generic_data_moment: GenericDataMoment::default(),
data: [69; MOMENT_LENGTH]
}
}
}
// --------
// PROCESSED STRUCTS
// These aren't set to 1-byte alignment as these are *separate* from
// structs meant to be written to from raw data
// --------
#[derive(Copy, Clone)]
struct Message31 {
header: Message31Header,
volume_data: VolumeData,
elevation_data: ElevationData,
radial_data: RadialData,
reflectivity_data: DataMoment,
velocity_data: DataMoment,
spectrum_width_data: DataMoment,
zdr_data: DataMoment,
phi_data: DataMoment,
rho_data: DataMoment,
cfp_data: DataMoment
}
impl Message31 {
fn new(header: Message31Header) -> Self {
Self {
header,
volume_data: VolumeData::default(),
elevation_data: ElevationData::default(),
radial_data: RadialData::default(),
reflectivity_data: DataMoment::default(),
velocity_data: DataMoment::default(),
spectrum_width_data: DataMoment::default(),
zdr_data: DataMoment::default(),
phi_data: DataMoment::default(),
rho_data: DataMoment::default(),
cfp_data: DataMoment::default()
}
}
}
// Structure to hold radial data (radials), which contains momentary data (bins)
struct ElevationScan {
radials: Vec<Message31>,
azimuth: [f32; 800], // Expecting 720 radials, but
// overscan/underscan can occur...
first_radial_idx: u16,
first_radial_az: f32,
num_radials: i32,
num_gates: i32,
mean_elev_angle: f32,
accum_elev_angle: f32,
gate_size: f32,
first_gate: f32
}
impl ElevationScan {
fn new() -> Self {
Self {
radials: Vec::new(),
azimuth: [0.0; 800],
first_radial_idx: 0,
first_radial_az: 0.0,
num_radials: 0,
num_gates: 0,
mean_elev_angle: 0.0,
accum_elev_angle: 0.0,
gate_size: 0.0,
first_gate: 0.0
}
}
}
// Aggregate structure for any and all data coming from message packets
struct VolumeScan {
volume_header: VolumeHeader,
elevation_scans: Vec<ElevationScan>
}
impl VolumeScan {
fn new(volume_header: VolumeHeader) -> Self {
Self { volume_header, elevation_scans: Vec::new() }
}
}
fn check_for_volume_scan<R: Read + Seek>(mut stream: R) -> Option<VolumeScan> {
let magic_str = b"AR2V00";
let mut magic_pos = 0;
let mut buf = [0u8; 1];
while magic_pos < magic_str.len() {
if stream.read_exact(&mut buf).is_err() {
println!("Stream error; EOF?");
return None;
}
if buf[0] == magic_str[magic_pos] {
magic_pos += 1;
} else {
magic_pos = 0;
}
}
println!("Found AR2V00xx identifier");
Some(process_volume_scan(stream))
}
fn process_volume_scan<R: Read>(mut stream: R) -> VolumeScan {
// Read volume header
// ------------------
let mut header_buffer = [0u8; mem::size_of::<VolumeHeader>()];
stream.read_exact(&mut header_buffer);
let volume_header: VolumeHeader = unsafe {
std::ptr::read(header_buffer.as_ptr() as *const _)
};
// Test the ICAO ID from the volume header real quick...
// TODO: We should include some kind of check against an expected
// ICAO ID via command line arguments
let icao = str::from_utf8(&volume_header.icao_id).unwrap();
println!("Detected ICAO ID: {}", icao);
// Read rest of the data
// ---------------------
let mut volume_scan = VolumeScan::new(volume_header);
let mut control_word: i32 = 0;
let mut verify_control_word: i32 = 0;
let mut bail = false;
while !bail {
let mut raw_buffer = Vec::new();
let mut raw_buffer_size = 0;
let mut decompressed_data = Vec::new();
// Read in the control word (size) of the LDM record. There's
// also a chance this is a new volume scan, so we need to
// check for the "magic string" we encountered earlier.
let mut control_word: i32 = 0;
match stream.read_i32::<BigEndian>() {
Ok(res) => {
control_word = res;
},
Err(e) => {
bail = true;
}
}
if control_word == 0x56325241 /* "AR2V" but big endian */ {
// We've got another volume scan. We're bailing here, but
// really we ought to just start the decoding process over
// again from this byte.
println!("Encountered magic string twice!");
bail = true;
continue;
}
// Decompress the data
// -------------------
// TODO: check if it's actually bzip2-compressed first...?
if control_word >= raw_buffer_size {
raw_buffer.resize(control_word as usize, 0);
raw_buffer_size = control_word;
}
stream.read_exact(&mut raw_buffer);
let mut decoder = BzDecoder::new(&raw_buffer[..]);
let mut decoding = true;
let mut buffer = [0u8; 8192];
while decoding {
match decoder.read(&mut buffer) {
Ok(bytes_read) => {
if bytes_read == 0 {
decoding = false;
break;
}
decompressed_data.extend_from_slice(&buffer[..bytes_read]);
},
Err(error) => match error.kind() {
DataMagic => {
// If we've made it here then the bzip2 magic
// string is missing. In this case, we'll
// assume it's uncompressed data, although I
// feel like it'd probably be better to check
// for the string ahead of time prior to
// decompressing it. Doing it this way allows
// us to deal with Archive II files with data
// blocks that may or may not be
// bzip2-compressed though...
decompressed_data.extend_from_slice(&raw_buffer[..raw_buffer_size as usize]);
},
_ => panic!("Problem decompressing bzip2 data: {error:?}")
}
}
// If we have no data left, we can assume we're done.
if decompressed_data.len() == 0 {
decoding = false;
bail = true;
}
}
// Process decompressed data
// -------------------------
let mut msg_length: usize = MSG_LENGTH;
let mut i: usize = 0;
while i < decompressed_data.len() {
let slice = &decompressed_data[i..i+mem::size_of::<Packet>()];
let packet: Packet = unsafe {
std::ptr::read(slice.as_ptr() as *const _)
};
// Message 31 packets are variable length; we must
// readjust. If this isn't a message 31 packet though,
// then return to standard length.
if packet.message_type == 31 {
msg_length = u16::from_be(packet.size) as usize * 2 + 12;
} else {
msg_length = MSG_LENGTH;
}
if packet.message_type != 0 {
// Handle messages accordingly. We'll break out
// message processing into separate routines since
// some of them can be pretty involved (e.g. 1, 31)
match packet.message_type {
31 => {
// println!("Processing message 31... length: {}", msg_length);
let raw_data = &decompressed_data[i + mem::size_of::<Packet>() + mem::size_of::<Message31Header>() - packet.data.len()..];
bail |= process_message_31(&mut volume_scan, packet, raw_data);
if bail {
// Bail as soon as the elevation (not
// volume!) scan is done.
// This works for my purposes, but isn't
// isn't ideal for loading a whole archive
// at once.
render_image(&mut volume_scan);
}
},
_ => {
println!("Message {} functionality unimplemented...", packet.message_type);
}
}
}
i += msg_length;
}
}
return volume_scan;
}
fn convert_float(float: f32) -> f32 {
return f32::from_bits(float.to_bits().swap_bytes());
}
fn process_message_31(volume_scan: &mut VolumeScan, packet: Packet, raw_data: &[u8]) -> bool {
let slice = &packet.data[..mem::size_of::<Message31Header>()];
let header: Message31Header = unsafe {
std::ptr::read(slice.as_ptr() as *const _)
};
// Set up our new elevation scan (really not doing this the best
// way...)
let elevation_number: usize = header.elevation_number as usize;
if volume_scan.elevation_scans.len() < elevation_number {
volume_scan.elevation_scans.push(ElevationScan::new());
}
let elevation_scan: &mut ElevationScan = &mut volume_scan.elevation_scans[elevation_number - 1];
// Bump up the radial counter. This should contain the # of
// radials when we're done.
elevation_scan.num_radials += 1;
// Add the elevation angle of this radial (it varies) to the
// elevation scan structure
elevation_scan.accum_elev_angle += convert_float(header.elevation_angle);
elevation_scan.mean_elev_angle = elevation_scan.accum_elev_angle / elevation_scan.num_radials as f32;
// Note that NEXRAD radials are numbered starting at 1, but for
// our purposes we will number them starting at zero.
let radial: usize = u16::from_be(header.radial_idx) as usize - 1;
// Record the azimuth of the start of this radial. 'angle' is
// actually the center of the radial so to get the start we will
// subtract off 1/2 of the radial spacing (which will currently
// either be 1 or 0.5 degrees).
elevation_scan.azimuth[radial] = match header.azimuth_res_spacing_code {
1 => convert_float(header.radial_azimuth) - 0.25,
2 => convert_float(header.radial_azimuth) - 0.5,
_ => convert_float(header.radial_azimuth)
};
// If we go negative from the above operation, then correct
if elevation_scan.azimuth[radial] < 0.0 {
elevation_scan.azimuth[radial] = 360.0 + elevation_scan.azimuth[radial];
}
// Set our azimuth and index of the northernmost radial here.
// We'll also check if the azimuth of this radial is further
// clockwise than the ones we've already recorded. It's nice to
// know where the first radial is.
let wrap_around = elevation_scan.azimuth[radial] < elevation_scan.first_radial_az;
if radial == 0 || wrap_around {
elevation_scan.first_radial_az = elevation_scan.azimuth[radial];
elevation_scan.first_radial_idx = radial as u16;
}
// If we've wrapped around (i.e. crossed 0deg north) and the
// azimuth of this radial is further clockwise than the first one
// we've processed, then this indicates an overscan condition. In
// other words, we have scanned past 360 degrees. In this case,
// we'll effectively nullify this and all subsequent radials by
// decrementing the number of radials, making these not count.
if wrap_around && elevation_scan.azimuth[radial] > elevation_scan.azimuth[0] {
elevation_scan.num_radials -= 1;
}
// Build 19 added CFP (clutter filter power removed) data, an
// ABI-breaking change. Since our Message31Header structure is
// still set up to have only 9 data blocks, we'll need to check
// the RDA software version and bump up our byte index should this
// data be present. For now though, we'll hardcode it (cool)
let cfp_present_hack: usize = 4;
// Build 20 introduced another ABI-breaking change with the "RPG
// weighted mean ZDR bias estimate" field. It's 2 bytes w/ 6
// leftover as spare memory. We'll address it here too.
let rpg_weighted_mean_zdr_bias_present_hack = 0;
// Copy over the data to our block structure
let mut byte_idx = cfp_present_hack;
for i in 0..u16::from_be(header.data_block_count) {
let slice = &raw_data[byte_idx..byte_idx + mem::size_of::<DataBlock>()];
let data_block: DataBlock = unsafe { std::ptr::read(slice.as_ptr() as *const _) };
let name = str::from_utf8(&data_block.moment_name).unwrap();
if radial >= elevation_scan.radials.len() {
// This will be where the read data goes. We'll insert this into
// the elevation scan's radial data vector once we're done.
elevation_scan.radials.push(Message31::new(header));
}
match name {
"VOL" => {
byte_idx += mem::size_of::<VolumeData>();
},
"ELV" => {
byte_idx += mem::size_of::<ElevationData>();
},
"RAD" => {
byte_idx += mem::size_of::<RadialData>();
},
"REF"|"VEL"|"SW "|"ZDR"|"PHI"|"RHO"|"CFP" => {
let slice = &raw_data[byte_idx..byte_idx + mem::size_of::<GenericDataMoment>()];
let generic_moment: GenericDataMoment = unsafe {
std::ptr::read(slice.as_ptr() as *const _)
};
byte_idx += mem::size_of::<GenericDataMoment>();
// Per the ICD:
//
// LDM is the amount of space in bytes required for a
// data moment array and equals ((NG * DWS) / 8) where
// NG is the number of gates at the gate spacing
// resolution specified and DWS is the number of bits
// stored for each gate (DWS is always a multiple of
// 8).
//
// Means "Length of Data Moment"; no relation to
// UCAR's related LDM software :P
let num_gates = u16::from_be(generic_moment.num_gates);
let ldm = {num_gates * generic_moment.word_size as u16 / 8} as usize;
// Record adjusted gate information
elevation_scan.num_gates = cmp::max(elevation_scan.num_gates, num_gates as i32);
elevation_scan.gate_size = u16::from_be(generic_moment.range_sample_interval) as f32;
elevation_scan.first_gate = u16::from_be(generic_moment.range) as f32;
// Record radial data (gross nested match here but whatever)
let mut radial_block: &mut Message31 = &mut elevation_scan.radials[radial];
let mut moment: DataMoment = DataMoment::default();
let radar = &raw_data[byte_idx..];
for (i, chunk) in radar.chunks_exact(1).enumerate().take(MOMENT_LENGTH) {
moment.data[i] = chunk[0] as u16;
// For 16-bit data types (e.g. ZDR) we'd probably want this instead:
// moment.data[i] = u16::from_le_bytes([chunk[0], chunk[1]]);
}
match name {
"REF" => { radial_block.reflectivity_data = moment; },
"VEL" => { radial_block.velocity_data = moment; },
"SW " => { radial_block.spectrum_width_data = moment; },
"ZDR" => { radial_block.zdr_data = moment; },
"PHI" => { radial_block.phi_data = moment; },
"RHO" => { radial_block.rho_data = moment; },
"CFP" => { radial_block.cfp_data = moment; }
_ => {}
}
byte_idx += ldm;
}
_ => {
// panic!("Unexpected block type!");
}
}
// Debug:
// println!("Block name: {}", name);
}
// For our purposes we only want to process one elevation scan for
// now... if we're at the end of the elevation (not volume) scan,
// we're done.
if header.radial_status == RadialStatus::EndElevationScan as u8 {
return true;
}
return false;
}
fn get_dbz(r: f32, theta: f32, elevation_scan: &ElevationScan, smooth: bool) -> f32 {
let num_radials: i32 = elevation_scan.num_radials as i32;
let last_radial: i32 = num_radials - 1;
let native_scale_r = elevation_scan.num_gates as f32 / 458.0; // Gates per km
let native_scale_theta = elevation_scan.num_radials as f32 / 360.0; // Radials per deg
// Determine the range to the first gate that contains data (this
// is in the metadata from the radar). There is no data for the
// first km or two close to the radar site. It's the so-called
// "cone of silence". If we're asked for data in this region we'll
// return 0 dBZ.
let range_to_first_gate_km = elevation_scan.first_gate / 1000.0 - elevation_scan.gate_size / 1000.0 / 2.0;
if r <= range_to_first_gate_km {
return 0.0;
}
// Translate r (which is in km) to a gate number by subtracting
// off the range to the first gate data (which is also in km),
// then multiplying by the native_scale_r value (which is in
// gates/km). We need two gate numbers, one 'below' the actual r
// (i.e. closer to the radar site) and one 'above' the actual r
// (i.e. further from the radar site). These top and bottom gates
// will be used for interpolation.
let bottom = ((r - range_to_first_gate_km) * native_scale_r) as usize;
let top = (bottom + 1) as usize;
// Not only do we need the gate numbers for the top and bottom,
// but we also need the relative distance of r from the top and
// bottom gates. This distance is used as an interpolation weight,
// or to determine the nearest neighbor. For now we're only doing
// nearest neighbor interpolation.
let dist_bottom: f32 = (r - range_to_first_gate_km) * native_scale_r - bottom as f32;
let dist_top: f32 = 1.0 - dist_bottom;
// Identify the radial immediately preceeding the northernmost radial. It crosses the 0 degree threshold
let mut trans_radial = elevation_scan.first_radial_idx as i32;
if trans_radial < 0 {
trans_radial = last_radial;
}
// Here's a cool trick. Estimate the proper radial num (0 -
// lastradial) by subtracting the azimuth of the northernmost
// radial then taking the integer portion of the remainder. This
// will provide an index, relative to the northernmost radial, of
// the radial that we need. Once we have an index, we have to
// resolve it to an actual radial number by adding it to the
// radial number of the northernmost radial. If we exceed
// lastradial, then we have simply wrapped around.
let mut radial_idx: i32 = ((theta - elevation_scan.first_radial_az) * native_scale_theta + elevation_scan.first_radial_idx as f32) as i32;
if radial_idx > last_radial {
radial_idx -= num_radials;
}
// Check to see if our estimated radial index was correct. Usually
// it will be, however in some cases it might not. We check for
// overestimates and underestimates and starting at the estimate
// we work back and forth trying to find the true value. This was
// done to improve interpolation accuracy for some sites that have
// very uneven spacing between radial azimuths. A safety check is
// provided to prevent hangs. Also, the transition radial (the one
// right before 0 degrees) is exempted from this check for obvious
// reasons. This could be made better, but meh.
let mut safety = 10;
while theta < elevation_scan.azimuth[radial_idx as usize] && radial_idx != trans_radial && safety != 0 {
safety -= 1;
if radial_idx > 0 {
radial_idx -= 1;
}
if radial_idx < 0 {
radial_idx = last_radial;
}
}
safety = 10;
while theta > elevation_scan.azimuth[radial_idx as usize + 1] && radial_idx != trans_radial && safety != 0 {
safety -= 1;
radial_idx += 1;
if radial_idx > last_radial {
radial_idx -= num_radials;
}
}
if safety == 0 {
return 0.0;
}
// Determine the left and right radial numbers to use for interpolation.
let mut left: usize = 0;
let mut right: usize = 0;
let mut dist_left: f32 = 0.0;
let mut dist_right: f32 = 0.0;
if radial_idx == trans_radial {
// We're processing the radial immediately before the northernmost one
left = trans_radial as usize;
right = elevation_scan.first_radial_idx as usize;
if theta < 180.0 {
dist_left = theta + 360.0 - elevation_scan.azimuth[left];
dist_right = elevation_scan.azimuth[right] - theta;
} else {
dist_left = theta - elevation_scan.azimuth[left];
dist_right = elevation_scan.azimuth[right] + 360.0 - theta;
}
} else {
left = radial_idx as usize;
right = left + 1;
if right > last_radial as usize {
right = 0;
}
dist_left = theta - elevation_scan.azimuth[left];
dist_right = elevation_scan.azimuth[right] - theta;
}
// Grab the dBZ values for all 4 neighboring points
let mut dbz_ll = elevation_scan.radials[left].reflectivity_data.data[bottom] as f32;
let mut dbz_ul = elevation_scan.radials[left].reflectivity_data.data[top] as f32;
let mut dbz_lr = elevation_scan.radials[right].reflectivity_data.data[bottom] as f32;
let mut dbz_ur = elevation_scan.radials[right].reflectivity_data.data[top] as f32;
// Apply a dBZ threshold (don't want anything below 0 for reasons)
dbz_ll = f32::max(dbz_ll, 66.0);
dbz_ul = f32::max(dbz_ul, 66.0);
dbz_lr = f32::max(dbz_lr, 66.0);
dbz_ur = f32::max(dbz_ur, 66.0);
// Perform interpolation
let mut dbz_top: f32 = 0.0;
let mut dbz_bottom: f32 = 0.0;
let mut dbz: f32 = 0.0;
if smooth {
dbz_top = ((dbz_ul * dist_right) + (dbz_ur * dist_left)) / (dist_left + dist_right);
dbz_bottom = ((dbz_ll * dist_right) + (dbz_lr * dist_left)) / (dist_left + dist_right);
dbz = ((dbz_top * dist_bottom) + (dbz_bottom * dist_top)) / (dist_top + dist_bottom);
} else {
if dist_left < dist_right {
dbz_top = dbz_ul;
dbz_bottom = dbz_ll;
} else {
dbz_top = dbz_ur;
dbz_bottom = dbz_lr;
}
if dist_top < dist_bottom {
dbz = dbz_top;
} else {
dbz = dbz_bottom;
}
}
return (dbz - 2.0) / 2.0 - 32.0;
}
const NUMINTERVALS: usize = 8;
fn render_image(volume_scan: &mut VolumeScan) {
let elevation_scan: &ElevationScan = &volume_scan.elevation_scans[0];
// We *could* use a 16-color scale, but this 128-color interpolated one is cooler
let mut color_scale: [[u8; 3]; 16 * NUMINTERVALS] = [[0; 3]; 16 * NUMINTERVALS];
build_smooth_color_scale(&mut color_scale);
// Render :D
let width = 1920;
let height = 1920;
let mut img: ImageBuffer<Rgba<u8>, Vec<u8>> = ImageBuffer::new(width, height);
let mut x1: f32 = 0.0;
let mut y1: f32 = 0.0;
for (x, y, pixel) in img.enumerate_pixels_mut() {
x1 = x as f32 - width as f32 / 2.0;
y1 = height as f32 / 2.0 - y as f32;
let r = f32::sqrt(x1 * x1 + y1 * y1) * 0.25;
let mut theta = (1.570795 - f32::atan2(y1, x1)) * 57.29583; // (90deg-atan2(y1,x1))*180/pi)
if theta < 0.0 {
theta = 360.0 + theta;
}
let mut dbz = get_dbz(r, theta, elevation_scan, false);
let color_idx: usize = (dbz * NUMINTERVALS as f32 / 5.0 + 0.5) as usize;
if color_idx >= 16 * NUMINTERVALS || dbz <= 0.0 {
continue;
}
let out_r = color_scale[color_idx][0];
let out_g = color_scale[color_idx][1];
let out_b = color_scale[color_idx][2];
*pixel = Rgba([out_r, out_g, out_b, 0xFF]);
}
img.save("output.png").expect("Failed to save image");
}
fn build_smooth_color_scale(scolorscale: &mut [[u8; 3]; 16 * NUMINTERVALS]) {
let colorscale: [[u8; 3]; 16] = [
[0, 0, 0],
[0, 72, 144],
[0, 96, 240],
[128, 224, 80],
[100, 184, 64],
[72, 144, 48],
[44, 104, 32],
[16, 64, 16],
[240, 160, 16],
[240, 118, 32],
[240, 16, 32],
[144, 0, 0],
[176, 32, 128],
[202, 64, 160],
[229, 96, 192],
[225, 128, 224],
];
let numintervals = NUMINTERVALS as f32;
let mut step = 0;
while step < 15 {
let (curr, curg, curb) = (colorscale[step][0], colorscale[step][1], colorscale[step][2]);
let (nextr, nextg, nextb) = (colorscale[step + 1][0], colorscale[step + 1][1], colorscale[step + 1][2]);
// Compute the deltas between the colors
let deltar = (nextr as f32 - curr as f32) / numintervals;
let deltag = (nextg as f32 - curg as f32) / numintervals;
let deltab = (nextb as f32 - curb as f32) / numintervals;
for t in 0..numintervals as usize {
scolorscale[step * numintervals as usize + t][0] = (curr as f32 + deltar * t as f32).round() as u8;
scolorscale[step * numintervals as usize + t][1] = (curg as f32 + deltag * t as f32).round() as u8;
scolorscale[step * numintervals as usize + t][2] = (curb as f32 + deltab * t as f32).round() as u8;
}
step += 1;
}
scolorscale[step * numintervals as usize][0] = colorscale[step][0];
scolorscale[step * numintervals as usize][1] = colorscale[step][1];
scolorscale[step * numintervals as usize][2] = colorscale[step][2];
// Last bit of color scale is all white
for t in 1..numintervals as usize {
scolorscale[step * numintervals as usize + t][0] = (248 + t) as u8;
scolorscale[step * numintervals as usize + t][1] = (248 + t) as u8;
scolorscale[step * numintervals as usize + t][2] = (248 + t) as u8;
}
}
fn main() -> io::Result<()> {
println!("Processing file...");
let file = File::open("test/KRAX.dat")?;
if let Some(scan) = check_for_volume_scan(file) {
println!("Successfully processed volume scan");
} else {
println!("Failed to find volume scan data");
}
Ok(())
}

BIN
test/KRAX.dat Normal file
View file

Binary file not shown.