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, "pathout", "/scr21/jsoc/aia/synoptic/mostrecent", "out path"},
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 TIME tbeg;
00030 char *module_name = "aia_most_recent", *imtype;
00031 int despike = 1, rescale = 1, regridtype = 1, verbose = 0, do_stretchmarks = 0;
00032 float scale_to = 0.6, scale_by;
00033
00034 int nice_intro(int help)
00035 {
00036 int usage = cmdparams_get_int(&cmdparams, "h", NULL) != 0;
00037 verbose = cmdparams_get_int(&cmdparams, "v", NULL) != 0;
00038 if (usage || help) {
00039 printf("aia_most_recent {-h} {-v} dsinp=series_record_spec dsout=output_series\n"
00040 " -h: print this message\n"
00041 " -v: verbose\n"
00042 "pathout=<output path>, default=/scr21/jsoc/aia/synoptic/mostrecent\n"
00043 "rescale=0, do not rescale to fixed plate scale, default=1\n"
00044 "regrid=0, nearest neighbor, default=1, bicubic\n"
00045 "scale_to=0.6, rescale to fixed plate scale, default=0.6\n"
00046 "do_stretchmarks=0, fill in empty pixels created, default=0\n"
00047 "dsinp=<recordset query> as <series>{[record specifier]} - required\n"
00048 "dsout=<series> - required if output to drms series\n");
00049 return(1);
00050 }
00051 return(0);
00052 }
00053
00054 void sprint_time_ISO (char *tstring, TIME t)
00055 {
00056 sprint_at(tstring,t);
00057 tstring[4] = tstring[7] = '-';
00058 tstring[10] = 'T';
00059 tstring[19] = '\0';
00060 }
00061
00062 int DoIt ()
00063 {
00064 int got_all_wl=0,iwl, irec, iseg, nrecs, nsegs, status, n_used=1, wl;
00065 int file, drms, oy, om, od, oh, on;
00066 const int wls[10] = { 94, 131, 171, 193, 211, 304, 335, 1600, 1700, 4500 };
00067 char c1, c2, c3, *date_str, *dsinp, *dsout, now_str[100], *outpathbase=NULL;
00068 char linkname[512], outpath[512], outfilename[512], *outcparms="compress";
00069 float crpix1, crpix2, cdelt1, cdelt2, crota2, x0, y0;
00070 float mag = 1.0, dx, dy;
00071 double gaex_obs, gaey_obs, gaez_obs, haex_obs, haey_obs, haez_obs, tr_step;
00072 long long tr_index;
00073 FILE *fimg_time;
00074 TIME t_obs, t_rec, tr_epoch, t0=0;
00075 DRMS_Record_t *inprec, *outrec;
00076 DRMS_RecordSet_t *inprs;
00077 DRMS_Keyword_t *inpkey = NULL, *outkey = NULL;
00078 DRMS_Array_t *inparr=NULL, *outarr=NULL;
00079 DRMS_Segment_t *inpseg, *outseg;
00080 char *selstr = "fsn = (select max(fsn) from ", query[512];
00081 char *whrstr = " where wavelnth=";
00082 if (nice_intro(0)) return(0);
00083
00084 dsinp = strdup(cmdparams_get_str(&cmdparams, "dsinp", NULL));
00085 dsout = strdup(cmdparams_get_str(&cmdparams, "dsout", NULL));
00086 outpathbase = strdup(cmdparams_get_str(&cmdparams, "pathout", NULL));
00087 if (strcmp(dsinp, NOT_SPECIFIED)==0) DIE("dsinp argument is required");
00088 if (strcmp(dsout,NOT_SPECIFIED)==0) DIE("dsout argument is required");
00089 rescale = cmdparams_get_int(&cmdparams, "rescale", NULL);
00090 scale_to = cmdparams_get_float(&cmdparams, "scale_to", NULL);
00091 regridtype = cmdparams_get_int(&cmdparams, "regrid", NULL);
00092 do_stretchmarks = cmdparams_get_int(&cmdparams, "do_stretchmarks", NULL);
00093
00094 for (iwl=0; iwl<10; iwl++) {
00095 int yr, mo, da, hr, mn, sc;
00096 void *output_array = NULL;
00097 int i, j, n, m, dtyp, nx, ny, outaxis[2];
00098 char *filename = NULL, *inpsegname = NULL;
00099 sprintf(query, "%s[? %s%s%s%d", dsinp, selstr, dsinp, whrstr, wls[iwl]);
00100 strcat(query, " and img_type='LIGHT') ?]");
00101 inprs = drms_open_records(drms_env, query, &status);
00102 if (status) DIE("cant open recordset query");
00103 nrecs = inprs->n;
00104 inprec = inprs->records[0];
00105 imtype = drms_getkey_string(inprec, "IMG_TYPE", &status);
00106 t_obs = drms_getkey_time(inprec, "T_OBS", &status);
00107 if (status) DIE("T_OBS not found!");
00108 date_str = drms_getkey_string(inprec, "T_OBS", &status);
00109 sscanf(date_str, "%d%c%d%c%d%c%d:%d:%d",
00110 &yr, &c1, &mo, &c2, &da, &c3, &hr, &mn, &sc);
00111 if (t_obs > t0) {
00112 oy = yr; om = mo; od = da; oh = hr; on = mn;
00113 t0 = t_obs;
00114 }
00115 wl = drms_getkey_int(inprec, "WAVELNTH", &status);
00116 if (status) DIE("WAVELNTH not found!");
00117 outrec = drms_create_record(drms_env, dsout, DRMS_PERMANENT, &status);
00118 if (status) DIE("cant create recordset");
00119 status = drms_copykeys(outrec, inprec, 0, kDRMS_KeyClass_Explicit);
00120 if (status) DIE("Error in drms_copykeys()");
00121
00122 sprint_time_ISO(now_str, CURRENT_SYSTEM_TIME);
00123 drms_setkey_string(outrec, "DATE", now_str);
00124 if ( 0 == drms_setkey_time(outrec, "T_REC", t_obs)) {
00125 tr_index = drms_getkey_longlong(outrec, "T_REC_index", &status);
00126 if (status) DIE("T_REC_index not found!");
00127 tr_step = drms_getkey_double(outrec, "T_REC_step", &status);
00128 if (status) DIE("T_REC_step not found!");
00129 tr_epoch = drms_getkey_time(outrec, "T_REC_epoch", &status);
00130 if (status) DIE("T_REC_epoch not found!");
00131 t_rec = tr_epoch + tr_index*tr_step;
00132 drms_setkey_time(outrec, "T_REC", t_rec);
00133 date_str = drms_getkey_string(outrec, "T_REC", &status);
00134 } else {
00135 tr_step = 180.0; tr_epoch = 0;
00136 tr_index = (t_obs - tr_epoch)/tr_step;
00137 t_rec = tr_epoch + tr_index*tr_step;
00138 }
00139 sscanf(date_str, "%d%c%d%c%d%c%d:%d",
00140 &yr, &c1, &mo, &c2, &da, &c3, &hr, &mn);
00141 sc = 0;
00142 crpix1 = drms_getkey_float(inprec, "CRPIX1", &status);
00143 if (status) DIE("CRPIX1 not found!");
00144 crpix2 = drms_getkey_float(inprec, "CRPIX2", &status);
00145 if (status) DIE("CRPIX2 not found!");
00146 cdelt1 = drms_getkey_float(inprec, "CDELT1", &status);
00147 if (status) DIE("CDELT1 not found!");
00148 cdelt2 = drms_getkey_float(inprec, "CDELT2", &status);
00149 if (status) DIE("CDELT2 not found!");
00150 if (rescale) {
00151 mag = fabs(cdelt1 / scale_to);
00152 cdelt1 = 4.0*cdelt1/mag;
00153 cdelt2 = 4.0*cdelt2/mag;
00154 drms_setkey_float(outrec, "CDELT1", cdelt1);
00155 drms_setkey_float(outrec, "CDELT2", cdelt2);
00156 }
00157 crota2 = drms_getkey_float(inprec, "CROTA2", &status);
00158 if (status) DIE("CROTA2 not found!");
00159 x0 = drms_getkey_float(inprec, "X0", &status);
00160
00161 y0 = drms_getkey_float(inprec, "Y0", &status);
00162
00163 nsegs = hcon_size(&inprec->segments);
00164 iseg=0;
00165 inpseg = drms_segment_lookupnum(inprec, iseg);
00166 inpsegname = inpseg->info->name;
00167 filename = inpseg->filename;
00168 outseg = drms_segment_lookupnum(outrec, iseg);
00169 if (!outseg) DIE("Cant get output segment");
00170 inparr = drms_segment_read(inpseg, DRMS_TYPE_FLOAT, &status);
00171 if (status) DIE("drms_segment_read failed!");
00172 dtyp = 3;
00173 n = inparr->axis[0]; m = inparr->axis[1];
00174
00175 dx = (n + 1.0)*0.5 - crpix1; dy = (m + 1.0)*0.5 - crpix2;
00176
00177 status = image_magrotate( inparr->data, n, m, dtyp, crota2,
00178 mag, dx, dy, &output_array, &nx, &ny,
00179 regridtype, do_stretchmarks);
00180 if (status) DIE("image_magrotate failed!");
00181 outaxis[0] = 1024; outaxis[1] = 1024;
00182 outarr = drms_array_create(DRMS_TYPE_FLOAT, 2, outaxis,
00183 NULL, &status);
00184 if (status) DIE("drms_array_create failed!");
00185
00186 outarr->bscale = 0.0625;
00187 outarr->bzero = 0.0;
00188 for (i=0; i<1024; i++) {
00189 int k0 = 4*i;
00190 for (j=0; j<1024; j++) {
00191 int k, l, nn=0, l0=4*j;
00192 float sum=0.0;
00193 for (k=k0; k<k0+4; k++) for (l=l0; l<l0+4; l++) {
00194 int ndx = 4096*k + l;
00195 float pv = *((float *)(output_array)+ndx);
00196 if (!isnan(pv)) {sum += pv; nn++; }
00197 }
00198 *((float *)(outarr->data)+(1024*i+j)) = sum/nn;
00199 }
00200 }
00201 drms_setkey_float(outrec, "CROTA2", 0.0);
00202 drms_setkey_float(outrec, "CRPIX1", (1024 + 1.0)*0.5);
00203 drms_setkey_float(outrec, "CRPIX2", (1024 + 1.0)*0.5);
00204 drms_setkey_float(outrec, "X0", 0.0);
00205 drms_setkey_float(outrec, "Y0", 0.0);
00206 drms_segment_write(outseg, outarr, 0);
00207 sprintf(outfilename, "%s/AIAsynoptic%4.4d.fits", outpathbase, wl);
00208 outcparms = "";
00209 if (DRMS_SUCCESS != fitsexport_export_tofile(outseg, outcparms,
00210 outfilename, NULL, NULL)) DIE("cant write FITS file");
00211 if (inparr) drms_free_array(inparr);
00212 if (output_array) free(output_array);
00213 if (outarr) drms_free_array(outarr);
00214 drms_close_record(outrec, DRMS_FREE_RECORD);
00215 }
00216 sprintf(outfilename, "%s/image_times", outpathbase);
00217 if (fimg_time = fopen(outfilename, "w")) {
00218 fprintf(fimg_time, "Time %d%2.2d%2.2d_%2.2d%2.2d00\n",
00219 oy, om, od, oh, on);
00220 fprintf(fimg_time,
00221 "mostrecent http://jsoc.stanford.edu/data/aia/synoptic/mostrecent/");
00222 fprintf(fimg_time, "AIAsynoptic\n");
00223 fclose(fimg_time);
00224 } else fprintf(stderr, "Can't open image_times file.\n");
00225
00226 return 0;
00227 }
