kopia lustrzana https://github.com/cyoung/stratux
316 wiersze
12 KiB
C
316 wiersze
12 KiB
C
//
|
|
// Copyright 2015, Oliver Jowett <oliver@mutability.co.uk>
|
|
//
|
|
|
|
// This file is free software: you may copy, redistribute and/or modify it
|
|
// under the terms of the GNU General Public License as published by the
|
|
// Free Software Foundation, either version 2 of the License, or (at your
|
|
// option) any later version.
|
|
//
|
|
// This file is distributed in the hope that it will be useful, but
|
|
// WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
// General Public License for more details.
|
|
//
|
|
// You should have received a copy of the GNU General Public License
|
|
// along with this program. If not, see <http://www.gnu.org/licenses/>.
|
|
|
|
#include <stdio.h>
|
|
#include <math.h>
|
|
|
|
#include "uat.h"
|
|
#include "uat_decode.h"
|
|
#include "reader.h"
|
|
|
|
#define BLOCK_WIDTH (48.0/60.0)
|
|
#define WIDE_BLOCK_WIDTH (96.0/60.0)
|
|
#define BLOCK_HEIGHT (4.0/60.0)
|
|
#define BLOCK_THRESHOLD 405000
|
|
#define BLOCKS_PER_RING 450
|
|
|
|
//
|
|
// This reads demodulated uplink messages and extracts NEXRAD global block representation formats - type 63 and 64
|
|
//
|
|
// The output format is a series of lines, one line per decoded block.
|
|
// Each line is space-separated and is formatted as:
|
|
//
|
|
// NEXRAD <type> <hour>:<minute> <scale> <north> <west> <height> <width> <data>
|
|
//
|
|
// where:
|
|
// <type> is Regional (for type 63) or CONUS (for type 64)
|
|
// <hour>:<minute> is the time from the PDU header - all blocks from one composite radar image will have the same time
|
|
// <scale> is the scale value of this block - 0 (high res), 1 (med res), or 2 (low res)
|
|
// <north> is the north edge of the block, in _integer arcminutes_. Divide by 60 to get degrees.
|
|
// <west> is the west edge of the block, in _positive integer arcminutes_. Divide by 60 to get degrees; subtract 360 if you want the conventional -180..+180 range.
|
|
// <height> is the height of the block, in integer arcminutes of latitude
|
|
// <width> is the width of the block, in integer arcminutes of longitude
|
|
//
|
|
// Each block contains 128 evenly spaced bins, in a grid of 32 (longitude) x 4 (latitude), working west-to-east then north-to-south.
|
|
// i.e. each bin represents a pixel that covers <width>/32 arcminutes of longitude by <height>/4 arcminutes of latitude.
|
|
//
|
|
// <data> is a string of 128 digits (no spaces); each character represents the intensity of one bin, in the order above.
|
|
//
|
|
|
|
// Given:
|
|
// bn - block number
|
|
// ns - north/south flag
|
|
// sf - scale factor
|
|
//
|
|
// compute the northwest corner of the referenced block and place it in (*latN, *lonW)
|
|
// and place the total height and width of the block in (*latSize, *lonSize)
|
|
void block_location(int bn, int ns, int sf, double *latN, double *lonW, double *latSize, double *lonSize)
|
|
{
|
|
// With SF=0:
|
|
|
|
// blocks are (48 arcminutes longitude) x (4 arcminute latitude) between 0 and 60 degrees latitude
|
|
// (450 blocks for each ring of latitude)
|
|
// blocks are (96 arcminutes longitude) x (4 arcminute latitude) between 60 and 90 degrees latitude
|
|
// (225 blocks for each ring of latitude) - but the block numbering continues to use
|
|
// a 48-arcminute spacing, so only even numbered blocks are meaningful.
|
|
// block zero is immediately northeast of (0,0), then blocks are numbered east-to-west, south-to-north.
|
|
//
|
|
// Southern hemisphere numbering is mirrored around the equator, and indicated by the "ns" flag.
|
|
|
|
// ^N
|
|
// | 405446 | 405448 | 405000 | 405002 |
|
|
// --------------------------------------------------------- 60 00' 00" N
|
|
// |404996|404997|404998|404999|404550|404551|404552|404553|
|
|
// --------------------------------------------------------- 59 56' 00" N
|
|
// ...
|
|
// | 896 | 897 | 898 | 899 | 450 | 451 | 452 | 453 |
|
|
// --------------------------------------------------------- 00 04' 00" N
|
|
// | 446 | 447 | 447 | 449 | 0 | 1 | 2 | 3 |
|
|
//W<------------------------------------------------------->E equator
|
|
// | 446* | 447* | 447* | 449* | 0* | 1* | 2* | 3* |
|
|
// --------------------------------------------------------- 00 04' 00" S
|
|
// | 896* | 897* | 898* | 899* | 450* | 451* | 452* | 453* |
|
|
// 2d24'W 1d36'W 0d48'W V 0d48'E 1d36'E 2d24'E
|
|
// (* = ns_flag set)
|
|
|
|
// Each block is subdivided into 32 (longitude) x 4 (latitude) bins.
|
|
// The bins are numbered starting at the northwest corner of the block,
|
|
// west-to-east then north-to-south.
|
|
|
|
// block 0:
|
|
//
|
|
// ------------------------------------ <- 0d04m00s N
|
|
// | 0 1 2 3 ... 28 29 30 31| <- each bin is 1 arcminute tall
|
|
// | 32 33 34 35 ... 60 61 62 63|
|
|
// | 64 65 66 67 ... 92 93 94 95|
|
|
// | 96 97 98 99 ... 124 125 126 127|
|
|
// ------------------------------------ <- 0N - equator
|
|
// ^ ^ each bin is ^
|
|
// 0E 1.5 arcminute wide 0d48m00s E
|
|
|
|
// With SF=1, an identical block numbering is used to locate the northwest corner of the block,
|
|
// but then each bin is 5x larger in both axes i.e. 240 x 20 or 480 x 20 arcminutes.
|
|
// this means that the block data will actually overlap 24 other block positions.
|
|
|
|
// With SF=2, it works like SF=1 but with a scale factor of 9x.
|
|
|
|
double raw_lat, raw_lon;
|
|
double scale;
|
|
|
|
if (sf == 1)
|
|
scale = 5.0;
|
|
else if (sf == 2)
|
|
scale = 9.0;
|
|
else
|
|
scale = 1.0;
|
|
|
|
if (bn >= BLOCK_THRESHOLD) {
|
|
// 60-90 degrees - even-numbered blocks only
|
|
bn = bn & ~1;
|
|
}
|
|
|
|
raw_lat = BLOCK_HEIGHT * trunc(bn / BLOCKS_PER_RING);
|
|
raw_lon = (bn % BLOCKS_PER_RING) * BLOCK_WIDTH;
|
|
|
|
*lonSize = (bn >= BLOCK_THRESHOLD ? WIDE_BLOCK_WIDTH : BLOCK_WIDTH) * scale;
|
|
*latSize = BLOCK_HEIGHT * scale;
|
|
|
|
// raw_lat/raw_lon points to the southwest corner in the northern hemisphere version
|
|
*lonW = raw_lon;
|
|
if (ns) {
|
|
// southern hemisphere, mirror along the equator
|
|
*latN = 0 - raw_lat;
|
|
} else {
|
|
// adjust to the northwest corner
|
|
*latN = raw_lat + BLOCK_HEIGHT;
|
|
}
|
|
}
|
|
|
|
void decode_nexrad(struct fisb_apdu *fisb)
|
|
{
|
|
// Header:
|
|
//
|
|
// byte/bit 7 6 5 4 3 2 1 0
|
|
// 0 |RLE|NS | Scale | MSB Block # |
|
|
// 1 | Block # |
|
|
// 2 | Block # LSB |
|
|
|
|
|
|
int rle_flag = (fisb->data[0] & 0x80) != 0;
|
|
int ns_flag = (fisb->data[0] & 0x40) != 0;
|
|
int block_num = ((fisb->data[0] & 0x0f) << 16) | (fisb->data[1] << 8) | (fisb->data[2]);
|
|
int scale_factor = (fisb->data[0] & 0x30) >> 4;
|
|
|
|
// now decode the bins
|
|
if (rle_flag) {
|
|
// One bin, 128 values, RLE-encoded
|
|
int i;
|
|
double latN = 0, lonW = 0, latSize = 0, lonSize = 0;
|
|
block_location(block_num, ns_flag, scale_factor, &latN, &lonW, &latSize, &lonSize);
|
|
|
|
fprintf(stdout, "NEXRAD %s %02d:%02d %d %.0f %.0f %.0f %.0f ",
|
|
fisb->product_id == 63 ? "Regional" : "CONUS",
|
|
fisb->hours,
|
|
fisb->minutes,
|
|
scale_factor,
|
|
latN * 60,
|
|
lonW * 60,
|
|
latSize * 60,
|
|
lonSize * 60);
|
|
|
|
// each byte following the header is:
|
|
// 7 6 5 4 3 2 1 0
|
|
// | runlength - 1 | intensity |
|
|
|
|
for (i = 3; i < fisb->length; ++i) {
|
|
int intensity = fisb->data[i] & 7;
|
|
int runlength = (fisb->data[i] >> 3) + 1;
|
|
|
|
while (runlength-- > 0)
|
|
fprintf(stdout, "%d", intensity);
|
|
}
|
|
fprintf(stdout, "\n");
|
|
} else {
|
|
int L = fisb->data[3] & 15;
|
|
int i;
|
|
int row_start, row_offset, row_size;
|
|
|
|
//
|
|
// Empty block representation, representing one
|
|
// or more blocks that are completely empty of
|
|
// data.
|
|
//
|
|
// 7 6 5 4 3 2 1 0
|
|
// 3 |b+4 |b+3 |b+2 |b+1 | length (L) |
|
|
// 4 |b+12|b+11|b+10|b+9 |b+8 |b+7 |b+6 |b+5 |
|
|
// ...
|
|
// 3+L |b+8L-3 ... b+8L+4|
|
|
|
|
// The block number from the header is always
|
|
// empty.
|
|
//
|
|
// If the bit for b+x is empty, then the block
|
|
// X to the right of the block from the header is
|
|
// empty. Note that the block is _always on the
|
|
// same row_ even if the offset would make the
|
|
// block cross the 0E meridian, so it is not simply
|
|
// a case of adding to the block number.
|
|
|
|
// find the lowest-numbered block of this row
|
|
if (block_num >= 405000) {
|
|
row_start = block_num - ((block_num - 405000) % 225);
|
|
row_size = 225;
|
|
} else {
|
|
row_start = block_num - (block_num % 450);
|
|
row_size = 450;
|
|
}
|
|
|
|
// find the offset of the first block in this row handled
|
|
// by this message
|
|
row_offset = block_num - row_start;
|
|
|
|
for (i = 0; i < L; ++i) {
|
|
int bb;
|
|
int j;
|
|
|
|
if (i == 0)
|
|
bb = (fisb->data[3] & 0xF0) | 0x08; // synthesize a first byte in the same format as all the other bytes
|
|
else
|
|
bb = (fisb->data[i+3]);
|
|
|
|
for (j = 0; j < 8; ++j) {
|
|
if (bb & (1 << j)) {
|
|
// find the relevant block for this bit, limited
|
|
// to the same row as the original block.
|
|
int row_x = (row_offset + 8*i + j - 3) % row_size;
|
|
int bn = row_start + row_x;
|
|
double latN = 0, lonW = 0, latSize = 0, lonSize = 0;
|
|
int k;
|
|
block_location(bn, ns_flag, scale_factor, &latN, &lonW, &latSize, &lonSize);
|
|
|
|
fprintf(stdout, "NEXRAD %s %02d:%02d %d %.0f %.0f %.0f %.0f ",
|
|
fisb->product_id == 63 ? "Regional" : "CONUS",
|
|
fisb->hours,
|
|
fisb->minutes,
|
|
scale_factor,
|
|
latN * 60,
|
|
lonW * 60,
|
|
latSize * 60,
|
|
lonSize * 60);
|
|
|
|
// seems to work best if we assume that
|
|
// CONUS empty blocks = intensity 1 (valid data, but no precipitation)
|
|
// regional empty blocks = intensity 0 (valid data <5dBz)
|
|
for (k = 0; k < 128; ++k)
|
|
fprintf(stdout, "%d", (fisb->product_id == 63 ? 0 : 1));
|
|
fprintf(stdout, "\n");
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
void handle_frame(frame_type_t type, uint8_t *frame, int len, void *extra)
|
|
{
|
|
if (type == UAT_UPLINK) {
|
|
struct uat_uplink_mdb mdb;
|
|
int i;
|
|
|
|
uat_decode_uplink_mdb(frame, &mdb);
|
|
if (!mdb.app_data_valid)
|
|
return;
|
|
|
|
for (i = 0; i < mdb.num_info_frames; ++i) {
|
|
struct fisb_apdu *fisb;
|
|
|
|
if (!mdb.info_frames[i].is_fisb)
|
|
continue;
|
|
|
|
fisb = &mdb.info_frames[i].fisb;
|
|
if (fisb->product_id != 63 && fisb->product_id != 64)
|
|
continue;
|
|
|
|
decode_nexrad(fisb);
|
|
}
|
|
}
|
|
|
|
fflush(stdout);
|
|
}
|
|
|
|
int main(int argc, char **argv)
|
|
{
|
|
struct dump978_reader *reader;
|
|
int framecount;
|
|
|
|
reader = dump978_reader_new(0,0);
|
|
if (!reader) {
|
|
perror("dump978_reader_new");
|
|
return 1;
|
|
}
|
|
|
|
while ((framecount = dump978_read_frames(reader, handle_frame, NULL)) > 0)
|
|
;
|
|
|
|
if (framecount < 0) {
|
|
perror("dump978_read_frames");
|
|
return 1;
|
|
}
|
|
|
|
return 0;
|
|
}
|
|
|