Aqua AIRS and AMSU projection

This commit is contained in:
Aang23 2021-10-22 18:02:59 +02:00
parent 699076e888
commit 07b52f5327
16 changed files with 336 additions and 40 deletions

View file

@ -139,10 +139,7 @@
"aqua_db_decoder": {}
},
"products": {
"eos_modis": {
"terra_mode": "0",
"correct_bowtie": "1"
},
"aqua_airs": {},
"aqua_ceres": {},
"aqua_amsu": {}

View file

@ -1,3 +1,3 @@
AQUA
1 27424U 02022A 21244.48848139 .00000150 00000-0 43339-4 0 9995
2 27424 98.2121 184.7351 0002252 69.2846 78.9339 14.57107061 28101
1 27424U 02022A 21295.15021017 .00000141 00000-0 41449-4 0 9993
2 27424 98.2196 234.7086 0000549 84.0205 338.1202 14.57091753 35529

View file

@ -27,4 +27,13 @@ namespace ccsds
return double(offset + days) * 86400.0 + double(milliseconds_of_day) / double(ms_scale) + double(microseconds_of_millisecond) / double(us_of_ms_scale);
}
double parseCCSDSTimeFullRawUnsegmented(uint8_t *data, int offset, double ms_scale)
{
uint8_t seconds_to_convert = data[1] & 0b1111111;
uint32_t seconds_since_epoch = data[2] << 24 | data[3] << 16 | data[4] << 8 | data[5];
uint16_t subsecond_time = data[6] << 8 | data[7];
return offset * 86400.0 + double(seconds_since_epoch - seconds_to_convert) + double(subsecond_time) * double(ms_scale);
}
}

View file

@ -5,7 +5,8 @@
namespace ccsds
{
time_t parseCCSDSTime(CCSDSPacket &pkt, int offset = 0, int ms_scale = 1000); // Parse timestamp with 1 second accuracy
double parseCCSDSTimeFull(CCSDSPacket &pkt, int offset, int ms_scale = 1000, int ns_of_ms_scale = 1000); // Parse timestamp with full accuracy
time_t parseCCSDSTime(CCSDSPacket &pkt, int offset = 0, int ms_scale = 1000); // Parse timestamp with 1 second accuracy
double parseCCSDSTimeFull(CCSDSPacket &pkt, int offset, int ms_scale = 1000, int ns_of_ms_scale = 1000); // Parse timestamp with full accuracy
double parseCCSDSTimeFullRaw(uint8_t *data, int offset, int ms_scale = 1000, int ns_of_ms_scale = 1000); // Parse timestamp with full accuracy, from a raw buffer
double parseCCSDSTimeFullRawUnsegmented(uint8_t *data, int offset, double ms_scale); // Parse Coarse timestamp. Used by some EOS instruments on Aqua
}

View file

@ -55,6 +55,8 @@ namespace geodetic
projected_image.draw_circle(map_cc1.first, map_cc1.second, ceil(circle_radius), color, 0.4 * opacity);
}
//projected_image.draw_point(map_cc1.first, map_cc1.second, color, opacity);
if (progress != nullptr)
*progress = float(currentScan) / float(image.height());
}

View file

