#include <math.h>
#include <stdio.h>
#include <string.h>
#include "read_gome.h"
#include "time_class.h"
#include "peteys_tmpl_lib.h"
//generates a single spectra in a geometric series:
float * default_spectra(float lambda0, float lambdan, int n) {
float *spec;
spec=new float[n+1];
for (long i=0; i<=n; i++) {
spec[i]=lambda0*pow(lambdan/lambda0, (float) i/(float) n);
}
return spec;
}
//generates a single spectra in a geometric series:
float * default_spectra_mid(float lambda0, float lambdan, int n) {
float *spec;
spec=new float[n];
for (long i=0; i<n; i++) {
spec[i]=lambda0*pow(lambdan/lambda0, (i+0.5)/(float) n);
}
return spec;
}
#define NHEAD 13
#define MAXLL 200
#define NPPERSCAN 3
//#define NCHAN NPMD+5
#define NBANDMAX 4
#define NRECBAND 5
#define NSPECRAWMAX NBAND*1024
time_class convert_date(char *date) {
time_class t;
int year;
int mon;
int day;
int hour;
int min;
float sec;
sscanf(date, "%d2", &day);
sscanf(date+7, "%d4", &year);
sscanf(date+12, "%d2", &hour);
sscanf(date+15, "%d2", &min);
sscanf(date+18, "%f", &sec);
date[6]=0;
if (strcmp(date+3, "JAN")==0) {
mon=1;
} else if (strcmp(date+3, "FEB")==0) {
mon=2;
} else if (strcmp(date+3, "MAR")==0) {
mon=3;
} else if (strcmp(date+3, "APR")==0) {
mon=4;
} else if (strcmp(date+3, "MAY")==0) {
mon=5;
} else if (strcmp(date+3, "JUN")==0) {
mon=6;
} else if (strcmp(date+3, "JUL")==0) {
mon=7;
} else if (strcmp(date+3, "AUG")==0) {
mon=8;
} else if (strcmp(date+3, "SEP")==0) {
mon=9;
} else if (strcmp(date+3, "OCT")==0) {
mon=10;
} else if (strcmp(date+3, "NOV")==0) {
mon=11;
} else if (strcmp(date+3, "DEC")==0) {
mon=12;
} else {
fprintf(stderr, "Unrecognized month: %s\n", date+3);
}
t.init(year, mon, day, hour, min, sec);
return t;
}
#define SUN 0
#define GROUND 1
#define DUMMY 1
int NBAND=4;
float glob_specmin[NBANDMAX]={285, 313, 405, 602};
float glob_specmax[NBANDMAX]={313, 405, 602, 790};
int glob_bind[NBANDMAX]={0, 1, 2, 3};
void rg_select_band(int b1b, int b2b, int b3, int b4) {
int bflag[NBANDMAX];
int k;
// float specmin[NBANDMAX]={285, 313, 405, 602};
// float specmax[NBANDMAX]={313, 405, 602, 790};
bflag[0]=b1b; bflag[1]=b2b; bflag[2]=b3; bflag[3]=b4;
k=0;
for (int i=0; i<NBANDMAX; i++) {
if (bflag[i]) {
//glob_specmin[k]=specmin[i];
//glob_specmax[k]=specmax[i];
glob_bind[k]=i;
k++;
}
}
NBAND=k;
printf("%d bands selected\n", NBAND);
}
float ** rg_get_band(FILE *in, int &n, char tp, char dum=0) {
char line[MAXLL];
char band[2];
float integration_time;
float spec0;
float specf;
float **data;
fgets(line, MAXLL, in);
//printf(line);
if (tp == SUN) {
sscanf(line+9, "%f %f %d", &spec0, &specf, &n);
//printf(line+9);
} else {
if (strncmp(line, "Band", 4) != 0) {
fseek(in, -strlen(line), SEEK_CUR);
n=0;
return NULL;
}
sscanf(line+7, "%f %f %f %d", &integration_time, &spec0, &specf, &n);
//printf("%f %f %f\n", integration_time, spec0, specf);
}
//printf("%d frequencies found in %6s\n", n, line);
if (dum != 0) {
for (long i=0; i<n; i++) fgets(line, MAXLL, in);
return NULL;
}
data=new float *[NRECBAND];
data[0]=new float[n*NRECBAND];
for (int i=1; i<NRECBAND; i++) data[i]=data[0]+i*n;
for (long i=0; i<n; i++) {
fgets(line, MAXLL, in);
sscanf(line, "%f %f %f %f %f", data[0]+i, data[1]+i, data[2]+i, data[3]+i, data[4]+i);
//printf("%f\n", data[0]+i);
//printf("%s", line);
}
return data;
}
int rg_get_raw(FILE *in, float *specraw, float *countsraw, char tp) {
//we use four bands and only take them within these limits:
//(we try to form a continuous spectra}
int sind, find; //bounds on the current band
int nraw_ttl;
float **specdata;
int nraw;
int nband2;
int *bind;
if (tp == SUN) {
nband2=NBANDMAX;
bind=new int[NBANDMAX];
for (int i=0; i<NBANDMAX; i++) bind[i]=i;
} else {
nband2=NBAND;
bind=glob_bind;
}
nraw_ttl=0;
for (long j=0; j<nband2; j++) {
specdata=rg_get_band(in, nraw, tp);
if (specdata==NULL) return 0; //empty pixel
sind=bin_search(specdata[0], nraw, glob_specmin[bind[j]], -1);
find=bin_search(specdata[0], nraw, glob_specmax[bind[j]], -1);
for (int h=sind; h<=find; h++) {
specraw[nraw_ttl]=specdata[0][h];
//printf("%f\n", specraw[nraw_ttl]);
countsraw[nraw_ttl]=specdata[1][h];
nraw_ttl++;
}
delete [] specdata[0];
delete [] specdata;
}
if (tp == SUN) delete [] bind;
return nraw_ttl;
}
void rg_interpolate_raw(float *specraw, float *countsraw, int nraw, float *spec, float *counts, int nchan) {
int l1, l2; //for interpolating the spectra
long lastind;
lastind=-1;
//interpolate the raw sunshine spectra to our grid:
l1=bin_search(specraw, nraw, spec[0], lastind);
for (long j=0; j<nchan; j++) {
l2=bin_search(specraw, nraw, spec[j+1], lastind);
//printf("%f %f %d %d\n", spec[j], spec[j+1], l1, l2);
counts[j]=0;
for (long h=l1; h<l2; h++) counts[j]+=countsraw[h];
counts[j]/=(l2-l1);
l1=l2;
}
}
gome_data * read_gome(char *fname,
float **spec,
int * nchan,
int nb,
int get_pmd) {
int npix;
int ipi;
int isc;
int k;
FILE *in;
char line[MAXLL];
double vald[NPPERSCAN];
float val[10];
char datestr[32];
gome_data *data;
double missing=9.99999e-99;
float lon;
float *sunshine; //sunshine spectra
float countsraw[NSPECRAWMAX];
float specraw[NSPECRAWMAX];
int nraw;
int ncum;
int nspec; //number of spectra, not including PMD
data=new gome_data;
in=fopen(fname, "r");
for (long i=0; i<NHEAD; i++) {
fgets(line, MAXLL, in);
//printf(line);
}
nspec=0;
for (int i=0; i<nb; i++) nspec+=nchan[i];
data->npix=0;
if (get_pmd==1) data->nchan=nspec+NAUX+NPMD; else data->nchan=nspec+NAUX;
data->missing=0./0.;
//read in the sunshine spectra:
//printf("Reading sun spectra\n");
sunshine=new float[nspec];
ncum=0;
if (nb > 0) {
nraw=rg_get_raw(in, specraw, countsraw, SUN);
for (int j=0; j<nb; j++) {
rg_interpolate_raw(specraw, countsraw, nraw, spec[j], sunshine+ncum, nchan[j]);
ncum+=nchan[j];
}
} else {
//do a dummy read:
for (int j=0; j<NBANDMAX; j++) {
rg_get_band(in, nraw, SUN, DUMMY);
}
}
fgets(line, MAXLL, in);
sscanf(line+46, "%d", &npix);
printf("%d ground pixels contained in %s\n", npix, fname);
data->data=new gome_rec[npix];
for (long i=0; i<npix; i++) {
data->data[data->npix].counts=new float[data->nchan];
fgets(line, MAXLL, in);
sscanf(line+12, "%d %d %d", &ipi, &k, &isc);
//printf(line);
//printf("Reading ground pixel: %d %d %d\n", i, ipi, isc);
fgets(line, MAXLL, in);
//printf("%s\n", line);
data->data[data->npix].date=convert_date(line);
data->data[data->npix].date.write_string(datestr);
fgets(line, MAXLL, in);
fgets(line, MAXLL, in);
fgets(line, MAXLL, in);
sscanf(line, "%f %f %f %f", val, val+1,
data->data[data->npix].counts,
data->data[data->npix].counts+1);
fgets(line, MAXLL, in);
sscanf(line, "%f %f %f %f", val, val+1,
data->data[data->npix].counts+2,
data->data[data->npix].counts+3);
fgets(line, MAXLL, in);
fgets(line, MAXLL, in);
//printf(line);
sscanf(line, "%f %f %f %f %f %f %f %f %f %f",
val, val+1, val+2, val+3, val+4, val+5, val+6, val+7,
&data->data[data->npix].lat, &lon);
if (lon > 180.) lon=lon-360.;
data->data[data->npix].lon=lon;
//store the pixel number:
data->data[data->npix].counts[4]=ipi;
data->data[data->npix].counts[5]=isc;
fgets(line, MAXLL, in);
//printf(line);
//printf("(%f, %f) %s; %ld\n", data.data[data.nscan].lon[isc],
// data.data[data.nscan].lat[isc], datestr, isc);
//read in and throw away (or not) the PMD data:
if (get_pmd==1) {
for (long j=0; j<16; j++) {
fgets(line, MAXLL, in);
fscanf(in, "%lg%lg%lg", vald, vald+1, vald+2);
if (vald[0] == missing) data->data[data->npix].counts[NAUX+j*3]=data->missing;
else data->data[data->npix].counts[NAUX+j*3]=vald[0];
if (vald[1] == missing) data->data[data->npix].counts[NAUX+j*3+1]=data->missing;
else data->data[data->npix].counts[NAUX+j*3+1]=vald[1];
if (vald[2] == missing) data->data[data->npix].counts[NAUX+j*3+2]=data->missing;
else data->data[data->npix].counts[NAUX+j*3+2]=vald[2];
//printf(line);
}
} else {
for (long j=0; j<16; j++) fgets(line, MAXLL, in);
}
if (nb > 0) {
int offset;
if (get_pmd==1) offset=NAUX+NPMD; else offset=NAUX;
nraw=rg_get_raw(in, specraw, countsraw, GROUND);
if (nraw==0) continue; //empty pixel
ncum=offset;
for (int j=0; j<nb; j++) {
rg_interpolate_raw(specraw, countsraw, nraw, spec[j], data->data[data->npix].counts+ncum, nchan[j]);
ncum+=nchan[j];
}
//divide by the sunshine spectra:
for (int j=0; j<nspec; j++) data->data[data->npix].counts[j+offset]/=sunshine[j];
} else {
//do a dummy read:
for (int j=0; j<NBAND; j++) {
rg_get_band(in, nraw, GROUND, DUMMY);
if (nraw==0) break;
}
if (nraw==0) continue; //empty pixel
}
data->npix++;
}
printf("%d pixels extracted from %s\n", data->npix, fname);
fclose(in);
delete [] sunshine;
return data;
}
void delete_gome_data(gome_data *data) {
for (int4_t i=0; i<data->npix; i++) {
delete [] data->data[i].counts;
}
delete [] data->data;
delete data;
}
#define COMMENT_LENGTH 500
//***note: the 'swap_end' parameter exists for compatibility reasons only
// and currently doesn't do anything
int4_t rawread_gome_head(char *filename, time_class & t1, time_class & t2,
int4_t &nchan, int swap_end) {
FILE *fs;
char comment[COMMENT_LENGTH];
int4_t nscan;
int4_t year;
double day;
//open the file:
fs=fopen(filename, "r");
if (fs == NULL) {
fprintf(stderr, "File, %s, could not be opened for reading\n", filename);
nchan=0;
return -1;
}
//read the comment header:
fread(comment, sizeof(char), COMMENT_LENGTH, fs);
fprintf(stderr, "%s", comment);
//read the attributes:
fread(&year, sizeof(year), 1, fs);
fread(&day, sizeof(day), 1, fs);
//this is awkward:
t1.init(year, 1, 1, 0, 0, 0);
t1.add(day);
fread(&year, sizeof(year), 1, fs);
fread(&day, sizeof(day), 1, fs);
t2.init(year, 1, 1, 0, 0, 0);
t2.add(day);
fread(&nscan, sizeof(nscan), 1, fs);
fread(&nchan, sizeof(nchan), 1, fs);
//fread(&missing, sizeof(missing), 1, fs);
fclose(fs);
return nscan;
}
int4_t rawwrite_gome(char * filename, gome_data *data, int swap_end) {
FILE *fs;
int4_t nwrit;
char comment[COMMENT_LENGTH];
char tstring1[30], tstring2[30];
int4_t year;
double day;
//open the file:
fs=fopen(filename, "w");
if (fs == NULL) {
fprintf(stderr, "File, %s, could not be opened for writing\n", filename);
return 1;
}
//write the headers:
//comment header:
data->data[0].date.write_string(tstring1);
data->data[data->npix-1].date.write_string(tstring2);
sprintf(comment, "GOME or GOME-derived data in raw binary format\nDate range:%30s-%30s\n%10d scan lines; %3d channels\n",
tstring1, tstring2, data->npix, data->nchan);
fwrite(comment, sizeof(char), COMMENT_LENGTH, fs);
//attributes:
year=(int4_t) data->data[0].date.year();
day=data->data[0].date.doy();
fwrite(&year, sizeof(year), 1, fs);
fwrite(&day, sizeof(day), 1, fs);
year=(int4_t) data->data[data->npix-1].date.year();
day=data->data[data->npix-1].date.doy();
fwrite(&year, sizeof(year), 1, fs);
fwrite(&day, sizeof(day), 1, fs);
fwrite(&data->npix, sizeof(data->npix), 1, fs);
fwrite(&data->nchan, sizeof(data->nchan), 1, fs);
fwrite(&data->missing, sizeof(data->missing), 1, fs);
//write out the data scan-line by scan-line:
for (int4_t i=0; i<data->npix; i++) {
year=(int4_t) data->data[i].date.year();
day=data->data[i].date.doy();
//try to align everything with the word boundaries
//for easy byte-swapping:
fwrite(&year, sizeof(year), 1, fs);
fwrite(&day, sizeof(day), 1, fs);
fwrite(&data->data[i].lon, sizeof(float), 1, fs);
fwrite(&data->data[i].lat, sizeof(float), 1, fs);
fwrite(data->data[i].counts, sizeof(float), data->nchan, fs);
}
fclose(fs);
return 0;
}
gome_data * rawread_gome(char * filename, int swap_end) {
FILE *fs;
time_class t1, t2;
char comment[COMMENT_LENGTH];
gome_data *data;
int4_t year;
double day;
//open the file:
fs=fopen(filename, "r");
if (fs == NULL) {
fprintf(stderr, "File, %s, could not be opened for reading\n", filename);
return NULL;
}
//read the comment header:
fread(comment, sizeof(char), COMMENT_LENGTH, fs);
fprintf(stderr, "%s", comment);
data=new gome_data;
//read the attributes:
fread(&year, sizeof(year), 1, fs);
fread(&day, sizeof(day), 1, fs);
fread(&year, sizeof(year), 1, fs);
fread(&day, sizeof(day), 1, fs);
fread(&data->npix, sizeof(data->npix), 1, fs);
fread(&data->nchan, sizeof(data->nchan), 1, fs);
fread(&data->missing, sizeof(data->missing), 1, fs);
data->data=new gome_rec[data->npix];
//read the data pixel by pixel:
for (long i=0; i<data->npix; i++) {
//allocate space:
fread(&year, sizeof(year), 1, fs);
fread(&day, sizeof(day), 1, fs);
data->data[i].date.init(year, 1, 1, 0, 0, 0);
data->data[i].date.add(day);
fread(&data->data[i].lon, sizeof(float), 1, fs);
fread(&data->data[i].lat, sizeof(float), 1, fs);
data->data[i].counts=new float[data->nchan];
fread(data->data[i].counts, sizeof(float), data->nchan, fs);
}
fclose(fs);
return data;
}