TC低菊渲染（污染）实验贴

wfIRstzhang · 发表于 2026-2-24 14:41

五彩斑斓的灰（不是）有色BD色阶测试

wfIRstzhang · 发表于 2026-2-24 14:41

本帖最后由 wfIRstzhang 于 2026-2-24 15:48 编辑

#include <algorithm>
#include <chrono>
#include <cmath>
#include <complex>
#include <cpgplot.h>
#include <iostream>
#include <limits>
#include <map>
#include <mfhdf.h>
#include <numbers>
#include <stdexcept>
#include <type_traits>
#include <vector>
#define PNG_OUTPUT
constexpr double Hc_k = 1.4387768775039337e4;
constexpr double H2c_2 = 1.1910429723971884140794892e8;
constexpr double Deg2rad = std::numbers::pi / 180;
constexpr double K2c = -273.15;
static int lowcol = 16, hicol = 55;
constexpr double brtemp(double bt, double lam_um) {
return Hc_k / lam_um / std::log(H2c_2 / pow(lam_um, 5) / bt + 1);
}
int segmented_linear_color(float val, float loclow, float lochigh,
int locmincol, int locmaxcol) {
return locmincol + (locmaxcol - locmincol) *
std::clamp<double>((val - loclow) / (lochigh - loclow), 0., 1.);
}
int color_brtemp(double brtempC) {
if (brtempC >= 9)
return segmented_linear_color(brtempC, 9, 33, 35, 16);
if (brtempC >= -31)
return segmented_linear_color(brtempC, -31, 9, 36, 45);
if (brtempC >= -42)
return 45;
if (brtempC >= -54)
return 46;
if (brtempC >= -64)
return 47;
if (brtempC >= -70)
return 48;
if (brtempC >= -76)
return 49;
if (brtempC >= -77.6)
return 50;
if (brtempC >= -79.3)
return 51;
if (brtempC >= -81)
return 52;
if (brtempC >= -85)
return 53;
if (brtempC >= -90)
return 54;
return 55;
}
void set_color() {
for (int i = 16; i < 36; i++)
{
float nsegp = float(i - 16) / 19;
float bg = nsegp * 0.6 + 0.4;
cpgscr(i, bg - 0.1,
bg, bg);
}
for (int i = 36; i < 45; i++)
{
float cof = pow(float(i - 35), 1.1) / pow(10, 1.1);
float use = pow(45 - i, 2) / pow(10, 2);
float gd = 0.5 * cof + 0.47;
// cpgscr(i, use + (1 - use) * (cof / 3 + 0.5),
// use + (1 - use) * (cof / 2.25 + 0.3),
// use + (1 - use) * (0.34 + cof / 7));
cpgscr(i, gd + 0.03, gd + 0.03, gd);
}
cpgscr(45, 0, 1, 0.5);
cpgscr(46, 0, 0.7, 0.35);
cpgscr(47, 0, 0.4, 0.2);
cpgscr(48, 1, 0.5, 0);
cpgscr(49, 1, 0, 0.5);
cpgscr(50, 0, 0.65, 1);
cpgscr(51, 0, 0.5, 1);
cpgscr(52, 0, 0.35, 1);
cpgscr(53, 1, 1, 0.2);
cpgscr(54, 1, 0.6, 0.2);
cpgscr(55, 1, 1, 1);
}
template<int32 Nt>
using SDtype =
std::conditional_t<Nt == 3, int8,
std::conditional_t<Nt == 4, char8,
std::conditional_t<Nt == 5, float32,
std::conditional_t<Nt == 6, float64,
std::conditional_t<Nt == 20, int8,
std::conditional_t<Nt == 21, uint8,
std::conditional_t<Nt == 22, int16,
std::conditional_t<Nt == 23, uint16,
std::conditional_t<Nt == 24, int32,
std::conditional_t<Nt == 25, uint32, void>>>>>>>>>>;
struct SDfile {
int32 sdid;
int32 ndataset, nattr;
template<class Tp>
struct Attribute {
int32 obj_id, index;
int32 type, count;
char name[256];
std::vector<Tp> data;
Attribute(int32 pobj_id, int32 pindex)
:obj_id(pobj_id), index(pindex) {
SDattrinfo(pobj_id, pindex, name, &type, &count);
data.assign(count, Tp());
SDreadattr(obj_id, index, data.data());
}
};
struct SDS {
int32 All_begin[MAX_VAR_DIMS]{};
template<class Tp>
struct Data_array {
int32 dimn;
int32 dims[MAX_VAR_DIMS];
Tp* data;
Tp operator[] (std::initializer_list<int32> dimensions) const {
auto it = dimensions.end();
size_t dimit = dimn - 1, pnow = 1, datanow = 0;
do
{
--it;
datanow += *it * pnow;
pnow *= dims[dimit--];
} while (it != dimensions.begin());
return data[datanow];
}
Data_array(const SDS* from_sds, const int32* vdims)
:dimn(from_sds->dimn) {
std::copy_n(vdims, from_sds->dimn, dims);
size_t tot = 1;
for (int i = 0; i < dimn; i++)
tot *= dims[i];
data = new Tp[tot];
}
~Data_array() { delete[] data; }
};
int32 sdsid;
char name[256];
int32 dimn;
int32 dims[MAX_VAR_DIMS];
int32 dattype;
int32 nattr;
template<class Tp>
Data_array<Tp> get() {
Data_array<Tp> res(this, dims);
SDreaddata(sdsid, All_begin, nullptr, dims, res.data);
return res;
}
template<class Tp>
Data_array<Tp> get(const std::vector<int32>& start, const std::vector<int32>& cntrd) {
auto mvas = start, mvac = cntrd;
Data_array<Tp> res(this, cntrd.data());
SDreaddata(sdsid, mvas.data(), nullptr, mvac.data(), res.data);
return res;
}
template<class Tp>
Attribute<Tp> attribute(int32 attr_index) {
return Attribute<Tp>(sdsid, attr_index);
}
template<class Tp>
Attribute<Tp> attribute(std::string_view attrname) {
return Attribute<Tp>(sdsid, SDfindattr(sdsid, attrname.data()));
}
SDS(int sds_id)
:sdsid(sds_id) {
SDgetinfo(sds_id, name, &dimn, dims, &dattype, &nattr);
}
~SDS() { SDendaccess(sdsid); }
};
SDS select(int32 index) {
return SDS(SDselect(sdid, index));
}
SDS select(std::string_view sdsname) {
return SDS(SDselect(sdid, SDnametoindex(sdid, sdsname.data())));
}
template<class Tp>
Attribute<Tp> attribute(int32 attr_index) {
return Attribute<Tp>(sdid, attr_index);
}
template<class Tp>
Attribute<Tp> attribute(std::string_view attrname) {
return Attribute<Tp>(sdid, SDfindattr(sdid, attrname.data()));
}
SDfile(std::string_view filename)
:sdid(SDstart(filename.data(), DFACC_READ)) {
SDfileinfo(sdid, &ndataset, &nattr);
}
~SDfile() { SDend(sdid); }
};
std::map<std::string, double> Center_Wl_in_um{
{"1", (0.62 + 0.67) / 2},
{"2", (0.841 + 0.876) / 2},
{"3", (0.459 + 0.479) / 2},
{"4", (0.545 + 0.565) / 2},
{"5", (1.23 + 1.25) / 2},
{"6", (1.628 + 1.652) / 2},
{"7", (2.105 + 2.155) / 2},
{"8", (0.405 + 0.42) / 2},
{"9", (0.438 + 0.448) / 2},
{"10", (0.483 + 0.493) / 2},
{"11", (0.526 + 0.536) / 2},
{"12", (0.546 + 0.556) / 2},
{"13hi", (0.662 + 0.672) / 2},
{"13lo", (0.662 + 0.672) / 2},
{"14hi", (0.673 + 0.683) / 2},
{"14lo", (0.673 + 0.683) / 2},
{"15", (0.743 + 0.753) / 2},
{"16", (0.862 + 0.877) / 2},
{"17", (0.89 + 0.92) / 2},
{"18", (0.931 + 0.941) / 2},
{"19", (0.915 + 0.965) / 2},
{"20", (3.66 + 3.84) / 2},
{"21", (3.929 + 3.989) / 2},
{"22", (3.929 + 3.989) / 2},
{"23", (4.02 + 4.08) / 2},
{"24", (4.433 + 4.498) / 2},
{"25", (4.482 + 4.549) / 2},
{"26", (1.36 + 1.39) / 2},
{"27", (6.535 + 6.895) / 2},
{"28", (7.175 + 7.475) / 2},
{"29", (8.4 + 8.7) / 2},
{"30", (9.58 + 9.88) / 2},
{"31", (10.78 + 11.28) / 2},
{"32", (11.77 + 12.27) / 2},
{"33", (13.185 + 13.485) / 2},
{"34", (13.485 + 13.785) / 2},
{"35", (13.785 + 14.085) / 2},
{"36", (14.085 + 14.385) / 2}
};
template<class Tp>
Tp* interp(const SDfile::SDS::Data_array<Tp>& init, int dsth, int dstw);
int main(int argc, char** argv) {
using namespace std;
if (argc < 2)
throw invalid_argument("必须在 argv[1] 提供文件名");
auto lasttime = chrono::high_resolution_clock::now();
SDfile file(argv[1]);
/*
Latitude,Longitude,EV_1KM_RefSB,EV_1KM_RefSB_Uncert_Indexes,EV_1KM_Emissive,
EV_1KM_Emissive_Uncert_Indexes,EV_250_Aggr1km_RefSB,EV_250_Aggr1km_RefSB_Uncert_Indexes,
EV_250_Aggr1km_RefSB_Samples_Used,EV_500_Aggr1km_RefSB,EV_500_Aggr1km_RefSB_Uncert_Indexes,
EV_500_Aggr1km_RefSB_Samples_Used,Height,SensorZenith,SensorAzimuth,
Range,SolarZenith,SolarAzimuth,gflags,EV_Band26,EV_Band26_Uncert_Indexes,Band_250M,
Band_500M,Band_1KM_RefSB,Band_1KM_Emissive,Noise in Thermal Detectors,
Change in relative responses of thermal detectors,DC Restore Change for Thermal Bands,
DC Restore Change for Reflective 250m Bands,DC Restore Change for Reflective 500m Bands,
DC Restore Change for Reflective 1km Bands
*/
string main_sdsname = "EV_1KM_Emissive";
auto sds = file.select(main_sdsname);
auto fulldat = sds.get<uint16>();
auto getbandname = [&] (int bdid) {
string bandname;
auto bds = sds.attribute<char>("band_names");
int cnt = 0;
for (auto ch : bds.data)
if (ch == ',')
{
if (cnt == bdid)
break;
bandname.clear();
++cnt;
}
else
bandname.push_back(ch);
return bandname;
};
auto readband = [&] (int bdid) {
auto emmmat = new float[sds.dims[1] * sds.dims[2]];
string bandname = getbandname(bdid);
float ro = sds.attribute<float>("radiance_offsets").data[bdid],
rs = sds.attribute<float>("radiance_scales").data[bdid];
for (int i = 0; i < sds.dims[1]; i++)
for (int j = 0; j < sds.dims[2]; j++)
emmmat[i * sds.dims[2] + j] = (fulldat[{bdid, i, j}] - ro) * rs;
for (int i = 0; i < sds.dims[1] * sds.dims[2]; i++)
emmmat[i] = brtemp(emmmat[i], Center_Wl_in_um.at(bandname));
return emmmat;
};
int bdid = 10;
string lbinfo = main_sdsname + " Band\\d" + getbandname(bdid) + "\\u";
//lbinfo += " Anomaly";
lbinfo += " " + to_string(Center_Wl_in_um.at(getbandname(bdid))) + "um";
auto emmmat = readband(bdid);
string outname = string("'") + argv[1] + "-";
for (auto ch : lbinfo)
if (ch == ' ')
outname.push_back('-');
else if (ch != '\\')
outname.push_back(ch);
outname += ".png'";
auto solarzsds = file.select("SolarZenith");
auto solarzl0 = solarzsds.get<int16>();
SDfile::SDS::Data_array<float> solarzl1(&solarzsds, solarzsds.dims);
for (int i = 0; i < solarzl0.dims[0] * solarzl0.dims[1]; i++)
solarzl1.data[i] = solarzl0.data[i] / 1e2;
auto solarzl2 = interp<float>(solarzl1, sds.dims[1], sds.dims[2]);
auto latsds = file.select("Latitude");
auto latl1 = latsds.get<float>();
auto latl2 = interp<float>(latl1, sds.dims[1], sds.dims[2]);
auto lonsds = file.select("Longitude");
auto lonl1 = lonsds.get<float>();
for (int i = 0; i < lonl1.dims[0] * lonl1.dims[1]; i++)
if (lonl1.data[i] < 0)
lonl1.data[i] += 360;
auto lonl2 = interp<float>(lonl1, sds.dims[1], sds.dims[2]);
// for (int i = 0; i < sds.dims[1] * sds.dims[2]; i++)
// emmmat[i] = brtemp(emmmat[i] / cos(Deg2rad * solarzl2[i]),
// Center_Wl_in_um.at(getbandname(bdid)));
auto tmp = file.attribute<char>("CoreMetadata.0").data;
string coremeta(tmp.begin(), tmp.end());
auto getcmt = [&] (string obj) {
auto posm = coremeta.find("OBJECT = " + obj);
auto pos = coremeta.find("VALUE = ", posm);
int offset = 23, backoff = 0;
if (coremeta[pos + offset] == '"')
{
++offset;
++backoff;
}
auto end = coremeta.find('\n', pos + offset);
return coremeta.substr(pos + offset, end - pos - offset - backoff);
};
string device = getcmt("ASSOCIATEDSENSORSHORTNAME") + "/" +
getcmt("ASSOCIATEDPLATFORMSHORTNAME");
string endtime = getcmt("RANGEENDINGDATE") + " " +
getcmt("RANGEENDINGTIME");
while (endtime.back() == '0')
endtime.pop_back();
endtime.pop_back();
endtime += "\\uUTC\\d";
cout << "读取数据用时： " <<
chrono::duration_cast<chrono::microseconds>(chrono::high_resolution_clock::now()
- lasttime).count() / 1e6 << " s\n";
lasttime = chrono::high_resolution_clock::now();
#ifdef PNG_OUTPUT
if (cpgopen("tmp.ps/VCPS") < 0)
#else
if (cpgopen("/XWINDOW") < 0)
#endif
throw runtime_error("打开PGPLOT设备错误");
cpgbbuf();
auto latmmx =
minmax_element(latl1.data, latl1.data + latl1.dims[0] * latl1.dims[1]),
lonmmx =
minmax_element(lonl1.data, lonl1.data + lonl1.dims[0] * lonl1.dims[1]);
cpgenv(*lonmmx.first, *lonmmx.second, *latmmx.first, *latmmx.second
, 1, 1);
//cpgenv(117, 135, 10, 25, 1, 1);
cpglab("Lon", "Lat", argv[1]);
cpgmtxt("R", 1, 0, 0, lbinfo.c_str());
cpgmtxt("B", 2.5, 0.15, 0, device.c_str());
cpgmtxt("B", 2.5, 0.6, 0, endtime.c_str());
set_color();
cpgscir(lowcol, hicol);
for (int i = 0; i < sds.dims[1] * sds.dims[2]; i++)
{
if (isnan(emmmat[i]))
continue;
int col = color_brtemp(emmmat[i] + K2c);
int rind = i / sds.dims[2], cind = i % sds.dims[2];
double colrng = (float)abs(cind - sds.dims[2] / 2) / sds.dims[2];
double siz = 8e-3 * exp(3.5 * colrng + 0.8 * sqrt(abs(latl2[i]) * 1.3e-4)) +
4e-4 * exp(4 * pow(colrng, 2));
cpgpixl(&col, 1, 1, 1, 1, 1, 1, lonl2[i],
lonl2[i] + siz, latl2[i], latl2[i] + siz);
}
cpgsch(0.8);
cpgsci(color_brtemp(33));
cpgmtxt("R", 2.5, 0, 0, "W");
cpgsci(color_brtemp(24));
cpgmtxt("R", 2.5, 0, -1, "M");
cpgsci(color_brtemp(9.5));
cpgmtxt("R", 2.5, 0, -2, "G");
cpgsci(color_brtemp(8.5));
cpgmtxt("R", 2.5, 0.1, 0, "O");
cpgsci(color_brtemp(-30.5));
cpgmtxt("R", 2.5, 0.1, -1, "W");
cpgsci(color_brtemp(-36.5));
cpgmtxt("R", 2.5, 0.2, 0, "DG");
cpgsci(color_brtemp(-48));
cpgmtxt("R", 2.5, 0.3, 0, "MG");
cpgsci(color_brtemp(-59));
cpgmtxt("R", 2.5, 0.4, 0, "LG");
cpgsci(color_brtemp(-67));
cpgmtxt("R", 2.5, 0.5, 0, "B");
cpgsci(color_brtemp(-73));
cpgmtxt("R", 2.5, 0.6, 0, "W");
cpgsci(color_brtemp(-76.5));
cpgmtxt("R", 2.5, 0.7, 0, "C");
cpgsci(color_brtemp(-78));
cpgmtxt("R", 2.5, 0.7, -1, "M");
cpgsci(color_brtemp(-80.5));
cpgmtxt("R", 2.5, 0.7, -2, "G");
cpgsci(color_brtemp(-81.5));
cpgmtxt("R", 2.5, 0.8, 0, "C");
cpgsci(color_brtemp(-87));
cpgmtxt("R", 2.5, 0.8, -1, "D");
cpgsci(color_brtemp(-100));
cpgmtxt("R", 2.5, 0.8, -2, "G");
cpgebuf();
cpgend();
cout << "PGPLOT输出用时： " <<
chrono::duration_cast<chrono::microseconds>(chrono::high_resolution_clock::now()
- lasttime).count() / 1e6 << " s\n";
lasttime = chrono::high_resolution_clock::now();
#ifdef PNG_OUTPUT
system(("magick -density 240 tmp.ps " + outname).c_str());
cout << "PNG转换用时： " <<
chrono::duration_cast<chrono::microseconds>(chrono::high_resolution_clock::now()
- lasttime).count() / 1e6 << " s\n";
cout << "输出到PNG文件： " << outname << endl;
system("rm tmp.ps");
#endif
delete[] solarzl2;
delete[] latl2;
delete[] lonl2;
delete[] emmmat;
return 0;
}
template<class Tp>
Tp* interp(const SDfile::SDS::Data_array<Tp>& init, int dsth, int dstw) {
auto res = new Tp[dsth * dstw];
res[0] = init[{0, 0}];
res[dsth * dstw - 1] = init[{init.dims[0] - 1, init.dims[1] - 1}];
double deltaX = double(init.dims[0] - 1) / (dsth - 1), deltaY = double(init.dims[1] - 1) / (dstw - 1);
for (int i = 0; i < dsth; i++)
for (int j = 0; j < dstw; j++)
if ((i || j) && (i != dsth - 1 || j != dstw - 1))
{
double mappedX = deltaX * i, mappedY = deltaY * j;
double cX = mappedX - floor(mappedX), cY = mappedY - floor(mappedY);
int x0 = mappedX, x1 = std::min<int>(ceil(mappedX), init.dims[0] - 1);
int y0 = mappedY, y1 = std::min<int>(ceil(mappedY), init.dims[1] - 1);
res[i * dstw + j] = (1 - cX) * (1 - cY) * init[{x0, y0}] +
(1 - cX) * cY * init[{x0, y1}] +
cX * (1 - cY) * init[{x1, y0}] +
cX * cY * init[{x1, y1}];
}
return res;
}

复制代码

wfIRstzhang · 发表于 2026-2-24 16:03

wfIRstzhang · 发表于 2026-2-26 07:27

Band01 VIS和真没有加海温反演修正的粗糙伪VIS

wfIRstzhang · 发表于 2026-2-26 07:30

		自动登录	找回密码
密码			立即注册

[耍宝娱乐] TC低菊渲染（污染）实验贴

本帖子中包含更多资源

本帖子中包含更多资源

本帖子中包含更多资源

本帖子中包含更多资源

本帖子中包含更多资源