@ -86,5 +86,85 @@ namespace image
template cimg_library::CImg<unsigned char> correctGenericBowTie(cimg_library::CImg<unsigned char> &, const int, const long, const float, const float);
template cimg_library::CImg<unsigned short> correctGenericBowTie(cimg_library::CImg<unsigned short> &, const int, const long, const float, const float);
/*
template <typename T>
cimg_library::CImg<T> correctSingleBowTie(cimg_library::CImg<T> &inputImage, const int channelCount, const long scanHeight, const float alpha, const float beta)
{
// Compute everything we'll need
const long height = inputImage.height();
const long width = inputImage.width();
const long scanCount = height / scanHeight;
const long halfWidth = width / 2;
// Create our output image
cimg_library::CImg<T> outputImage = cimg_library::CImg<T>(width, height, 1, channelCount);
// Reserve our buffers
T *scan_buffer_input = new T[height * width];
T *scan_buffer_output = new T[height * width];
T *col_buffer_input = new T[scanHeight];
T *col_buffer_output = new T[scanHeight];
for (int channel = 0; channel < channelCount; channel++)
{
for (int scanNumber = 0; scanNumber < scanCount; scanNumber++)
{
//std::cout << "Processing scan " << scanNumber << std::endl;
// Load out input buffer
for (int lineNumber = 0; lineNumber < scanHeight; lineNumber++)
{
for (int pixelNumber = 0; pixelNumber < width; pixelNumber++)
{
scan_buffer_input[lineNumber * width + pixelNumber] = inputImage[(channel * width * height) + ((scanNumber * scanHeight) + lineNumber) * width + pixelNumber];
}
}
for (int rowNumber = 0; rowNumber < width; rowNumber++)
{
// Load our column
for (int lineNumber = 0; lineNumber < scanHeight; lineNumber++)
{
col_buffer_input[lineNumber] = scan_buffer_input[lineNumber * width + rowNumber];
}
int topPixelCount = int((((halfWidth - abs(rowNumber - halfWidth)) / (float)halfWidth) * alpha + beta) * scanHeight);
topPixelCount = std::min<int>(topPixelCount, scanHeight);
for (int i = 0; i < scanHeight; i++)
{
col_buffer_output[i] = col_buffer_input[int(((float)i / (float)scanHeight) * topPixelCount)];
}
// Offload our columm
for (int lineNumber = 0; lineNumber < scanHeight; lineNumber++)
{
scan_buffer_output[lineNumber * width + rowNumber] = col_buffer_output[lineNumber];
}
}
// Offload out output buffer
for (int lineNumber = 0; lineNumber < scanHeight; lineNumber++)
{
for (int pixelNumber = 0; pixelNumber < width; pixelNumber++)
{
outputImage[(channel * width * height) + ((scanNumber * scanHeight) + lineNumber) * width + pixelNumber] = scan_buffer_output[lineNumber * width + pixelNumber];
}
}
}
}
delete[] scan_buffer_input;
delete[] scan_buffer_output;
delete[] col_buffer_input;
delete[] col_buffer_output;
return outputImage;
}
template cimg_library::CImg<unsigned char> correctSingleBowTie(cimg_library::CImg<unsigned char> &, const int, const long, const float, const float);
template cimg_library::CImg<unsigned short> correctSingleBowTie(cimg_library::CImg<unsigned short> &, const int, const long, const float, const float);
*/
}
}

View file

@ -10,5 +10,8 @@ namespace image
{
template <typename T>
cimg_library::CImg<T> correctGenericBowTie(cimg_library::CImg<T> &inputImage, const int channelCount, const long scanHeight, const float alpha, const float beta);
//template <typename T>
//cimg_library::CImg<T> correctSingleBowTie(cimg_library::CImg<T> &inputImage, const int channelCount, const long scanHeight, const float alpha, const float beta);
}
}

View file

@ -1,5 +1,6 @@
#include "airs_reader.h"
#include "utils.h"
#include "common/ccsds/ccsds_time.h"
namespace aqua
{
@ -16,6 +17,8 @@ namespace aqua
hd_channels[i] = new unsigned short[10000 * 90 * 9];
}
lines = 0;
timestamps_ifov.push_back(std::vector<double>(90));
std::fill(timestamps_ifov[lines].begin(), timestamps_ifov[lines].end(), -1);
}
AIRSReader::~AIRSReader()
@ -88,9 +91,17 @@ namespace aqua
}
}
// Timestamp
double timestamp = ccsds::parseCCSDSTimeFullRawUnsegmented(&packet.payload[1], -4383, 15.3e-6);
timestamps_ifov[lines][pix_pos] = timestamp;
// Frame counter
if (counter == 22 || counter == 278 || counter == 534)
{
lines++;
timestamps_ifov.push_back(std::vector<double>(90));
std::fill(timestamps_ifov[lines].begin(), timestamps_ifov[lines].end(), -1);
}
}
cimg_library::CImg<unsigned short> AIRSReader::getChannel(int channel)

View file

