00001 #include <string.h>
00002 #include "jsoc_main.h"
00003 #include "drms.h"
00004 #include "fitsexport.h"
00005 #include "drms_names.h"
00006
00007 #define NOT_SPECIFIED "***Not Specified***"
00008 #define DIE(msg) {fprintf(stderr,"$$$$ %s: %s\n", module_name, msg); return 1;}
00009
00010 int image_magrotate(void *, int, int, int, float, float, float, float,
00011 void **, int *, int *, int, int);
00012
00013 ModuleArgs_t module_args[] =
00014 {
00015 {ARG_STRING, "dsinp", NOT_SPECIFIED, "Input series query"},
00016 {ARG_STRING, "dsout", NOT_SPECIFIED, "Output series"},
00017 {ARG_STRING, "pathout", "/scr21/jsoc/aia/synoptic", "out base path"},
00018 {ARG_STRING, "drms", "0", "write to drms series"},
00019 {ARG_STRING, "file", "1", "write to hdir FITS"},
00020 {ARG_STRING, "rescale", "1", "rescale to fixed plate scale"},
00021 {ARG_STRING, "regrid", "1", "regrid type 0: nearest neighbor, 1: bicubic"},
00022 {ARG_STRING, "scale_to", "0.6", "rescale to fixed plate scale"},
00023 {ARG_STRING, "do_stretchmarks", "0", "fill in empty pixels created"},
00024 {ARG_FLAG, "h", "0", "Print usage message and quit"},
00025 {ARG_FLAG, "v", "0", "verbose flag"},
00026 {ARG_END}
00027 };
00028
00029 char *module_name = "aia_synoptic", *imtype;
00030 int rescale = 1, regridtype = 1, verbose = 0, do_stretchmarks = 0;
00031 float scale_to = 0.6, scale_by;
00032
00033 int nice_intro(int help)
00034 {
00035 int usage = cmdparams_get_int(&cmdparams, "h", NULL) != 0;
00036 verbose = cmdparams_get_int(&cmdparams, "v", NULL) != 0;
00037 if (usage || help) {
00038 printf("aia_lev1p5 {-h} {-v} dsinp=series_record_spec dsout=output_series\n"
00039 " -h: print this message\n"
00040 " -v: verbose\n"
00041 "pathout=<output directory path>, default=/scr21/jsoc/aia/synoptic\n"
00042 "drms=1, write output to DRMS series, default=0\n"
00043 "file=1, write output to hdir FITS file, default=1\n"
00044 "rescale=0, do not rescale to fixed plate scale, default=1\n"
00045 "regrid=0, nearest neighbor, default=1, bicubic\n"
00046 "scale_to=0.6, rescale to fixed plate scale, default=0.6\n"
00047 "do_stretchmarks=0, fill in empty pixels created, default=0\n"
00048 "dsinp=<recordset query> as <series>{[record specifier]} - required\n"
00049 "dsout=<series> - required if output to drms series\n");
00050 return(1);
00051 }
00052 return(0);
00053 }
00054
00055 void sprint_time_ISO (char *tstring, TIME t)
00056 {
00057 sprint_at(tstring,t);
00058 tstring[4] = tstring[7] = '-';
00059 tstring[10] = 'T';
00060 tstring[19] = '\0';
00061 }
00062
00063 int check_outpath(int yr, int mo, int da, int hr, char *outpathbase)
00064 {
00065 struct stat sbuf;
00066 char outfilename[512];
00067
00068 if (stat(outpathbase, &sbuf)) {
00069 fprintf(stderr, "Cant stat output base path\n"); return 1;
00070 }
00071 sprintf(outfilename, "%s/%4.4d", outpathbase, yr);
00072 if (stat(outfilename, &sbuf)) {
00073 if (mkdir(outfilename, 0777)) {
00074 fprintf(stderr, "Cant make year directory '%s'\n", outfilename);
00075 return 1;
00076 }
00077 }
00078 sprintf(outfilename, "%s/%4.4d/%2.2d", outpathbase, yr, mo);
00079 if (stat(outfilename, &sbuf)) {
00080 if (mkdir(outfilename, 0777)) {
00081 fprintf(stderr, "Cant make month directory '%s'\n", outfilename);
00082 return 1;
00083 }
00084 }
00085 sprintf(outfilename, "%s/%4.4d/%2.2d/%2.2d", outpathbase, yr, mo, da);
00086 if (stat(outfilename, &sbuf)) {
00087 if (mkdir(outfilename, 0777)) {
00088 fprintf(stderr, "Cant make day directory '%s'\n", outfilename);
00089 return 1;
00090 }
00091 }
00092 sprintf(outfilename, "%s/%4.4d/%2.2d/%2.2d/H%2.2d00",
00093 outpathbase, yr, mo, da, hr);
00094 if (stat(outfilename, &sbuf)) {
00095 if (mkdir(outfilename, 0777)) {
00096 fprintf(stderr, "Cant make day directory '%s'\n", outfilename);
00097 return 1;
00098 }
00099 }
00100 return 0;
00101 }
00102
00103 int DoIt ()
00104 {
00105 int irec, iseg, nrecs, nsegs, status, n_used=1, wl, explim, fsn, quality;
00106 int file, drms, cmdexp, nogood, oy, om, od, oh, on;
00107 char c1, c2, c3, *date_str, *dsinp, *dsout, now_str[100], *outpathbase=NULL;
00108 char querystr[512], outpath[512], outfilename[512], *outcparms="compress";
00109 float crpix1, crpix2, cdelt1, cdelt2, crota2, x0, y0;
00110 float mag = 1.0, dx, dy;
00111 double gaex_obs, gaey_obs, gaez_obs, haex_obs, haey_obs, haez_obs, tr_step;
00112 long long tr_index;
00113 FILE *fimg_time;
00114 TIME t_obs, t_rec, tr_epoch, t_obs_alt;
00115 DRMS_Record_t *inprec, *outrec;
00116 DRMS_RecordSet_t *inprs, *aecrs;
00117 DRMS_Keyword_t *inpkey = NULL, *outkey = NULL;
00118 DRMS_Array_t *inparr=NULL, *outarr=NULL;
00119 DRMS_Segment_t *inpseg, *outseg;
00120 HIterator_t *hit = NULL;
00121 if (nice_intro(0)) return(0);
00122
00123 dsinp = strdup(cmdparams_get_str(&cmdparams, "dsinp", NULL));
00124 dsout = strdup(cmdparams_get_str(&cmdparams, "dsout", NULL));
00125 outpathbase = strdup(cmdparams_get_str(&cmdparams, "pathout", NULL));
00126 drms = cmdparams_get_int(&cmdparams, "drms", NULL);
00127 file = cmdparams_get_int(&cmdparams, "file", NULL);
00128 if (strcmp(dsinp, NOT_SPECIFIED)==0) DIE("dsinp argument is required");
00129 if (strcmp(dsout,NOT_SPECIFIED)==0) DIE("dsout argument is required");
00130 rescale = cmdparams_get_int(&cmdparams, "rescale", NULL);
00131 scale_to = cmdparams_get_float(&cmdparams, "scale_to", NULL);
00132 regridtype = cmdparams_get_int(&cmdparams, "regrid", NULL);
00133 do_stretchmarks = cmdparams_get_int(&cmdparams, "do_stretchmarks", NULL);
00134
00135 inprs = drms_open_records(drms_env, dsinp, &status);
00136 if (status) DIE("cant open recordset query");
00137 drms_stage_records(inprs, 1, 0);
00138 nrecs = inprs->n;
00139 printf("%d records\n", nrecs);
00140 for (irec=0; irec<nrecs; irec++) {
00141 int yr, mo, da, hr, mn, sc;
00142 inprec = inprs->records[irec];
00143 wl = drms_getkey_int(inprec, "WAVELNTH", &status);
00144 if (status) DIE("WAVELNTH not found!");
00145 switch (wl) {
00146 case 94:
00147 case 131:
00148 case 211:
00149 case 304:
00150 case 335:
00151 case 1600: explim = 2800; break;
00152 case 171:
00153 case 193: explim = 1931; break;
00154 case 1700: explim = 966; break;
00155 case 4500: explim = 483; break;
00156 default: DIE("Bad wavelength");
00157 }
00158 imtype = drms_getkey_string(inprec, "IMG_TYPE", &status);
00159 if (status) DIE("cant get IMG_TYPE");
00160 fsn = drms_getkey_int(inprec, "FSN", &status);
00161 if (status) DIE("cant get FSN");
00162 cmdexp = drms_getkey_int(inprec, "AIMGSHCE", &status);
00163 if (status) DIE("cant get commanded exposure AIMGSHCE");
00164 quality = drms_getkey_int(inprec, "QUALITY", &status);
00165 if (status) DIE("cant get QUALITY");
00166 t_obs = drms_getkey_time(inprec, "T_OBS", &status);
00167 if (status) DIE("T_OBS not found!");
00168 if (cmdexp < explim) {
00169 char *pos;
00170 strncpy(querystr, dsinp, 512);
00171 t_obs_alt = t_obs - 30.0;
00172 pos = index(querystr, '[');
00173 snprintf(pos, 512, "[? T_REC between %10f and %10f ?][? WAVELNTH = %d ?]",
00174 t_obs - 30.0, t_obs + 30.0, wl);
00175 nogood = 1;
00176 printf("%s\n", querystr);
00177 } else nogood = 0;
00178 while (nogood) {
00179 nogood = 0;
00180 }
00181 if (nogood) continue;
00182 date_str = drms_getkey_string(inprec, "T_OBS", &status);
00183 sscanf(date_str, "%d%c%d%c%d%c%d:%d:%d",
00184 &yr, &c1, &mo, &c2, &da, &c3, &hr, &mn, &sc);
00185 if (file) {
00186 if (check_outpath(yr, mo, da, hr, outpathbase)) DIE("Cant write to out");
00187 }
00188 outrec = drms_create_record(drms_env, dsout, DRMS_PERMANENT, &status);
00189 if (status) DIE("cant create recordset");
00190 status = drms_copykeys(outrec, inprec, 0, kDRMS_KeyClass_Explicit);
00191 if (status) DIE("Error in drms_copykeys()");
00192
00193 sprint_time_ISO(now_str, CURRENT_SYSTEM_TIME);
00194 drms_setkey_string(outrec, "DATE", now_str);
00195 t_rec = t_obs + 30.0;
00196 if ( 0 == drms_setkey_time(outrec, "T_REC", t_rec)) {
00197 date_str = drms_getkey_string(outrec, "T_REC", &status);
00198 t_rec = (long long) (t_obs+0.5);
00199 drms_setkey_time(outrec, "T_REC", t_rec);
00200 } else {
00201 tr_step = 120.0; tr_epoch = 0;
00202 tr_index = (t_obs - tr_epoch)/tr_step;
00203 t_rec = tr_epoch + tr_index*tr_step;
00204 }
00205 sscanf(date_str, "%d%c%d%c%d%c%d:%d",
00206 &yr, &c1, &mo, &c2, &da, &c3, &hr, &mn);
00207 if (file) {
00208 if (check_outpath(yr, mo, da, hr, outpathbase)) DIE("Cant write to out");
00209 }
00210 crpix1 = drms_getkey_float(inprec, "CRPIX1", &status);
00211 if (status) DIE("CRPIX1 not found!");
00212 crpix2 = drms_getkey_float(inprec, "CRPIX2", &status);
00213 if (status) DIE("CRPIX2 not found!");
00214 cdelt1 = drms_getkey_float(inprec, "CDELT1", &status);
00215 if (status) DIE("CDELT1 not found!");
00216 cdelt2 = drms_getkey_float(inprec, "CDELT2", &status);
00217 if (status) DIE("CDELT2 not found!");
00218 if (rescale) {
00219 mag = fabs(cdelt1 / scale_to);
00220 cdelt1 = 4.0*cdelt1/mag;
00221 cdelt2 = 4.0*cdelt2/mag;
00222 drms_setkey_float(outrec, "CDELT1", cdelt1);
00223 drms_setkey_float(outrec, "CDELT2", cdelt2);
00224 }
00225 crota2 = drms_getkey_float(inprec, "CROTA2", &status);
00226 if (status) DIE("CROTA2 not found!");
00227 x0 = drms_getkey_float(inprec, "X0", &status);
00228
00229 y0 = drms_getkey_float(inprec, "Y0", &status);
00230
00231 if (hit) hiter_destroy(&hit);
00232 nsegs = hcon_size(&inprec->segments);
00233 for (iseg=0; iseg<1; iseg++) {
00234 void *output_array = NULL;
00235 int i, j, n, m, dtyp, nx, ny, outaxis[2];
00236 char *filename = NULL, *inpsegname = NULL;
00237 inpseg = drms_segment_lookupnum(inprec, iseg);
00238 inpsegname = inpseg->info->name;
00239 filename = inpseg->filename;
00240 if (0 == iseg) {
00241 outseg = drms_segment_lookupnum(outrec, iseg);
00242 if (!outseg) DIE("Cant get output segment");
00243 inparr = drms_segment_read(inpseg, DRMS_TYPE_FLOAT, &status);
00244 if (status) DIE("drms_segment_read failed!");
00245 dtyp = 3;
00246 n = inparr->axis[0]; m = inparr->axis[1];
00247
00248 dx = (n + 1.0)*0.5 - crpix1; dy = (m + 1.0)*0.5 - crpix2;
00249
00250 status = image_magrotate( inparr->data, n, m, dtyp, crota2,
00251 mag, dx, dy, &output_array, &nx, &ny,
00252 regridtype, do_stretchmarks);
00253 if (status) DIE("image_magrotate failed!");
00254 outaxis[0] = 1024; outaxis[1] = 1024;
00255 outarr = drms_array_create(DRMS_TYPE_FLOAT, 2, outaxis,
00256 NULL, &status);
00257 if (status) DIE("drms_array_create failed!");
00258
00259 outarr->bscale = 0.0625;
00260 outarr->bzero = 0.0;
00261 for (i=0; i<1024; i++) {
00262 int k0 = 4*i;
00263 for (j=0; j<1024; j++) {
00264 int k, l, nn=0, l0=4*j;
00265 float sum=0.0;
00266 for (k=k0; k<k0+4; k++) for (l=l0; l<l0+4; l++) {
00267 int ndx = 4096*k + l;
00268 float pv = *((float *)(output_array)+ndx);
00269 if (!isnan(pv)) {sum += pv; nn++; }
00270 }
00271 *((float *)(outarr->data)+(1024*i+j)) = sum/nn;
00272 }
00273 }
00274 drms_setkey_float(outrec, "CROTA2", 0.0);
00275 drms_setkey_float(outrec, "CRPIX1", (1024 + 1.0)*0.5);
00276 drms_setkey_float(outrec, "CRPIX2", (1024 + 1.0)*0.5);
00277 drms_setkey_float(outrec, "X0", 0.0);
00278 drms_setkey_float(outrec, "Y0", 0.0);
00279 drms_segment_write(outseg, outarr, 0);
00280 if (file) {
00281 struct stat sbuf;
00282 sprintf(outpath, "%s/%4.4d/%2.2d/%2.2d/H%2.2d00",
00283 outpathbase, yr, mo, da, hr);
00284 sprintf(outfilename,
00285 "%s/AIA%4.4d%2.2d%2.2d_%2.2d%2.2d_%4.4d.fits",
00286 outpath, yr, mo, da, hr, mn, wl);
00287 oy = yr; om = mo; od = da; oh = hr; on = mn;
00288 if (stat(outfilename, &sbuf)) {
00289 outcparms = "compress";
00290 if (DRMS_SUCCESS != fitsexport_export_tofile(outseg, outcparms,
00291 outfilename, NULL, NULL)) DIE("cant write FITS file");
00292 }
00293 }
00294 if (inparr) drms_free_array(inparr);
00295 if (output_array) free(output_array);
00296 if (outarr) drms_free_array(outarr);
00297 }
00298 }
00299 if (drms) drms_close_record(outrec, DRMS_INSERT_RECORD);
00300 else drms_close_record(outrec, DRMS_FREE_RECORD);
00301 }
00302 if (wl > 1500) return 0;
00303 sprintf(outpath, "%s/image_times", outpathbase);
00304 if (fimg_time = fopen(outpath, "w")) {
00305 fprintf(fimg_time, "Time %d%2.2d%2.2d_%2.2d%2.2d00\n",
00306 oy, om, od, oh, on);
00307 fprintf(fimg_time, "synoptic http://jsoc.stanford.edu/data/aia/synoptic/");
00308 fprintf(fimg_time, "%d/%2.2d/%2.2d/H%2.2d00/AIA%d%2.2d%2.2d_%2.2d%2.2d00\n"
00309 , oy, om, od, oh, oy, om, od, oh, on);
00310 fclose(fimg_time);
00311 } else fprintf(stderr, "Can't open image_times file.\n");
00312 sprintf(outpath, "%s/image_times.json", outpathbase);
00313 if (fimg_time = fopen(outpath, "w")) {
00314 fprintf(fimg_time, "{\"first\":\"20100513_000000\",\"last\":\"");
00315 fprintf(fimg_time, "%d%2.2d%2.2d_%2.2d%2.2d00\"}\n", oy, om, od, oh, on);
00316 fclose(fimg_time);
00317 } else fprintf(stderr, "Can't open image_times.json file.\n");
00318
00319 return 0;
00320 }
