気象衛星ひまわり物理量算出プログラム
プログラム
hband_pq.c
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<unistd.h>
#include<stdint.h>
uint8_t convert_to8bit(int bit, int value) {
return (uint8_t)((value * 255) / (pow(2, bit)-1));
}
uint16_t expand_to16bit(uint8_t value) {
return (uint16_t)(value * 257);
}
int read_band(FILE *fi) {
int band;
fseek(fi, 598+3, SEEK_SET);
fread(&band, sizeof(short int), 1, fi);
return band;
}
double read_cwl(FILE *fi) {
double cwl;
fseek(fi, 598+5, SEEK_SET);
fread(&cwl, sizeof(double), 1, fi);
return cwl;
}
int read_segsum(FILE *fi) {
int segsum;
fseek(fi, 1004+3, SEEK_SET);
segsum = fgetc(fi);
return segsum;
}
int read_colnum(FILE *fi) {
int colnum;
fseek(fi, 282+5, SEEK_SET);
fread(&colnum, sizeof(short int), 1, fi);
return colnum;
}
int read_rownum(FILE *fi) {
int rownum;
fseek(fi, 282+7, SEEK_SET);
fread(&rownum, sizeof(short int), 1, fi);
return rownum;
}
int read_validbit(FILE *fi) {
int validbit;
fseek(fi, 598+13, SEEK_SET);
fread(&validbit, sizeof(short int), 1, fi);
return validbit;
}
double read_radiance_coef(FILE *fi) {
double coef;
fseek(fi, 598+19, SEEK_SET);
fread(&coef, sizeof(double), 1, fi);
return coef;
}
double read_radiance_const(FILE *fi) {
double rconst;
fseek(fi, 598+27, SEEK_SET);
fread(&rconst, sizeof(double), 1, fi);
return rconst;
}
double read_albedo_coef(FILE *fi) {
double coef;
fseek(fi, 598+35, SEEK_SET);
fread(&coef, sizeof(double), 1, fi);
return coef;
}
double read_btemp_c0(FILE *fi) {
double c0;
fseek(fi, 598+35, SEEK_SET);
fread(&c0, sizeof(double), 1, fi);
return c0;
}
double read_btemp_c1(FILE *fi) {
double c1;
fseek(fi, 598+43, SEEK_SET);
fread(&c1, sizeof(double), 1, fi);
return c1;
}
double read_btemp_c2(FILE *fi) {
double c2;
fseek(fi, 598+51, SEEK_SET);
fread(&c2, sizeof(double), 1, fi);
return c2;
}
double read_lspeed(FILE *fi) {
double speed;
fseek(fi, 598+75, SEEK_SET);
fread(&speed, sizeof(double), 1, fi);
return speed;
}
double read_plank(FILE *fi) {
double plank;
fseek(fi, 598+83, SEEK_SET);
fread(&plank, sizeof(double), 1, fi);
return plank;
}
double read_boltzmann(FILE *fi) {
double bolt;
fseek(fi, 598+91, SEEK_SET);
fread(&bolt, sizeof(double), 1, fi);
return bolt;
}
int calc_blocklen(FILE *fi) {
int i, blocklen = 0;
short int tmp_blocklen;
for (i=0; i < 11; i++) {
fseek(fi, blocklen+1, SEEK_SET);
fread(&tmp_blocklen, sizeof(short int), 1, fi);
blocklen = blocklen + tmp_blocklen;
}
return blocklen;
}
double calc_linear_func(double coef, double cons, int x) {
return x * coef + cons;
}
double calc_linear_func_double(double coef, double cons, double x) {
return x * coef + cons;
}
double calc_btemp(double c0, double c1, double c2, double c, double h, double k, double cwl, double rad) {
double te, tel, ter, data;
data = rad * pow(10.0, 6.0);
tel = (h * c) / (k * cwl);
ter = 1.0 / (log(((2.0*h*c*c)/(pow(cwl,5.0)*data))+1.0));
te = tel * ter;
return c0 + c1 * te + c2 * pow(te, 2.0);
}
void write_data(double *phy_data, int x_len, int y_len, int dn, char *output_file) {
FILE *fo;
fo = fopen(output_file, "wb");
fwrite(&x_len, sizeof(int), 1, fo);
fwrite(&y_len, sizeof(int), 1, fo);
fwrite(phy_data, sizeof(double), x_len * y_len, fo);
fclose(fo);
return;
}
int main(int argc, char *argv[]) {
int i, fi_num=0, opt, bg, dn;
char *b = "2", *fo_name = "", *d = "1", *fi_name = "";
while ((opt = getopt(argc, argv, "r:b:o:d:")) != -1) {
switch (opt) {
case 'b':
b = optarg;
bg = atoi(b);
break;
case 'o':
fo_name = optarg;
break;
case 'd':
d = optarg;
dn = atoi(d);
break;
default:
printf("Usage: %s [-b number] [-d number] [-o string] arg1 ...\n", argv[0]);
return -1;
}
}
fi_name = argv[optind];
int validbit, colnum, rownum, blocklen, band;
double r_coef, r_cons;
FILE *fi;
fi = fopen(fi_name, "r");
validbit = read_validbit(fi);
band = read_band(fi);
colnum = read_colnum(fi);
rownum = read_rownum(fi);
r_coef = read_radiance_coef(fi);
r_cons = read_radiance_const(fi);
blocklen = calc_blocklen(fi);
fseek(fi, blocklen, SEEK_SET);
int *image_data = (int *)malloc((colnum/dn) * (rownum/dn) * sizeof(int));
unsigned short (*image)[colnum] = malloc(rownum * colnum * 2);
fread(image, 2, rownum*colnum, fi);
for (int y = 0; y < rownum; y++) {
for (int x = 0; x < colnum; x++) {
if (image[y][x] == 65535 || image[y][x] == 65534) {
image_data[(y/dn) * (colnum/dn) + (x/dn)] = 65535;
} else {
image_data[(y/dn) * (colnum/dn) + (x/dn)] = (int)image[y][x];
}
x += (dn-1);
}
y += (dn-1);
}
free(image);
double *radiances = (double *) malloc((colnum/dn) * (rownum/dn) * sizeof(double));
for (int y = 0; y < rownum/dn; y++) {
for (int x = 0; x < colnum/dn; x++) {
if (image_data[y * (colnum/dn) + x] == 65535) {
radiances[y * (colnum/dn) + x] = 65535;
} else {
radiances[y * (colnum/dn) + x] = calc_linear_func(r_coef, r_cons, image_data[y * (colnum/dn) + x]);
}
}
}
double *phy_data = (double *)malloc((colnum/dn) * (rownum/dn) * sizeof(double));
if (band < 7) {
double acoef;
acoef = read_albedo_coef(fi);
for (int y = 0; y < rownum/dn; y++) {
for (int x = 0; x < colnum/dn; x++) {
if (radiances[y * colnum/dn + x] == 65535) {
phy_data[y * (colnum/dn) + x] = (double)bg;
} else {
phy_data[y * (colnum/dn) + x] = calc_linear_func_double(acoef, 0, radiances[y * (colnum/dn) + x]);
}
}
}
} else {
double cwl, tc0, tc1, tc2, te, tb;
double k = 1.380 * pow(10.0,-23.0),h = 6.626*pow(10.0,-34.0), c = 2.998*pow(10.0,8.0);
cwl = read_cwl(fi);
cwl *= pow(10.0, -6.0);
tc0 = read_btemp_c0(fi);
tc1 = read_btemp_c1(fi);
tc2 = read_btemp_c2(fi);
for (int y = 0; y < rownum/dn; y++) {
for (int x = 0; x < colnum/dn; x++) {
if (radiances[y * colnum/dn + x] == 65535) {
phy_data[y * (colnum/dn) + x] = (double)bg;
} else {
phy_data[y * (colnum/dn) + x] = calc_btemp(tc0, tc1, tc2, c, h, k, cwl, radiances[y * (colnum/dn) + x]);
}
}
}
}
fclose(fi);
write_data(phy_data, colnum/dn, rownum/dn, dn, fo_name);
free(image_data);
free(radiances);
free(phy_data);
return 0;
}
コンパイル
$ cc -o hband_pq hband_pq.c -lm
使い方
$ ./hband_pq -d 1 -b 0 -o HS_H09_20251001_0000_B01_JP01_R10_S0101.DAT.raw HS_H09_20251001_0000_B01_JP01_R10_S0101.DAT
オプション
-d: 分解能 (例: 2 -> 1/2、3 -> 1/3)
-b: 背景色 (0〜255)
-o: 出力先ファイル