@ -20,6 +20,7 @@ namespace aqua
AIRSReader();
~AIRSReader();
int lines;
std::vector<std::vector<double>> timestamps_ifov;
void work(ccsds::CCSDSPacket &packet);
cimg_library::CImg<unsigned short> getChannel(int channel);
cimg_library::CImg<unsigned short> getHDChannel(int channel);

View file

@ -6,6 +6,9 @@
#include "logger.h"
#include <filesystem>
#include "imgui/imgui.h"
#include "common/geodetic/projection/satellite_reprojector.h"
#include "common/geodetic/projection/proj_file.h"
#include "modules/eos/eos.h"
#define BUFFER_SIZE 8192
@ -101,6 +104,16 @@ namespace aqua
WRITE_IMAGE(airs_reader.getHDChannel(i), directory + "/AIRS-HD-" + std::to_string(i + 1) + ".png");
}
// There nearly 3000 channels... So we write that in a specific folder not to fill up the main one
if (!std::filesystem::exists(directory + "/Channels"))
std::filesystem::create_directory(directory + "/Channels");
for (int i = 0; i < 2666; i++)
{
logger->info("Channel " + std::to_string(i + 1) + "...");
WRITE_IMAGE(airs_reader.getChannel(i), directory + "/Channels/AIRS-" + std::to_string(i + 1) + ".png");
}
// Output a few nice composites as well
logger->info("Global Composite...");
int all_width_count = 100;
@ -148,6 +161,74 @@ namespace aqua
image321.draw_image(0, 0, 0, 2, airs_reader.getHDChannel(0));
}
WRITE_IMAGE(image321, directory + "/AIRS-HD-RGB-321.png");
// Normal res projecition
if (airs_reader.lines > 0)
{
int norad = EOS_AQUA_NORAD;
// Setup Projecition
std::shared_ptr<geodetic::projection::LEOScanProjectorSettings_IFOV> proj_settings = std::make_shared<geodetic::projection::LEOScanProjectorSettings_IFOV>(
98.8, // Scan angle
1.1, // IFOV X scan angle
1.4, // IFOV Y scan angle
-0.5, // Roll offset
0.0, // Pitch offset
-2.5, // Yaw offset
-1.1, // Time offset
90, // Number of IFOVs
1, // IFOV Width
1, // IFOV Height
airs_reader.getChannel(0).width(), // Image width
true, // Invert scan
tle::getTLEfromNORAD(norad), // TLEs
airs_reader.timestamps_ifov // Timestamps
);
geodetic::projection::LEOScanProjector projector(proj_settings);
{
geodetic::projection::proj_file::LEO_GeodeticReferenceFile geofile = geodetic::projection::proj_file::leoRefFileFromProjector(norad, proj_settings);
geodetic::projection::proj_file::writeReferenceFile(geofile, directory + "/AIRS.georef");
}
logger->info("Projected 62 channel...");
cimg_library::CImg<unsigned char> projected_image = geodetic::projection::projectLEOToEquirectangularMapped(airs_reader.getChannel(62), projector, 2048 * 4, 1024 * 4, 1);
WRITE_IMAGE(projected_image, directory + "/AIRS-62-PROJ.png");
}
// HD Projection
if (airs_reader.lines > 0)
{
int norad = EOS_AQUA_NORAD;
// Setup Projecition
std::shared_ptr<geodetic::projection::LEOScanProjectorSettings_IFOV> proj_settings = std::make_shared<geodetic::projection::LEOScanProjectorSettings_IFOV>(
98.8, // Scan angle
1.1, // IFOV X scan angle
1.4, // IFOV Y scan angle
-0.5, // Roll offset
0.0, // Pitch offset
-2.5, // Yaw offset
-1.1, // Time offset
90, // Number of IFOVs
8, // IFOV Width
9, // IFOV Height
airs_reader.getHDChannel(0).width(), // Image width
true, // Invert scan
tle::getTLEfromNORAD(norad), // TLEs
airs_reader.timestamps_ifov // Timestamps
);
geodetic::projection::LEOScanProjector projector(proj_settings);
{
geodetic::projection::proj_file::LEO_GeodeticReferenceFile geofile = geodetic::projection::proj_file::leoRefFileFromProjector(norad, proj_settings);
geodetic::projection::proj_file::writeReferenceFile(geofile, directory + "/AIRS-HD.georef");
}
logger->info("Projected HD RGB 321 channel...");
cimg_library::CImg<unsigned char> projected_image = geodetic::projection::projectLEOToEquirectangularMapped(image321, projector, 2048 * 4, 1024 * 4, 3);
WRITE_IMAGE(projected_image, directory + "/AIRS-HD-RGB-321-PROJ.png");
}
}
void AquaAIRSDecoderModule::drawUI(bool window)

