|
|
楼主 |
发表于 2026-2-22 21:27
|
显示全部楼层
本帖最后由 wfIRstzhang 于 2026-2-22 21:38 编辑
这几张图观感不好，不过主要是看看不同波段对低云的体现
- #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 float low = -60, high = 40;
- static int lowcol = 16, hicol = 100;
- constexpr double brtemp(double bt, double lam_um) {
- return Hc_k / lam_um / std::log(H2c_2 / pow(lam_um, 5) / bt + 1);
- }
- int color_brtemp(double brtempC) {
- double zonnap = std::clamp((brtempC - low) / (high - low), 0., 1.);
- int mincol, maxcol;
- cpgqcir(&mincol, &maxcol);
- return mincol + (maxcol - mincol) * zonnap;
- }
- void set_color() {
- // cpgscr(16, 0.5, 0.5, 0.75);
- // cpgscr(15, 0.4, 0.4, 0.6);
- // cpgscr(14, 0.3, 0.3, 0.45);
- // cpgscr(13, 0.2, 0.2, 0.3);
- // cpgscr(12, 0.1, 0.1, 0.15);
- int midlen = (hicol - lowcol) / 2;
- for (int i = lowcol; i <= lowcol + midlen; i++)
- {
- float base = float(i - lowcol) / midlen;
- cpgscr(i, base, base, base);
- }
- for (int i = lowcol + midlen + 1; i <= lowcol + midlen * 2; i++)
- {
- float base = float(lowcol + midlen * 2 - i) / midlen;
- cpgscr(i, base * 0.6, base * 0.6, base);
- }
- }
- 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 = 14;
- 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)));
- float avg = 0;
- int nnancnt = 0;
- for (int i = 0; i < sds.dims[1] * sds.dims[2]; i++)
- if (!isnan(emmmat[i]))
- {
- avg += emmmat[i];
- ++nnancnt;
- }
- avg /= nnancnt;
- 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] - avg);
- 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);
- }
- cpgwedg("RI", 2, 3, low, high, "^C");
- 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;
- }
|
本帖子中包含更多资源
您需要 登录 才可以下载或查看，没有账号？立即注册
×
|
发表于 2026-2-20 21:05
