radar-g/loaders/etws_loader/src/raw.rs
2024-08-28 19:25:11 +08:00

754 lines
22 KiB
Rust

use abi_stable::std_types::vec;
use bytemuck::cast_slice;
use nom::{
bytes::complete::take,
combinator::map,
number::complete::{le_f32, le_i16, le_i64, le_u16, le_u32},
sequence::tuple,
IResult,
};
use zip::read::ZipArchive;
use std::{
collections::HashMap,
io::{self, Cursor, Read, Seek, SeekFrom},
process::id,
};
use std::{fs::File, path::Path};
const SUPPORTED_TYPE: [u32; 7] = [2, 3, 4, 7, 9, 10, 11];
#[derive(Debug)]
pub struct RadarData {
pub datetime: i64,
pub els: Vec<f64>,
pub azs: Vec<f64>,
pub ranges: Vec<f64>,
pub datas: HashMap<&'static str, Vec<f64>>,
}
#[derive(Debug)]
struct SiteConfig {
site_code: String,
site_name: String,
latitude: f32,
longitude: f32,
antenna_height: u32,
}
#[derive(Debug)]
struct TaskConfig {
task_name: String,
task_description: String,
polarization_type: u32,
scan_type: u32,
scan_beam_number: u32,
cut_number: u32,
ray_order: u32,
scan_start_time: i64,
}
#[derive(Debug)]
struct BeamConfig {
cut_index: i16,
tx_beam_index: i16,
rx_beam_elevation: f32,
tx_beam_gain: f32,
rx_beam_width_h: f32,
rx_beam_width_v: f32,
rx_beam_gain: f32,
process_mode: String,
wave_form: String,
n1_prf1: f32,
n1_prf2: f32,
n2_prf1: f32,
n2_prf2: f32,
dealiasing_mode: u32,
azimuth: f32,
start_angle: f32,
end_angle: f32,
angular_resolution: f32,
scan_speed: f32,
log_resolution: f32,
doppler_resolution: f32,
maximum_range1: u32,
maximum_range2: u32,
start_range: u32,
sample1: u32,
sample2: u32,
phase_mode: u32,
atmospheric_loss: f32,
nyquist_speed: f32,
moments_mask: i64,
moments_size_mask: i64,
misc_filter_mask: u32,
sqi_threshold: f32,
sig_threshold: f32,
csr_threshold: f32,
log_threshold: f32,
cpa_threshold: f32,
pmi_threshold: f32,
dplog_threshold: f32,
thresholds_r: [u8; 4],
dbt_mask: u32,
dbz_mask: u32,
velocity_mask: u32,
spectrum_width_mask: u32,
dp_mask: u32,
direction: u32,
ground_clutter_classifier_type: i16,
ground_clutter_filter_type: i16,
ground_clutter_filter_notch_width: i16,
ground_clutter_filter_window: i16,
}
#[derive(Debug)]
struct TypeHeader {
data_type: u32,
scale: u32,
offset: u32,
bin_length: u16,
flags: u16,
length: u32,
}
#[derive(Debug)]
struct RadialData {
radial_state: u32,
spot_blank: u32,
sequence_number: u32,
radial_number: u32,
elevation_number: u32,
azimuth: f32,
elevation: f32,
seconds: i64,
microseconds: u32,
length_of_data: u32,
moment_number: u32,
scan_beam_index: i16,
horizontal_estimated_noise: i16,
vertical_estimated_noise: i16,
prf_flag: u32,
}
#[derive(Debug)]
struct GenericHeader {
magic_number: u32,
major_version: u16,
minor_version: u16,
generic_type: u32,
product_type: u32,
reserved: [u8; 16],
}
#[derive(Debug)]
struct SubPulseConfig {
strategy: u32,
modulation: u32,
frequency: f32,
bandwidth: f32,
width: u32,
horizontal_noise: f32,
vertical_noise: f32,
horizontal_calibration: f32,
vertical_calibration: f32,
horizontal_noise_temperature: f32,
vertical_noise_temperature: f32,
zdr_calibration: f32,
phidp_calibration: f32,
ldr_calibration: f32,
pulse_points: i16,
}
#[derive(Debug)]
struct ScanConfig {
beam_index: u32,
beam_type: u32,
sub_pulse_number: u32,
tx_beam_direction: f32,
tx_beam_width_h: f32,
tx_beam_width_v: f32,
tx_beam_gain: f32,
sub_pulse_configs: Vec<SubPulseConfig>,
}
fn parse_scan_config(input: &[u8]) -> IResult<&[u8], ScanConfig> {
let (input, beam_index) = le_u32(input)?;
let (input, beam_type) = le_u32(input)?;
let (input, sub_pulse_number) = le_u32(input)?;
let (input, tx_beam_direction) = le_f32(input)?;
let (input, tx_beam_width_h) = le_f32(input)?;
let (input, tx_beam_width_v) = le_f32(input)?;
let (input, tx_beam_gain) = le_f32(input)?;
let (input, _) = take(100usize)(input)?; // Skip reserved bytes
// Parse the SubPulseConfigs
let mut sub_pulse_configs = Vec::new();
let mut input = input;
for _ in 0..4 {
let (new_input, sub_pulse_config) = parse_sub_pulse_config(input)?;
input = new_input;
sub_pulse_configs.push(sub_pulse_config);
}
let scan_config = ScanConfig {
beam_index,
beam_type,
sub_pulse_number,
tx_beam_direction,
tx_beam_width_h,
tx_beam_width_v,
tx_beam_gain,
sub_pulse_configs,
};
Ok((input, scan_config))
}
fn parse_sub_pulse_config(input: &[u8]) -> IResult<&[u8], SubPulseConfig> {
let (input, strategy) = le_u32(input)?;
let (input, modulation) = le_u32(input)?;
let (input, frequency) = le_f32(input)?;
let (input, bandwidth) = le_f32(input)?;
let (input, width) = le_u32(input)?;
let (input, horizontal_noise) = le_f32(input)?;
let (input, vertical_noise) = le_f32(input)?;
let (input, horizontal_calibration) = le_f32(input)?;
let (input, vertical_calibration) = le_f32(input)?;
let (input, horizontal_noise_temperature) = le_f32(input)?;
let (input, vertical_noise_temperature) = le_f32(input)?;
let (input, zdr_calibration) = le_f32(input)?;
let (input, phidp_calibration) = le_f32(input)?;
let (input, ldr_calibration) = le_f32(input)?;
let (input, pulse_points) = nom::number::complete::le_i16(input)?;
let (input, _) = take(70usize)(input)?; // Skip reserved bytes
let sub_pulse_config = SubPulseConfig {
strategy,
modulation,
frequency,
bandwidth,
width,
horizontal_noise,
vertical_noise,
horizontal_calibration,
vertical_calibration,
horizontal_noise_temperature,
vertical_noise_temperature,
zdr_calibration,
phidp_calibration,
ldr_calibration,
pulse_points,
};
Ok((input, sub_pulse_config))
}
pub fn parse_raw_data<P: AsRef<Path>>(path: P) -> io::Result<RadarData> {
let path = path.as_ref();
if path.extension().unwrap() == "zip" {
let file = File::open(path).unwrap();
let mut archive = ZipArchive::new(file).unwrap();
let mut file = archive.by_index(0).unwrap();
let mut buf = Vec::new();
file.read_to_end(&mut buf).unwrap();
if &buf[0..4] == b"SSTM" || &buf[16..20] == b"STDM" || &buf[0..4] == b"RSTM" {
return get_radar_data(buf);
} else {
panic!("Invalid file format");
}
} else {
let mut file = File::open(path).unwrap();
let mut buf = Vec::new();
file.read_to_end(&mut buf).unwrap();
if &buf[0..4] == b"SSTM" || &buf[16..20] == b"STDM" || &buf[0..4] == b"RSTM" {
return get_radar_data(buf);
} else {
panic!("Invalid file format");
}
}
}
fn get_radar_data(data: Vec<u8>) -> io::Result<RadarData> {
let total_length = data.len() as u64;
// 使用 Cursor 来处理 Vec<u8>
let mut cursor = Cursor::new(data);
// 读取 GENERIC_HEADER
let mut buffer = vec![0u8; 32];
cursor.read_exact(&mut buffer)?;
let (_, generic_header) = parse_generic_header(&buffer).unwrap();
// 读取 SITE_CONFIG
let mut buffer = vec![0u8; 128];
cursor.read_exact(&mut buffer)?;
let (_, site_config) = parse_site_config(&buffer).unwrap();
// 读取 TASK_CONFIG
let mut buffer = vec![0u8; 256];
cursor.read_exact(&mut buffer)?;
let (_, task_config) = parse_task_config(&buffer).unwrap();
// 读取 SCAN_CONFIG
let mut scan_configs = Vec::new();
for _ in 0..task_config.scan_beam_number {
let mut buffer = vec![0u8; 640];
cursor.read_exact(&mut buffer)?;
let (_, scan_config) = parse_scan_config(&buffer).unwrap();
scan_configs.push(scan_config);
}
// 读取 BEAM_CONFIG
let mut beam_configs = Vec::new();
for _ in 0..task_config.cut_number {
let mut buffer = vec![0u8; 256];
cursor.read_exact(&mut buffer)?;
let (_, beam_config) = parse_beam_config(&buffer).unwrap();
beam_configs.push(beam_config);
}
// 读取 RADIAL_DATA
let (datas, els, azs, ranges) =
parse_final_data(total_length, task_config.cut_number as usize, cursor)?;
// 构建并返回 RadarData 结构体
let radar_data = RadarData {
datetime: task_config.scan_start_time,
els,
azs,
ranges,
datas,
};
Ok(radar_data)
}
fn parse_generic_header(input: &[u8]) -> IResult<&[u8], GenericHeader> {
let (input, magic_number) = le_u32(input)?;
let (input, major_version) = le_u16(input)?;
let (input, minor_version) = le_u16(input)?;
let (input, generic_type) = le_u32(input)?;
let (input, product_type) = le_u32(input)?;
let (input, reserved) = take(16usize)(input)?;
let generic_header = GenericHeader {
magic_number,
major_version,
minor_version,
generic_type,
product_type,
reserved: reserved.try_into().unwrap(),
};
Ok((input, generic_header))
}
fn parse_site_config(input: &[u8]) -> IResult<&[u8], SiteConfig> {
let (input, site_code) = take(8usize)(input)?;
let (input, site_name) = take(32usize)(input)?;
let (input, latitude) = le_f32(input)?;
let (input, longitude) = le_f32(input)?;
let (input, antenna_height) = le_u32(input)?;
let (input, _) = take(4usize)(input)?; // Skip reserved bytes
let site_config = SiteConfig {
site_code: String::from_utf8(site_code.to_vec())
.unwrap()
.trim_end_matches('\0')
.to_string(),
site_name: String::from_utf8(site_name.to_vec())
.unwrap()
.trim_end_matches('\0')
.to_string(),
latitude,
longitude,
antenna_height,
};
Ok((input, site_config))
}
fn parse_task_config(input: &[u8]) -> IResult<&[u8], TaskConfig> {
let (input, task_name) = take(32usize)(input)?;
let (input, task_description) = take(128usize)(input)?;
let (input, polarization_type) = le_u32(input)?;
let (input, scan_type) = le_u32(input)?;
let (input, scan_beam_number) = le_u32(input)?;
let (input, cut_number) = le_u32(input)?;
let (input, ray_order) = le_u32(input)?;
let (input, scan_start_time) = le_i64(input)?;
let (input, _) = take(68usize)(input)?; // Skip reserved bytes
let task_config = TaskConfig {
task_name: String::from_utf8_lossy(task_name).to_string(),
task_description: String::from_utf8_lossy(task_description).to_string(),
polarization_type,
scan_type,
scan_beam_number,
cut_number,
ray_order,
scan_start_time,
};
Ok((input, task_config))
}
fn single_beam_block(raw: &[u8], bin_length: usize) -> io::Result<Vec<f64>> {
let raw = match bin_length {
1 => {
let new_raw: Vec<f64> = raw.into_iter().map(|v| *v as f64).collect();
new_raw
}
2 => {
let new_raw: Vec<f64> = cast_slice(&raw)
.into_iter()
.map(|v: &u16| *v as f64)
.collect();
new_raw
}
4 => {
let new_raw: Vec<f64> = cast_slice(&raw)
.into_iter()
.map(|v: &u32| *v as f64)
.collect();
new_raw
}
8 => {
let new_raw: Vec<f64> = cast_slice(&raw)
.into_iter()
.map(|v: &u64| *v as f64)
.collect();
new_raw
}
_ => {
let new_raw: Vec<f64> = raw.into_iter().map(|v| *v as f64).collect();
new_raw
}
};
Ok(raw)
}
fn beam_block_exact(cursor: &mut Cursor<Vec<u8>>, header: &TypeHeader) -> io::Result<Vec<u8>> {
let mut buf = vec![0u8; header.length as usize];
cursor.read_exact(&mut buf)?;
Ok(buf)
}
fn data_type(typ: u32) -> &'static str {
match typ {
2 => "DBZ",
3 => "VEL",
4 => "SW",
7 => "ZDR",
9 => "CC",
10 => "PHIDP",
11 => "KDP",
_ => "Unknown",
}
}
fn parse_final_data(
all_length: u64,
cut_number: usize,
mut cursor: Cursor<Vec<u8>>,
) -> io::Result<(
HashMap<&'static str, Vec<f64>>,
Vec<f64>,
Vec<f64>,
Vec<f64>,
)> {
let position = cursor.position();
let mut radial_buffer = vec![0u8; 128];
let mut header_buffer = vec![0u8; 32];
cursor.read_exact(&mut radial_buffer)?;
// First Block
let (_, radial_data) = parse_radial_data(&radial_buffer).unwrap();
// Radial_number: I don't know why not use it.
// let radial_number = radial_data.radial_number as usize;
// 变量类型
let type_number = radial_data.moment_number as usize;
let mut gate_len = vec![0; type_number];
let mut types = vec!["Unknown"; type_number];
// final datas
let mut datas = HashMap::new();
for idx in 0..type_number {
cursor.read_exact(&mut header_buffer)?;
let (_, header) = parse_type_header(&header_buffer).unwrap();
let first_beam_block = beam_block_exact(&mut cursor, &header)?;
let first_beam = single_beam_block(&first_beam_block, header.bin_length as usize)?;
gate_len[idx] = first_beam.len();
let data_type = data_type(header.data_type);
if data_type == "Unknown" {
continue;
}
types[idx] = data_type;
}
let first_block_position = cursor.position();
let roll_num = (all_length - position) / (first_block_position - position);
let azm_number = (roll_num / cut_number as u64) as usize;
for idx in 0..type_number {
if types[idx] == "Unknown" {
continue;
}
datas.insert(
types[idx],
vec![0f64; cut_number * azm_number * gate_len[idx]],
);
}
// Reset Cursor
cursor.set_position(position);
let mut els = vec![0f64; cut_number];
let mut azs = vec![0f64; azm_number];
for e in 0..cut_number {
for a in 0..azm_number {
cursor.read_exact(&mut radial_buffer)?;
let (_, radial_data) = parse_radial_data(&radial_buffer).unwrap();
els[e] = radial_data.elevation as f64;
azs[a] = radial_data.azimuth as f64;
let type_number = radial_data.moment_number as usize;
for typ in 0..type_number {
cursor.read_exact(&mut header_buffer)?;
let (_, header) = parse_type_header(&header_buffer).unwrap();
let beam_block = beam_block_exact(&mut cursor, &header)?;
let mut beam = single_beam_block(&beam_block, header.bin_length as usize)?;
beam.iter_mut()
.for_each(|v| *v = (*v - header.offset as f64) / header.scale as f64);
let typ_name = types.get(typ).unwrap();
if typ_name == &"Unknown" {
continue;
}
let e = radial_data.elevation_number as usize - 1;
datas.get_mut(typ_name).unwrap()[e * (azm_number * gate_len[typ])
+ a * gate_len[typ]
..e * (azm_number * gate_len[typ]) + (a + 1) * gate_len[typ]]
.copy_from_slice(&beam);
}
}
}
let mut ranges = vec![0f64; gate_len[0]];
ranges
.iter_mut()
.enumerate()
.for_each(|(idx, r)| *r = idx as f64 * 30.0);
Ok((datas, els, azs, ranges))
}
fn parse_beam_config(input: &[u8]) -> IResult<&[u8], BeamConfig> {
let (input, cut_index) = le_i16(input)?;
let (input, tx_beam_index) = le_i16(input)?;
let (input, rx_beam_elevation) = le_f32(input)?;
let (input, tx_beam_gain) = le_f32(input)?;
let (input, rx_beam_width_h) = le_f32(input)?;
let (input, rx_beam_width_v) = le_f32(input)?;
let (input, rx_beam_gain) = le_f32(input)?;
let (input, process_mode) = map(le_u32, |v| match v {
1 => "PPP".to_string(),
2 => "FFT".to_string(),
_ => "Unknown".to_string(),
})(input)?;
let (input, wave_form) = map(le_u32, |v| match v {
0 => "CS连续监测".to_string(),
1 => "CD连续多普勒".to_string(),
2 => "CDX多普勒扩展".to_string(),
3 => "Rx Test".to_string(),
4 => "BATCH批模式".to_string(),
5 => "Dual PRF双PRF".to_string(),
6 => "Staggered PRT 参差PRT".to_string(),
_ => "Unknown".to_string(),
})(input)?;
let (input, n1_prf1) = le_f32(input)?;
let (input, n1_prf2) = le_f32(input)?;
let (input, n2_prf1) = le_f32(input)?;
let (input, n2_prf2) = le_f32(input)?;
let (input, dealiasing_mode) = le_u32(input)?;
let (input, azimuth) = le_f32(input)?;
let (input, start_angle) = le_f32(input)?;
let (input, end_angle) = le_f32(input)?;
let (input, angular_resolution) = le_f32(input)?;
let (input, scan_speed) = le_f32(input)?;
let (input, log_resolution) = le_f32(input)?;
let (input, doppler_resolution) = le_f32(input)?;
let (input, maximum_range1) = le_u32(input)?;
let (input, maximum_range2) = le_u32(input)?;
let (input, start_range) = le_u32(input)?;
let (input, sample1) = le_u32(input)?;
let (input, sample2) = le_u32(input)?;
let (input, phase_mode) = le_u32(input)?;
let (input, atmospheric_loss) = le_f32(input)?;
let (input, nyquist_speed) = le_f32(input)?;
let (input, moments_mask) = le_i64(input)?;
let (input, moments_size_mask) = le_i64(input)?;
let (input, misc_filter_mask) = le_u32(input)?;
let (input, sqi_threshold) = le_f32(input)?;
let (input, sig_threshold) = le_f32(input)?;
let (input, csr_threshold) = le_f32(input)?;
let (input, log_threshold) = le_f32(input)?;
let (input, cpa_threshold) = le_f32(input)?;
let (input, pmi_threshold) = le_f32(input)?;
let (input, dplog_threshold) = le_f32(input)?;
let (input, thresholds_r) = take(4usize)(input)?;
let (input, dbt_mask) = le_u32(input)?;
let (input, dbz_mask) = le_u32(input)?;
let (input, velocity_mask) = le_u32(input)?;
let (input, spectrum_width_mask) = le_u32(input)?;
let (input, dp_mask) = le_u32(input)?;
let (input, direction) = le_u32(input)?;
let (input, ground_clutter_classifier_type) = le_i16(input)?;
let (input, ground_clutter_filter_type) = le_i16(input)?;
let (input, ground_clutter_filter_notch_width) = le_i16(input)?;
let (input, ground_clutter_filter_window) = le_i16(input)?;
let (input, _) = take(44usize)(input)?;
let beam_config = BeamConfig {
cut_index,
tx_beam_index,
rx_beam_elevation,
tx_beam_gain,
rx_beam_width_h,
rx_beam_width_v,
rx_beam_gain,
process_mode,
wave_form,
n1_prf1,
n1_prf2,
n2_prf1,
n2_prf2,
dealiasing_mode,
azimuth,
start_angle,
end_angle,
angular_resolution,
scan_speed,
log_resolution,
doppler_resolution,
maximum_range1,
maximum_range2,
start_range,
sample1,
sample2,
phase_mode,
atmospheric_loss,
nyquist_speed,
moments_mask,
moments_size_mask,
misc_filter_mask,
sqi_threshold,
sig_threshold,
csr_threshold,
log_threshold,
cpa_threshold,
pmi_threshold,
dplog_threshold,
thresholds_r: thresholds_r.try_into().unwrap(),
dbt_mask,
dbz_mask,
velocity_mask,
spectrum_width_mask,
dp_mask,
direction,
ground_clutter_classifier_type,
ground_clutter_filter_type,
ground_clutter_filter_notch_width,
ground_clutter_filter_window,
};
Ok((input, beam_config))
}
fn parse_radial_data(input: &[u8]) -> IResult<&[u8], RadialData> {
let (input, radial_state) = le_u32(input)?;
let (input, spot_blank) = le_u32(input)?;
let (input, sequence_number) = le_u32(input)?;
let (input, radial_number) = le_u32(input)?;
let (input, elevation_number) = le_u32(input)?;
let (input, azimuth) = le_f32(input)?;
let (input, elevation) = le_f32(input)?;
let (input, seconds) = le_i64(input)?;
let (input, microseconds) = le_u32(input)?;
let (input, length_of_data) = le_u32(input)?;
let (input, moment_number) = le_u32(input)?;
let (input, scan_beam_index) = le_i16(input)?;
let (input, horizontal_estimated_noise) = le_i16(input)?;
let (input, vertical_estimated_noise) = le_i16(input)?;
let (input, prf_flag) = le_u32(input)?;
let (input, _) = take(70usize)(input)?; // Skip reserved bytes
let radial_data = RadialData {
radial_state,
spot_blank,
sequence_number,
radial_number,
elevation_number,
azimuth,
elevation,
seconds,
microseconds,
length_of_data,
moment_number,
scan_beam_index,
horizontal_estimated_noise,
vertical_estimated_noise,
prf_flag,
};
Ok((input, radial_data))
}
fn parse_type_header(input: &[u8]) -> IResult<&[u8], TypeHeader> {
let (input, data_type) = le_u32(input)?;
let (input, scale) = le_u32(input)?;
let (input, offset) = le_u32(input)?;
let (input, bin_length) = le_u16(input)?;
let (input, flags) = le_u16(input)?;
let (input, length) = le_u32(input)?;
let (input, _) = take(12usize)(input)?; // Skip reserved bytes
let type_header = TypeHeader {
data_type,
scale,
offset,
bin_length,
flags,
length,
};
Ok((input, type_header))
}
mod test {
use crate::raw::parse_raw_data;
use super::get_radar_data;
#[test]
fn test() {
use std::fs::File;
let radar_data = parse_raw_data(
"/Users/tsuki/Desktop/Z_RADR_I_X5775_20230726180000_O_DOR-XPD-CAP-FMT.BIN.zip",
)
.unwrap();
}
}