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