View file

@ -1,4 +1,5 @@
#include "amsu_a1_reader.h"
#include "common/ccsds/ccsds_time.h"
namespace aqua
{
@ -44,6 +45,9 @@ namespace aqua
}
}
double timestamp = ccsds::parseCCSDSTimeFullRawUnsegmented(&packet.payload[1], -4383, 15.3e-6);
timestamps.push_back(timestamp);
// Frame counter
lines++;
}

View file

@ -20,6 +20,7 @@ namespace aqua
AMSUA1Reader();
~AMSUA1Reader();
int lines;
std::vector<double> timestamps;
void work(ccsds::CCSDSPacket &packet);
cimg_library::CImg<unsigned short> getChannel(int channel);
};

View file

@ -1,4 +1,5 @@
#include "amsu_a2_reader.h"
#include "common/ccsds/ccsds_time.h"
namespace aqua
{
@ -41,6 +42,9 @@ namespace aqua
}
}
double timestamp = ccsds::parseCCSDSTimeFullRawUnsegmented(&packet.payload[1], -4383, 15.3e-6);
timestamps.push_back(timestamp);
// Frame counter
lines++;
}

View file

@ -20,6 +20,7 @@ namespace aqua
AMSUA2Reader();
~AMSUA2Reader();
int lines;
std::vector<double> timestamps;
void work(ccsds::CCSDSPacket &packet);
cimg_library::CImg<unsigned short> getChannel(int channel);
};

View file

@ -7,6 +7,9 @@
#include "logger.h"
#include <filesystem>
#include "imgui/imgui.h"
#include "common/geodetic/projection/satellite_reprojector.h"
#include "common/geodetic/projection/proj_file.h"
#include "modules/eos/eos.h"
#define BUFFER_SIZE 8192
@ -156,6 +159,84 @@ namespace aqua
imageAll.draw_image(30 * 6, height, 0, 0, a1reader.getChannel(12));
}
WRITE_IMAGE(imageAll, directory + "/AMSU-ALL.png");
// Reproject to an equirectangular proj
if (a1reader.lines > 0 || a2reader.lines > 0)
{
// Get satellite info
int norad = EOS_AQUA_NORAD;
// Setup Projecition. Twice with the same parameters except we need different timestamps
// There is no "real" guarantee the A1 / A2 output will always be identical
std::shared_ptr<geodetic::projection::LEOScanProjectorSettings_SCANLINE> proj_settings_a1 =
a1reader.lines > 0 ? std::make_shared<geodetic::projection::LEOScanProjectorSettings_SCANLINE>(
98, // Scan angle
-5, // Roll offset
0, // Pitch offset
0, // Yaw offset
10, // Time offset
a1reader.getChannel(0).width(), // Image width
true, // Invert scan
tle::getTLEfromNORAD(norad), // TLEs
a1reader.timestamps // Timestamps
)
: nullptr;
std::shared_ptr<geodetic::projection::LEOScanProjectorSettings_SCANLINE> proj_settings_a2 =
a2reader.lines > 0 ? std::make_shared<geodetic::projection::LEOScanProjectorSettings_SCANLINE>(
98, // Scan angle
-5, // Roll offset
0, // Pitch offset
0, // Yaw offset
10, // Time offset
a2reader.getChannel(0).width(), // Image width
true, // Invert scan
tle::getTLEfromNORAD(norad), // TLEs
a2reader.timestamps // Timestamps
)
: nullptr;
if (a1reader.lines > 0)
{
geodetic::projection::LEOScanProjector projector_a1(proj_settings_a1);
{
geodetic::projection::proj_file::LEO_GeodeticReferenceFile geofile_a1 = geodetic::projection::proj_file::leoRefFileFromProjector(norad, proj_settings_a1);
geodetic::projection::proj_file::writeReferenceFile(geofile_a1, directory + "/AMSU-A1.georef");
}
for (int i = 0; i < 13; i++)
{
cimg_library::CImg<unsigned short> image = a1reader.getChannel(i).equalize(1000).normalize(0, 65535);
image.equalize(1000);
image.normalize(0, 65535);
logger->info("Projected channel A1 " + std::to_string(i + 3) + "...");
cimg_library::CImg<unsigned char> projected_image = geodetic::projection::projectLEOToEquirectangularMapped(image, projector_a1, 1024, 512);
WRITE_IMAGE(projected_image, directory + "/AMSU-A1-" + std::to_string(i + 3) + "-PROJ.png");
}
}
if (a2reader.lines > 0)
{
geodetic::projection::LEOScanProjector projector_a2(proj_settings_a2);
{
geodetic::projection::proj_file::LEO_GeodeticReferenceFile geofile_a1 = geodetic::projection::proj_file::leoRefFileFromProjector(norad, proj_settings_a1);
geodetic::projection::proj_file::writeReferenceFile(geofile_a1, directory + "/AMSU-A1.georef");
geodetic::projection::proj_file::LEO_GeodeticReferenceFile geofile_a2 = geodetic::projection::proj_file::leoRefFileFromProjector(norad, proj_settings_a2);
geodetic::projection::proj_file::writeReferenceFile(geofile_a2, directory + "/AMSU-A2.georef");
}
for (int i = 0; i < 2; i++)
{
cimg_library::CImg<unsigned short> image = a2reader.getChannel(i).equalize(1000).normalize(0, 65535);
image.equalize(1000);
image.normalize(0, 65535);
logger->info("Projected channel A2 " + std::to_string(i + 1) + "...");
cimg_library::CImg<unsigned char> projected_image = geodetic::projection::projectLEOToEquirectangularMapped(image, projector_a2, 1024, 512);
WRITE_IMAGE(projected_image, directory + "/AMSU-A2-" + std::to_string(i + 1) + "-PROJ.png");
}
}
}
}
void AquaAMSUDecoderModule::drawUI(bool window)

View file

@ -13,6 +13,7 @@
#include "modules/eos/eos.h"
#include "nlohmann/json_utils.h"
#include "common/geodetic/projection/satellite_reprojector.h"
#include "common/geodetic/projection/proj_file.h"
#define BUFFER_SIZE 8192
@ -199,6 +200,53 @@ namespace eos
const long scanHeight_500 = 20;
const long scanHeight_1000 = 10;
// Setup GEO projectors
std::shared_ptr<geodetic::projection::LEOScanProjectorSettings_SCANLINE> proj_settings_1000 = std::make_shared<geodetic::projection::LEOScanProjectorSettings_SCANLINE>(
109.9, // Scan angle
-0.0, // Roll offset
0, // Pitch offset
-3.0, // Yaw offset
-3.5, // Time offset
1354, // Image width
true, // Invert scan
tle::getTLEfromNORAD(norad), // TLEs
reader.timestamps_1000 // Timestamps
);
geodetic::projection::LEOScanProjector projector_1000(proj_settings_1000);
std::shared_ptr<geodetic::projection::LEOScanProjectorSettings_SCANLINE> proj_settings_500 = std::make_shared<geodetic::projection::LEOScanProjectorSettings_SCANLINE>(
109.9, // Scan angle
-0.0, // Roll offset
0, // Pitch offset
-3.0, // Yaw offset
-3.5, // Time offset
1354 * 2, // Image width
true, // Invert scan
tle::getTLEfromNORAD(norad), // TLEs
reader.timestamps_500 // Timestamps
);
geodetic::projection::LEOScanProjector projector_500(proj_settings_500);
std::shared_ptr<geodetic::projection::LEOScanProjectorSettings_SCANLINE> proj_settings_250 = std::make_shared<geodetic::projection::LEOScanProjectorSettings_SCANLINE>(
109.9, // Scan angle
-0.0, // Roll offset
0, // Pitch offset
-3.0, // Yaw offset
-3.5, // Time offset
1354 * 4, // Image width
true, // Invert scan
tle::getTLEfromNORAD(norad), // TLEs
reader.timestamps_250 // Timestamps
);
geodetic::projection::LEOScanProjector projector_250(proj_settings_250);
{
geodetic::projection::proj_file::LEO_GeodeticReferenceFile geofile_1000 = geodetic::projection::proj_file::leoRefFileFromProjector(norad, proj_settings_1000);
geodetic::projection::proj_file::LEO_GeodeticReferenceFile geofile_500 = geodetic::projection::proj_file::leoRefFileFromProjector(norad, proj_settings_500);
geodetic::projection::proj_file::LEO_GeodeticReferenceFile geofile_250 = geodetic::projection::proj_file::leoRefFileFromProjector(norad, proj_settings_250);
geodetic::projection::proj_file::writeReferenceFile(geofile_1000, directory + "/MODIS-1000.georef");
geodetic::projection::proj_file::writeReferenceFile(geofile_500, directory + "/MODIS-500.georef");
geodetic::projection::proj_file::writeReferenceFile(geofile_250, directory + "/MODIS-250.georef");
}
for (int i = 0; i < 2; i++)
{
cimg_library::CImg<unsigned short> image = reader.getImage250m(i);
@ -364,22 +412,8 @@ namespace eos
{
//pre_wb.resize(pre_wb.width() / 4, pre_wb.height() / 4);
// Setup Projecition
std::shared_ptr<geodetic::projection::LEOScanProjectorSettings_SCANLINE> proj_settings = std::make_shared<geodetic::projection::LEOScanProjectorSettings_SCANLINE>(
109.7, // Scan angle
-0.1, // Roll offset
0, // Pitch offset
-2.5, // Yaw offset
-2.0, // Time offset
pre_wb.width(), // Image width
true, // Invert scan
tle::getTLEfromNORAD(norad), // TLEs
reader.timestamps_250 // Timestamps
);
geodetic::projection::LEOScanProjector projector(proj_settings);
logger->info("Projected Channel 143 EQURAW...");
cimg_library::CImg<unsigned char> projected_image = geodetic::projection::projectLEOToEquirectangularMapped(pre_wb, projector, 2048 * 4, 1024 * 4, 3);
cimg_library::CImg<unsigned char> projected_image = geodetic::projection::projectLEOToEquirectangularMapped(pre_wb, projector_250, 2048 * 4, 1024 * 4, 3);
WRITE_IMAGE(projected_image, directory + "/MODIS-143-EQURAW-PROJ.png");
}
}
@ -396,27 +430,13 @@ namespace eos
EOS_ORBIT_HEIGHT,
EOS_MODIS_SWATH,
EOS_MODIS_RES1000);
//WRITE_IMAGE(corrected23, directory + "/MODIS-29-EQU-CORRECTED.png");
WRITE_IMAGE(corrected23, directory + "/MODIS-29-EQU-CORRECTED.png");
// Reproject to an equirectangular proj
{
// Setup Projecition
std::shared_ptr<geodetic::projection::LEOScanProjectorSettings_SCANLINE> proj_settings = std::make_shared<geodetic::projection::LEOScanProjectorSettings_SCANLINE>(
109.7, // Scan angle
-0.1, // Roll offset
0, // Pitch offset
-2.5, // Yaw offset
-2.0, // Time offset
image23.width(), // Image width
true, // Invert scan
tle::getTLEfromNORAD(norad), // TLEs
reader.timestamps_1000 // Timestamps
);
geodetic::projection::LEOScanProjector projector(proj_settings);
logger->info("Projected Channel 29...");
cimg_library::CImg<unsigned char> projected_image = geodetic::projection::projectLEOToEquirectangularMapped(image23, projector, 2048 * 4, 1024 * 4, 1);
//projected_image.crop(20778, 2853, 20778 + 5145, 2853 + 3573);
cimg_library::CImg<unsigned char> projected_image = geodetic::projection::projectLEOToEquirectangularMapped(image23, projector_1000, 2048 * 4, 1024 * 4, 1);
//projected_image.crop(18800, 2668, 18800 + 5516, 2668 + 3700);
WRITE_IMAGE(projected_image, directory + "/MODIS-29-EQU-PROJ.png");
}
}