108 mbobra 1.1 // Mapping method
109 enum projection {
110 carree,
111 cassini,
112 mercator,
113 cyleqa,
114 sineqa,
115 gnomonic,
116 postel,
117 stereographic,
118 orthographic,
119 lambert
120 };
121
122 // WSC code
123 char *wcsCode[] = {"CAR", "CAS", "MER", "CEA", "GLS", "TAN", "ARC", "STG",
124 "SIN", "ZEA"};
125
126 // Ephemeris information
127 struct ephemeris {
128 double disk_lonc, disk_latc;
129 mbobra 1.1 double disk_xc, disk_yc;
130 double rSun, asd, pa;
131 };
132
133 // Mapping information
134 struct mapInfo {
135 float xc, yc; // reference point: center
136 int nrow, ncol; // size
137 float xscale, yscale; // scale
138 int nbin;
139 enum projection proj; // projection method
140 struct ephemeris ephem; // ephemeris info
141 float *xi_out, *zeta_out; // coordinate on full disk image to sample at
142 };
143
144 /* ========================================================================================================== */
145
146 /* Get all input data series */
147 int getInputRS(DRMS_RecordSet_t **mharpRS_ptr, char *mharpQuery);
148
149 /* Get other data series */
150 mbobra 1.1 int getInputRS_aux(DRMS_RecordSet_t **inRS_ptr, char *inQuery, DRMS_RecordSet_t *harpRS);
151
152 /* Find record from record set with given T_rec */
153 int getInputRec_aux(DRMS_Record_t **inRec_ptr, DRMS_RecordSet_t *inRS, TIME trec);
154
155 /* Create CEA record */
156 int createCeaRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *contRec, DRMS_Record_t *sharpRec, struct swIndex *swKeys_ptr);
157
158 /* Mapping single segment, wrapper */
159 int mapScaler(DRMS_Record_t *sharpRec, DRMS_Record_t *inRec, DRMS_Record_t *harpRec,
160 struct mapInfo *mInfo, char *segName);
161
162 /* Determine reference point coordinate and patch size according to input */
163 int findPosition(DRMS_Record_t *inRec, struct mapInfo *mInfo);
164
165 /* Get ephemeris information */
166 int getEphemeris(DRMS_Record_t *inRec, struct ephemeris *ephem);
167
168 /* Compute the coordinates at which the full disk image is sampled */
169 void findCoord(struct mapInfo *mInfo);
170
171 mbobra 1.1 /* Mapping function */
172 int performSampling(float *outData, float *inData, struct mapInfo *mInfo, int interpOpt);
173
174 // ===================
175
176 /* Create Cutout record */
177 int createCutRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *contRec, DRMS_Record_t *sharpRec, struct swIndex *swKeys_ptr);
178
179 /* Get cutout and write segment */
180 int writeCutout(DRMS_Record_t *outRec, DRMS_Record_t *inRec, DRMS_Record_t *harpRec, char *SegName);
181
182 // ===================
183
184 /* Compute space weather indices */
185 void computeSWIndex(struct swIndex *swKeys_ptr, DRMS_Record_t *inRec, struct mapInfo *mInfo);
186
187 /* Set space weather indices */
188 void setSWIndex(DRMS_Record_t *outRec, struct swIndex *swKeys_ptr);
189
190 /* Set all keywords */
191 void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *mharpRec, struct mapInfo *mInfo);
192 mbobra 1.1
193 // ===================
194
195 /* Nearest neighbor interpolation */
196 float nnb (float *f, int nx, int ny, double x, double y);
197
198 /* Wrapper for Jesper's rebin code */
199 void frebin (float *image_in, float *image_out, int nx, int ny, int nbin, int gauss);
200
201 /* ========================================================================================================== */
202
203 /* Cutout segment names, input identical to output */
204 char *MharpSegs[] = {"magnetogram", "bitmap"};
205 char *CutSegs[] = {"magnetogram", "bitmap", "continuum"};
206 char *CEASegs[] = {"magnetogram", "bitmap", "continuum"};
207 // For BUNIT
208 char *CutBunits[] = {"Mx/cm^2", " ", "DN/s"};
209 char *CEABunits[] = {"Mx/cm^2", " ", "DN/s"};
210 /* ========================================================================================================== */
211
212 char *module_name = "smarp";
213 mbobra 1.1
214 ModuleArgs_t module_args[] =
215 {
216 {ARG_STRING, "mharp", kNotSpecified, "Input Mharp series."},
217 {ARG_STRING, "cont", kNotSpecified, "Input Continuum series."},
218 {ARG_STRING, "sharp_cea", kNotSpecified, "Output Sharp CEA series."},
219 {ARG_STRING, "sharp_cut", kNotSpecified, "Output Sharp cutout series."},
220 {ARG_END}
221 };
222
223 int DoIt(void)
224 {
225 int errbufstat = setvbuf(stderr, NULL, _IONBF, BUFSIZ);
226 int outbufstat = setvbuf(stdout, NULL, _IONBF, BUFSIZ);
227 int status = DRMS_SUCCESS;
228 int nrecs, irec;
229 char *mharpQuery;
230 char *contQuery;
231 char *sharpCeaQuery, *sharpCutQuery;
232 DRMS_RecordSet_t *mharpRS = NULL;
233 DRMS_RecordSet_t *contRS = NULL;
234 mbobra 1.1
235 /* Get parameters */
236
237 mharpQuery = (char *) params_get_str(&cmdparams, "mharp");
238 sharpCeaQuery = (char *) params_get_str(&cmdparams, "sharp_cea");
239 sharpCutQuery = (char *) params_get_str(&cmdparams, "sharp_cut");
240 contQuery = (char *) params_get_str(&cmdparams, "cont");
241
242 /* Get input data, check everything */
243 if (getInputRS(&mharpRS, mharpQuery))
244 DIE("Input harp data error.");
245 nrecs = mharpRS->n;
246
247 if (getInputRS_aux(&contRS, contQuery, mharpRS))
248 DIE("Input continuum data error.");
249
250 /* Start */
251
252 printf("==============\nStart. %d image(s) in total.\n", nrecs);
253
254 for (irec = 0; irec < nrecs; irec++) {
255 mbobra 1.1
256 /* Records in work */
257
258 DRMS_Record_t *mharpRec = NULL;
259 DRMS_Record_t *contRec = NULL;
260
261 mharpRec = mharpRS->records[irec];
262 TIME trec = drms_getkey_time(mharpRec, "T_REC", &status);
263
264 struct swIndex swKeys;
265
266 if (getInputRec_aux(&contRec, contRS, trec)) {
267 printf("Fetching Continuum failed, image #%d skipped.\n", irec);
268 continue;
269 }
270
271 printf("Obtained all the data \n");
272
273 /* Create CEA record */
274
275 DRMS_Record_t *sharpCeaRec = drms_create_record(drms_env, sharpCeaQuery, DRMS_PERMANENT, &status);
276 mbobra 1.1 if (status) { // if failed
277 printf("Creating CEA failed, image #%d skipped.\n", irec);
278 continue;
279 }
280 if (createCeaRecord(mharpRec, contRec, sharpCeaRec, &swKeys)) { // do the work
281 printf("Creating CEA failed, image #%d skipped.\n", irec);
282 drms_close_record(sharpCeaRec, DRMS_FREE_RECORD);
283 continue;
284 } // swKeys updated here
285
286 drms_close_record(sharpCeaRec, DRMS_INSERT_RECORD);
287
288 printf("Created CEA record \n");
289
290 /* Create Cutout record */
291
292 DRMS_Record_t *sharpCutRec = drms_create_record(drms_env, sharpCutQuery, DRMS_PERMANENT, &status);
293 if (status) { // if failed
294 printf("Creating cutout failed, image #%d skipped.\n", irec);
295 continue;
296 }
297 mbobra 1.1
298 if (createCutRecord(mharpRec, contRec, sharpCutRec, &swKeys)) { // do the work
299 printf("Creating cutout failed, image #%d skipped.\n", irec);
300 drms_close_record(sharpCutRec, DRMS_FREE_RECORD);
301 continue;
302 } // swKeys used here
303 drms_close_record(sharpCutRec, DRMS_INSERT_RECORD);
304 printf("Created CUT record \n");
305 /* Done */
306
307 printf("Image #%d done.\n", irec);
308
309 } // irec
310
311
312 drms_close_records(mharpRS, DRMS_FREE_RECORD);
313 drms_close_records(contRS, DRMS_FREE_RECORD);
314
315 return 0;
316
317 } // DoIt
318 mbobra 1.1
319 // ===================================================================
320 // ===================================================================
321 // ===================================================================
322
323 /*
324 * Get input data series, including mHarp and bharp
325 * Need all records to match, otherwise quit
326 *
327 */
328
329 int getInputRS(DRMS_RecordSet_t **mharpRS_ptr, char *mharpQuery)
330 {
331 int status = 0;
332 *mharpRS_ptr = drms_open_records(drms_env, mharpQuery, &status);
333 if (status || (*mharpRS_ptr)->n == 0) return 1;
334 return 0;
335 }
336
337 /*
338 * Get other data series, check all T_REC are available
339 mbobra 1.1 */
340
341 int getInputRS_aux(DRMS_RecordSet_t **inRS_ptr, char *inQuery, DRMS_RecordSet_t *harpRS)
342 {
343
344 int status = 0;
345
346 *inRS_ptr = drms_open_records(drms_env, inQuery, &status);
347 if (status || (*inRS_ptr)->n == 0) return status;
348
349 // Check if all T_rec are available, need to match both ways
350 int n = harpRS->n, n0 = (*inRS_ptr)->n;
351
352 for (int i0 = 0; i0 < n0; i0++) {
353 DRMS_Record_t *inRec = (*inRS_ptr)->records[i0];
354 TIME trec0 = drms_getkey_time(inRec, "T_REC", &status);
355 TIME trec = 0;
356 for (int i = 0; i < n; i++) {
357 DRMS_Record_t *harpRec = harpRS->records[i];
358 trec = drms_getkey_time(harpRec, "T_REC", &status);
359 if (fabs(trec0 - trec) < 10) break;
360 mbobra 1.1 }
361 if (fabs(trec0 - trec) >= 10) return 1;
362 }
363
364 for (int i = 0; i < n; i++) {
365 DRMS_Record_t *harpRec = harpRS->records[i];
366 TIME trec = drms_getkey_time(harpRec, "T_REC", &status);
367 TIME trec0 = 0;
368 for (int i0 = 0; i0 < n0; i0++) {
369 DRMS_Record_t *inRec = (*inRS_ptr)->records[i0];
370 trec0 = drms_getkey_time(inRec, "T_REC", &status);
371 if (fabs(trec0 - trec) < 10) break;
372 }
373 if (fabs(trec0 - trec) >= 10) return 1;
374 }
375
376 return 0;
377
378 }
379
380 /*
381 mbobra 1.1 * Find record from record set with given T_rec
382 */
383
384 int getInputRec_aux(DRMS_Record_t **inRec_ptr, DRMS_RecordSet_t *inRS, TIME trec)
385 {
386
387 int status = 0;
388
389 int n = inRS->n;
390 for (int i = 0; i < n; i++) {
391 *inRec_ptr = inRS->records[i];
392 TIME trec0 = drms_getkey_time((*inRec_ptr), "T_REC", &status);
393 if (fabs(trec0 - trec) < 10) return 0;
394 }
395
396 return 1;
397
398 }
399
400
401
402 mbobra 1.1
403 /*
404 * Create CEA record: top level subroutine
405 * Also compute all the space weather keywords here
406 */
407
408 int createCeaRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *contRec, DRMS_Record_t *sharpRec, struct swIndex *swKeys_ptr)
409 {
410
411 int status = 0;
412 DRMS_Segment_t *inSeg;
413 DRMS_Array_t *inArray;
414 int val;
415
416 struct mapInfo mInfo;
417 mInfo.proj = (enum projection) cyleqa; // projection method
418 mInfo.xscale = XSCALE;
419 mInfo.yscale = YSCALE;
420
421 int ncol0, nrow0; // oversampled map size
422
423 mbobra 1.1 // Get ephemeris
424
425 if (getEphemeris(mharpRec, &(mInfo.ephem))) {
426 SHOW("CEA: get ephemeris error\n");
427 return 1;
428 }
429
430 // Find position
431
432 if (findPosition(mharpRec, &mInfo)) {
433 SHOW("CEA: find position error\n");
434 return 1;
435 }
436
437 // ========================================
438 // Do this for all bitmaps, Aug 12 2013 XS
439 // ========================================
440
441 mInfo.nbin = 1; // for bitmaps. suppress anti-aliasing
442 ncol0 = mInfo.ncol;
443 nrow0 = mInfo.nrow;
444 mbobra 1.1
445 mInfo.xi_out = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
446 mInfo.zeta_out = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
447
448 findCoord(&mInfo); // compute it here so it could be shared by the following 4 functions
449
450 if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "bitmap")) {
451 SHOW("CEA: mapping bitmap error\n");
452 return 1;
453 }
454 printf("Bitmap mapping done.\n");
455
456 free(mInfo.xi_out);
457 free(mInfo.zeta_out);
458
459 // ========================================
460 // Do this again for floats, Aug 12 2013 XS
461 // ========================================
462 // Create xi_out, zeta_out array in mInfo:
463 // Coordinates to sample in original full disk image
464
465 mbobra 1.1 mInfo.nbin = NBIN;
466 ncol0 = mInfo.ncol * mInfo.nbin + (mInfo.nbin / 2) * 2; // pad with nbin/2 on edge to avoid NAN
467 nrow0 = mInfo.nrow * mInfo.nbin + (mInfo.nbin / 2) * 2;
468
469 mInfo.xi_out = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
470 mInfo.zeta_out = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
471
472 findCoord(&mInfo); // compute it here so it could be shared by the following 4 functions
473
474 // Mapping single segment: Mharp, etc.
475
476 if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "magnetogram")) {
477 SHOW("CEA: mapping magnetogram error\n");
478 return 1;
479 }
480 printf("Magnetogram mapping done.\n");
481
482 if (mapScaler(sharpRec, contRec, mharpRec, &mInfo, "continuum")) {
483 SHOW("CEA: mapping continuum error\n");
484 return 1;
485 }
486 mbobra 1.1 printf("Intensitygram mapping done.\n");
487
488 // Keywords & Links
489 copy_patch_keys(mharpRec, sharpRec);
490 copy_geo_keys(mharpRec, sharpRec);
491 // rename HARPNUM to TARPNUM
492 val = drms_getkey_double(mharpRec, "HARPNUM", &status);
493 drms_setkey_double(sharpRec, "TARPNUM", val);
494 // copy everything else
495 drms_copykey(sharpRec, mharpRec, "T_REC");
496 drms_copykey(sharpRec, mharpRec, "CDELT1");
497 drms_copykey(sharpRec, mharpRec, "RSUN_OBS");
498 drms_copykey(sharpRec, mharpRec, "DSUN_OBS");
499 drms_copykey(sharpRec, mharpRec, "OBS_VR");
500 drms_copykey(sharpRec, mharpRec, "OBS_VW");
501 drms_copykey(sharpRec, mharpRec, "OBS_VN");
502 drms_copykey(sharpRec, mharpRec, "CRLN_OBS");
503 drms_copykey(sharpRec, mharpRec, "CRLT_OBS");
504 drms_copykey(sharpRec, mharpRec, "CAR_ROT");
505 drms_copykey(sharpRec, mharpRec, "SIZE_SPT");
506 drms_copykey(sharpRec, mharpRec, "AREA_SPT");
507 mbobra 1.1 drms_copykey(sharpRec, mharpRec, "DATE__OBS");
508 drms_copykey(sharpRec, mharpRec, "T_OBS");
509 drms_copykey(sharpRec, mharpRec, "T_MAXPIX");
510 drms_copykey(sharpRec, mharpRec, "QUALITY");
511 drms_copykey(sharpRec, mharpRec, "NPIX_SPT");
512 drms_copykey(sharpRec, mharpRec, "ARS_NCLN");
513 drms_copykey(sharpRec, mharpRec, "ARS_MODL");
514 drms_copykey(sharpRec, mharpRec, "ARS_EDGE");
515 drms_copykey(sharpRec, mharpRec, "ARS_BETA");
516 drms_copykey(sharpRec, mharpRec, "T_MID1");
517 drms_copykey(sharpRec, mharpRec, "T_CMPASS");
518
519 DRMS_Link_t *mHarpLink = hcon_lookup_lower(&sharpRec->links, "MTARP");
520 if (mHarpLink) {
521 drms_link_set("MTARP", sharpRec, mharpRec);
522 }
523
524 // set other keywords
525 setKeys(sharpRec, mharpRec, &mInfo);
526
527 // set space weather keywords
528 mbobra 1.1 computeSWIndex(swKeys_ptr, sharpRec, &mInfo);
529 printf("Space weather indices done.\n");
|
579 mbobra 1.1
580 float *inData;
581 int xsz = inArray->axis[0], ysz = inArray->axis[1];
582 if ((xsz != FOURK) || (ysz != FOURK)) { // for bitmap, make tmp full disk
583 float *inData0 = (float *) inArray->data;
584 inData = (float *) (calloc(FOURK2, sizeof(float)));
585 int x0 = (int) drms_getkey_float(harpRec, "CRPIX1", &status) - 1;
586 int y0 = (int) drms_getkey_float(harpRec, "CRPIX2", &status) - 1;
587 int ind_map;
588 for (int row = 0; row < ysz; row++) {
589 for (int col = 0; col < xsz; col++) {
590 ind_map = (row + y0) * FOURK + (col + x0);
591 inData[ind_map] = inData0[row * xsz + col];
592 }
593 }
594 drms_free_array(inArray); inArray = NULL;
595 } else {
596 inData = (float *) inArray->data;
597 }
598
599 // Mapping
600 mbobra 1.1
601 float *map = (float *) (malloc(nxny * sizeof(float)));
602 if (performSampling(map, inData, mInfo, interpOpt)) // Add interpOpt for different types, Aug 12 XS
603 {if (inArray) drms_free_array(inArray); free(map); return 1;}
604
605 // Write out
606
607 DRMS_Segment_t *outSeg = NULL;
608 outSeg = drms_segment_lookup(sharpRec, segName);
609 if (!outSeg) return 1;
610
611 // DRMS_Type_t arrayType = outSeg->info->type;
612 DRMS_Array_t *outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, map, &status);
613 if (status) {if (inArray) drms_free_array(inArray); free(map); return 1;}
614
615 // convert to needed data type
616
617 // drms_array_convert_inplace(outSeg->info->type, 0, 1, outArray); // Jan 02 2013
618
619 outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1];
620 // outArray->parent_segment = outSeg;
621 mbobra 1.1 outArray->israw = 0; // always compressed
622 outArray->bzero = outSeg->bzero;
623 outArray->bscale = outSeg->bscale;
624
625 status = drms_segment_write(outSeg, outArray, 0);
626 if (status) return 0;
627
628 if (inArray) drms_free_array(inArray);
629 if ((xsz != FOURK) || (ysz != FOURK)) free(inData); // Dec 18 2012
630 if (outArray) drms_free_array(outArray);
631 return 0;
632
633 }
634
635 /*
636 * Determine reference point coordinate and patch size according to keywords
637 * xc, yc are the coordinate of patch center, in degrees
638 * ncol and nrow are the final size
639 */
640
641 int findPosition(DRMS_Record_t *inRec, struct mapInfo *mInfo)
642 mbobra 1.1 {
643
644 int status = 0;
645 int harpnum = drms_getkey_int(inRec, "TARPNUM", &status);
646 TIME trec = drms_getkey_time(inRec, "T_REC", &status);
647 float disk_lonc = drms_getkey_float(inRec, "CRLN_OBS", &status);
648
649 /* Center coord */
650 // Changed into double Jun 16 2014 XS
651
652 double minlon = drms_getkey_double(inRec, "LONDTMIN", &status); if (status) return 1; // Stonyhurst lon
653 double maxlon = drms_getkey_double(inRec, "LONDTMAX", &status); if (status) return 1;
654 double minlat = drms_getkey_double(inRec, "LATDTMIN", &status); if (status) return 1;
655 double maxlat = drms_getkey_double(inRec, "LATDTMAX", &status); if (status) return 1;
656
657 // A bug fixer for HARP (per M. Turmon)
658 // When AR is below threshold, "LONDTMIN", "LONDTMAX" will be wrong
659 // Also keywords such as "SIZE" will be NaN
660 // We compute minlon & minlat then by
661 // LONDTMIN(t) = LONDTMIN(t0) + (t - t0) * OMEGA_DT
662
663 mbobra 1.1 // float psize = drms_getkey_float(inRec, "SIZE", &status);
664 // if (psize != psize) {
665
666 if (minlon != minlon || maxlon != maxlon) { // check lons instead of SIZE
667 TIME t0 = drms_getkey_time(inRec, "T_FRST1", &status); if (status) return 1; // changed from T_FRST to T_FRST1, T_FRST may not exist
668 double omega = drms_getkey_double(inRec, "OMEGA_DT", &status); if (status) return 1;
669 char firstRecQuery[100], t0_str[100];
670 sprint_time(t0_str, t0, "TAI", 0);
671 snprintf(firstRecQuery, 100, "%s[%d][%s]", inRec->seriesinfo->seriesname, harpnum, t0_str);
672 DRMS_RecordSet_t *tmpRS = drms_open_records(drms_env, firstRecQuery, &status);
673 if (status || tmpRS->n != 1) return 1;
674 DRMS_Record_t *tmpRec = tmpRS->records[0];
675 double minlon0 = drms_getkey_double(tmpRec, "LONDTMIN", &status); if (status) return 1;
676 double maxlon0 = drms_getkey_double(tmpRec, "LONDTMAX", &status); if (status) return 1;
677 minlon = minlon0 + (trec - t0) * omega / SECINDAY;
678 maxlon = maxlon0 + (trec - t0) * omega / SECINDAY;
679 printf("%s, %f, %f\n", firstRecQuery, minlon, maxlon);
680 }
681
682 mInfo->xc = (maxlon + minlon) / 2. + disk_lonc;
683 mInfo->yc = (maxlat + minlat) / 2.;
684 mbobra 1.1
685 /* Size */
686 // Rounded to 1.d3 precision first. Jun 16 2014 XS
687 // The previous fix does not work. LONDTMAX-LONDTMIN varies from frame to frame
688 // Need to find out the maximum possible difference, MAXLONDIFF (1.2e-4)
689 // Now, ncol = (maxlon-minlon)/xscale, if the decimal part is outside 0.5 \pm (MAXLONDIFF/xscale)
690 // proceed as it is. else, all use floor on ncol
691
692 float dpix = (MAXLONDIFF / mInfo->xscale) * 1.5; // "danger zone"
693 float ncol = (maxlon - minlon) / mInfo->xscale;
694 float d_ncol = fabs(ncol - floor(ncol) - 0.5); // distance to 0.5
695 if (d_ncol < dpix) {
696 mInfo->ncol = floor(ncol);
697 } else {
698 mInfo->ncol = round(ncol);
699 }
700
701 mInfo->nrow = round((maxlat - minlat) / mInfo->yscale);
702
703 return 0;
704
705 mbobra 1.1 }
706
707
708 /*
709 * Fetch ephemeris info from a DRMS record
710 */
711
712 int getEphemeris(DRMS_Record_t *inRec, struct ephemeris *ephem)
713 {
714
715 int status = 0;
716
717 float crota2 = drms_getkey_float(inRec, "CROTA2", &status); // rotation
718 double sina = sin(crota2 * RADSINDEG);
719 double cosa = cos(crota2 * RADSINDEG);
720
721 ephem->pa = - crota2 * RADSINDEG;
722 ephem->disk_latc = drms_getkey_float(inRec, "CRLT_OBS", &status) * RADSINDEG;
723 ephem->disk_lonc = drms_getkey_float(inRec, "CRLN_OBS", &status) * RADSINDEG;
724
725 float crvalx = 0.0;
726 mbobra 1.1 float crvaly = 0.0;
727 float crpix1 = drms_getkey_float(inRec, "IMCRPIX1", &status);
728 float crpix2 = drms_getkey_float(inRec, "IMCRPIX2", &status);
729 float cdelt = drms_getkey_float(inRec, "CDELT1", &status); // in arcsec, assumimg dx=dy
730 ephem->disk_xc = PIX_X(0.0,0.0) - 1.0; // Center of disk in pixel, starting at 0
731 ephem->disk_yc = PIX_Y(0.0,0.0) - 1.0;
732
733 float dSun = drms_getkey_float(inRec, "DSUN_OBS", &status);
734 float rSun_ref = drms_getkey_float(inRec, "RSUN_REF", &status);
735 if (status) rSun_ref = 6.96e8;
736
737 ephem->asd = asin(rSun_ref/dSun);
738 ephem->rSun = asin(rSun_ref / dSun) * RAD2ARCSEC / cdelt;
739
740 return 0;
741
742 }
743
744
745 /*
746 * Compute the coordinates to be sampled on full disk image
747 mbobra 1.1 * mInfo->xi_out & mInfo->zeta_out
748 * This is oversampled, its size is ncol0 & nrow0 as shown below
749 */
750
751 void findCoord(struct mapInfo *mInfo)
752 {
753
754 int ncol0 = mInfo->ncol * mInfo->nbin + (mInfo->nbin / 2) * 2; // pad with nbin/2 on edge to avoid NAN
755 int nrow0 = mInfo->nrow * mInfo->nbin + (mInfo->nbin / 2) * 2;
756
757 float xscale0 = mInfo->xscale / mInfo->nbin * RADSINDEG; // oversampling resolution
758 float yscale0 = mInfo->yscale / mInfo->nbin * RADSINDEG; // in rad
759
760 double lonc = mInfo->xc * RADSINDEG; // in rad
761 double latc = mInfo->yc * RADSINDEG;
762
763 double disk_lonc = (mInfo->ephem).disk_lonc;
764 double disk_latc = (mInfo->ephem).disk_latc;
765
766 double rSun = (mInfo->ephem).rSun;
767 double disk_xc = (mInfo->ephem).disk_xc / rSun;
768 mbobra 1.1 double disk_yc = (mInfo->ephem).disk_yc / rSun;
769 double pa = (mInfo->ephem).pa;
770
771 // Temp pointers
772
773 float *xi_out = mInfo->xi_out;
774 float *zeta_out = mInfo->zeta_out;
775
776 // start
777
778 double x, y; // map coord
779 double lat, lon; // helio coord
780 double xi, zeta; // image coord (for one point)
781
782 int ind_map;
783
784 for (int row0 = 0; row0 < nrow0; row0++) {
785 for (int col0 = 0; col0 < ncol0; col0++) {
786
787 ind_map = row0 * ncol0 + col0;
788
789 mbobra 1.1 x = (col0 + 0.5 - ncol0/2.) * xscale0; // in rad
790 y = (row0 + 0.5 - nrow0/2.) * yscale0;
791
792 /* map grid [x, y] corresponds to the point [lon, lat] in the heliographic coordinates.
793 * the [x, y] are in radians with respect of the center of the map [xcMap, ycMap].
794 * projection methods could be Mercator, Lambert, and many others. [maplonc, mapLatc]
795 * is the heliographic longitude and latitude of the map center. Both are in degree.
796 */
797
798 if (plane2sphere (x, y, latc, lonc, &lat, &lon, (int) mInfo->proj)) {
799 xi_out[ind_map] = -1;
800 zeta_out[ind_map] = -1;
801 continue;
802 }
803
804 /* map the grid [lon, lat] in the heliographic coordinates to [xi, zeta], a point in the
805 * image coordinates. The image properties, xCenter, yCenter, rSun, pa, ecc and chi are given.
806 */
807
808 if (sphere2img (lat, lon, disk_latc, disk_lonc, &xi, &zeta,
809 disk_xc, disk_yc, 1.0, pa, 0., 0., 0., 0.)) {
810 mbobra 1.1 xi_out[ind_map] = -1;
811 zeta_out[ind_map] = -1;
812 continue;
813 }
814
815 xi_out[ind_map] = xi * rSun;
816 zeta_out[ind_map] = zeta * rSun;
817
818 }
819 }
820
821 }
822
823
824 /*
825 * Sampling function
826 * oversampling by nbin, then binning using a Gaussian
827 * save results in outData, always of float type
828 */
829
830 int performSampling(float *outData, float *inData, struct mapInfo *mInfo, int interpOpt)
831 mbobra 1.1 {
832
833 int status = 0;
834 int ind_map;
835
836 int ncol0 = mInfo->ncol * mInfo->nbin + (mInfo->nbin / 2) * 2; // pad with nbin/2 on edge to avoid NAN
837 int nrow0 = mInfo->nrow * mInfo->nbin + (mInfo->nbin / 2) * 2;
838
839 // Changed Aug 12 2013, XS, for bitmaps
840 float *outData0;
841 if (interpOpt == 3 && mInfo->nbin == 1) {
842 outData0 = outData;
843 } else {
844 outData0 = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
845 }
846
847 float *xi_out = mInfo->xi_out;
848 float *zeta_out = mInfo->zeta_out;
849
850 // Interpolation
851
852 mbobra 1.1 struct fint_struct pars;
853 // Aug 12 2013, passed in as argument now
854
855 switch (interpOpt) {
856 case 0: // Wiener, 6 order, 1 constraint
857 init_finterpolate_wiener(&pars, 6, 1, 6, 2, 1, 1, NULL, dpath);
858 break;
859 case 1: // Cubic convolution
860 init_finterpolate_cubic_conv(&pars, 1., 3.);
861 break;
862 case 2: // Bilinear
863 init_finterpolate_linear(&pars, 1.);
864 break;
865 case 3: // Near neighbor
866 break;
867 default:
868 return 1;
869 }
870
871 //printf("interpOpt = %d, nbin = %d ", interpOpt, mInfo->nbin);
872 if (interpOpt == 3) { // Aug 6 2013, Xudong
873 mbobra 1.1 for (int row0 = 0; row0 < nrow0; row0++) {
874 for (int col0 = 0; col0 < ncol0; col0++) {
875 ind_map = row0 * ncol0 + col0;
876 outData0[ind_map] = nnb(inData, FOURK, FOURK, xi_out[ind_map], zeta_out[ind_map]);
877 }
878 }
879 } else {
880 finterpolate(&pars, inData, xi_out, zeta_out, outData0,
881 FOURK, FOURK, FOURK, ncol0, nrow0, ncol0, DRMS_MISSING_FLOAT);
882 }
883
884 // Rebinning, smoothing
885
886 if (interpOpt == 3 && mInfo->nbin == 1) {
887 return 0;
888 } else {
889 frebin(outData0, outData, ncol0, nrow0, mInfo->nbin, 1); // Gaussian
890 free(outData0);
891 }
892
893 return 0;
894 mbobra 1.1
895 }
896
897 /*
898 * Create Cutout record: top level subroutine
899 * Do the loops on segments and set the keywords here
900 * Work is done in writeCutout routine below
901 */
902
903 int createCutRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *contRec, DRMS_Record_t *sharpRec, struct swIndex *swKeys_ptr)
904 {
905
906 int status = 0;
907 int val;
908 int iHarpSeg;
909 int nMharpSegs = ARRLENGTH(MharpSegs);
910
911 // Cutout Mharp
912
913 for (iHarpSeg = 0; iHarpSeg < nMharpSegs; iHarpSeg++) {
914 if (writeCutout(sharpRec, mharpRec, mharpRec, MharpSegs[iHarpSeg])) {
915 mbobra 1.1 printf("Mharp cutout fails for %s\n", MharpSegs[iHarpSeg]);
916 printf("iHarpSeg nMharpSegs %d %d \n",iHarpSeg,nMharpSegs);
917 break;
918 }
919 }
920 if (iHarpSeg != nMharpSegs) {
921 SHOW("Cutout: segment number mismatch\n");
922 return 1; // if failed
923 }
924 printf("Magnetogram cutout done.\n");
925
926 // Cutout Continuum
927
928 if (writeCutout(sharpRec, contRec, mharpRec, "continuum")) {
929 printf("Continuum cutout failed\n");
930 return 1;
931 }
932 printf("Intensitygram cutout done.\n");
933
934 // Keywords & Links
935 copy_patch_keys(mharpRec, sharpRec);
936 mbobra 1.1 copy_geo_keys(mharpRec, sharpRec);
937 // rename HARPNUM to TARPNUM
938 val = drms_getkey_double(mharpRec, "HARPNUM", &status);
939 drms_setkey_double(sharpRec, "TARPNUM", val);
940 // copy everything else
941 drms_copykey(sharpRec, mharpRec, "T_REC");
942 drms_copykey(sharpRec, mharpRec, "CDELT1");
943 drms_copykey(sharpRec, mharpRec, "RSUN_OBS");
944 drms_copykey(sharpRec, mharpRec, "DSUN_OBS");
945 drms_copykey(sharpRec, mharpRec, "OBS_VR");
946 drms_copykey(sharpRec, mharpRec, "OBS_VW");
947 drms_copykey(sharpRec, mharpRec, "OBS_VN");
948 drms_copykey(sharpRec, mharpRec, "CRLN_OBS");
949 drms_copykey(sharpRec, mharpRec, "CRLT_OBS");
950 drms_copykey(sharpRec, mharpRec, "CAR_ROT");
951 drms_copykey(sharpRec, mharpRec, "SIZE_SPT");
952 drms_copykey(sharpRec, mharpRec, "AREA_SPT");
953 drms_copykey(sharpRec, mharpRec, "DATE__OBS");
954 drms_copykey(sharpRec, mharpRec, "T_OBS");
955 drms_copykey(sharpRec, mharpRec, "T_MAXPIX");
956 drms_copykey(sharpRec, mharpRec, "QUALITY");
957 mbobra 1.1 drms_copykey(sharpRec, mharpRec, "NPIX_SPT");
958 drms_copykey(sharpRec, mharpRec, "ARS_NCLN");
959 drms_copykey(sharpRec, mharpRec, "ARS_MODL");
960 drms_copykey(sharpRec, mharpRec, "ARS_EDGE");
961 drms_copykey(sharpRec, mharpRec, "ARS_BETA");
962 drms_copykey(sharpRec, mharpRec, "T_MID1");
963 drms_copykey(sharpRec, mharpRec, "T_CMPASS");
964
965 DRMS_Link_t *mHarpLink = hcon_lookup_lower(&sharpRec->links, "MTARP");
966 if (mHarpLink) {
967 drms_link_set("MTARP", sharpRec, mharpRec);
968 }
969
970 setSWIndex(sharpRec, swKeys_ptr); // Set space weather indices
971 setKeys(sharpRec, mharpRec, NULL); // Set all other keywords, NULL specifies cutout
972
973 // Stats
974
975 int nCutSegs = 3;
976 for (int iSeg = 0; iSeg < 3; iSeg++) {
977 DRMS_Segment_t *outSeg = drms_segment_lookupnum(sharpRec, iSeg);
978 mbobra 1.1 DRMS_Array_t *outArray = drms_segment_read(outSeg, DRMS_TYPE_FLOAT, &status);
979 set_statistics(outSeg, outArray, 1);
980 drms_free_array(outArray);
981 }
982
983 return 0;
984
985 }
986
987 /*
988 * Get cutout and write segment
989 */
990
991 int writeCutout(DRMS_Record_t *outRec, DRMS_Record_t *inRec, DRMS_Record_t *harpRec, char *SegName)
992 {
993
994 int status = 0;
995
996 DRMS_Segment_t *inSeg = NULL, *outSeg = NULL;
997 DRMS_Array_t *cutoutArray = NULL;
998 // DRMS_Type_t arrayType;
999 mbobra 1.1
1000 int ll[2], ur[2], nx, ny, nxny; // lower-left and upper right coords
1001
1002 /* Info */
1003
1004 inSeg = drms_segment_lookup(inRec, SegName);
1005 if (!inSeg) return 1;
1006 //printf("SegName=%s\n",SegName); fflush(stdout);
1007 nx = (int) drms_getkey_float(harpRec, "CRSIZE1", &status);
1008 ny = (int) drms_getkey_float(harpRec, "CRSIZE2", &status);
1009 nxny = nx * ny;
1010 ll[0] = (int) drms_getkey_float(harpRec, "CRPIX1", &status) - 1; if (status) return 1;
1011 ll[1] = (int) drms_getkey_float(harpRec, "CRPIX2", &status) - 1; if (status) return 1;
1012 ur[0] = ll[0] + nx - 1; if (status) return 1;
1013 ur[1] = ll[1] + ny - 1; if (status) return 1;
1014 if (inSeg->axis[0] == nx && inSeg->axis[1] == ny) { // for bitmaps, infomaps, etc.
1015 cutoutArray = drms_segment_read(inSeg, DRMS_TYPE_DOUBLE, &status);
1016 if (status) return 1;
1017 } else if (inSeg->axis[0] == FOURK && inSeg->axis[1] == FOURK) { // for full disk ones
1018 cutoutArray = drms_segment_readslice(inSeg, DRMS_TYPE_DOUBLE, ll, ur, &status);
1019 if (status) return 1;
1020 mbobra 1.1 } else {
1021 return 1;
1022 }
1023 /* Write out */
1024 outSeg = drms_segment_lookup(outRec, SegName);
1025 if (!outSeg) return 1;
1026 outSeg->axis[0] = cutoutArray->axis[0];
1027 outSeg->axis[1] = cutoutArray->axis[1];
1028 cutoutArray->israw = 0; // always compressed
1029 cutoutArray->bzero = outSeg->bzero;
1030 cutoutArray->bscale = outSeg->bscale; // Same as inArray's
1031 status = drms_segment_write(outSeg, cutoutArray, 0);
1032 drms_free_array(cutoutArray);
1033 if (status) return 1;
1034 //printf("line1068\n"); fflush(stdout);
1035 return 0;
1036
1037 }
1038
1039
1040 /*
1041 mbobra 1.1 * Compute space weather indices
1042 */
1043
1044 void computeSWIndex(struct swIndex *swKeys_ptr, DRMS_Record_t *inRec, struct mapInfo *mInfo)
1045 {
1046
1047 int status = 0;
1048 int nx = mInfo->ncol, ny = mInfo->nrow;
1049 int nxny = nx * ny;
1050 int dims[2] = {nx, ny};
1051
|
1144 mbobra 1.1 };
1145
1146 /*
1147 * Set all keywords, no error checking for now
1148 */
1149
1150 void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *mharpRec, struct mapInfo *mInfo)
1151 {
1152
1153 int status = 0;
1154
1155 // Change a few geometry keywords for CEA & cutout records
1156 if (mInfo != NULL) { // CEA
1157 printf("Calculating CEA keys\n");
1158 drms_setkey_float(outRec, "CRPIX1", mInfo->ncol/2. + 0.5);
1159 drms_setkey_float(outRec, "CRPIX2", mInfo->nrow/2. + 0.5);
1160 drms_setkey_float(outRec, "CRVAL1", mInfo->xc);
1161 drms_setkey_float(outRec, "CRVAL2", mInfo->yc);
1162 drms_setkey_float(outRec, "CDELT1", mInfo->xscale);
1163 drms_setkey_float(outRec, "CDELT2", mInfo->yscale);
1164 drms_setkey_string(outRec, "CUNIT1", "degree");
1165 mbobra 1.1 drms_setkey_string(outRec, "CUNIT2", "degree");
1166 char key[64];
1167 snprintf (key, 64, "CRLN-%s", wcsCode[(int) mInfo->proj]);
1168 drms_setkey_string(outRec, "CTYPE1", key);
1169 snprintf (key, 64, "CRLT-%s", wcsCode[(int) mInfo->proj]);
1170 drms_setkey_string(outRec, "CTYPE2", key);
1171 drms_setkey_float(outRec, "CROTA2", 0.0);
1172 // Set BUNIT for each segment
1173 int nSeg = 3;
1174 for (int iSeg = 0; iSeg < nSeg; iSeg++) {
1175 DRMS_Segment_t *outSeg = NULL;
1176 outSeg = drms_segment_lookup(outRec, CEASegs[iSeg]);
1177 if (!outSeg) continue;
1178 char bunit_xxx[20];
1179 sprintf(bunit_xxx, "BUNIT_%03d", iSeg);
1180 //printf("%s, %s\n", bunit_xxx, CEABunits[iSeg]);
1181 drms_setkey_string(outRec, bunit_xxx, CEABunits[iSeg]);
1182 }
1183
1184 } else { // Cutout
1185
1186 mbobra 1.1 float disk_xc, disk_yc;
1187 disk_xc = drms_getkey_float(mharpRec, "IMCRPIX1", &status);
1188 disk_yc = drms_getkey_float(mharpRec, "IMCRPIX2", &status);
1189 float x_ll = drms_getkey_float(mharpRec, "CRPIX1", &status);
1190 float y_ll = drms_getkey_float(mharpRec, "CRPIX2", &status);
1191 // Defined as disk center's pixel address wrt lower-left of cutout
1192 drms_setkey_float(outRec, "CRPIX1", disk_xc - x_ll + 1.);
1193 drms_setkey_float(outRec, "CRPIX2", disk_yc - y_ll + 1.);
1194 // Always 0.
1195 drms_setkey_float(outRec, "CRVAL1", 0);
1196 drms_setkey_float(outRec, "CRVAL2", 0);
1197
1198 // Jan 2 2014 XS
1199 int nSeg = ARRLENGTH(CutSegs);
1200 for (int iSeg = 0; iSeg < nSeg; iSeg++)
1201 {
1202 DRMS_Segment_t *outSeg = NULL;
1203 outSeg = drms_segment_lookup(outRec, CutSegs[iSeg]);
1204 if (!outSeg) continue;
1205 // Set Bunit
1206 char bunit_xxx[20];
1207 mbobra 1.1 sprintf(bunit_xxx, "BUNIT_%03d", iSeg);
1208 //printf("%s, %s\n", bunit_xxx, CutBunits[iSeg]);
1209 drms_setkey_string(outRec, bunit_xxx, CutBunits[iSeg]);
1210 }
1211 }
1212
1213
1214 TIME val, trec, tnow, UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */
1215 tnow = (double)time(NULL);
1216 tnow += UNIX_epoch;
1217
1218 val = drms_getkey_time(mharpRec, "DATE", &status);
1219 drms_setkey_time(outRec, "DATE", tnow);
1220
1221 // set cvs commit version into keyword CODEVER7
1222 char *cvsinfo = strdup("$Id");
1223 char *cvsinfo2 = smarp_functions_version();
1224 char cvsinfoall[2048];
1225 strcat(cvsinfoall,cvsinfo);
1226 strcat(cvsinfoall,"\n");
1227 strcat(cvsinfoall,cvsinfo2);
1228 mbobra 1.1 status = drms_setkey_string(outRec, "CODEVER7", cvsinfoall);
1229
1230 };
1231
1232 /* ############# Nearest neighbour interpolation ############### */
1233
1234 float nnb (float *f, int nx, int ny, double x, double y)
1235 {
1236
1237 if (x <= -0.5 || y <= -0.5 || x > nx - 0.5 || y > ny - 0.5)
1238 return DRMS_MISSING_FLOAT;
1239 int ilow = floor (x);
1240 int jlow = floor (y);
1241 int i = ((x - ilow) > 0.5) ? ilow + 1 : ilow;
1242 int j = ((y - jlow) > 0.5) ? jlow + 1 : jlow;
1243 return f[j * nx + i];
1244
1245 }
1246
1247 /* ################## Wrapper for Jesper's rebin code ################## */
1248
1249 mbobra 1.1 void frebin (float *image_in, float *image_out, int nx, int ny, int nbin, int gauss)
1250 {
1251
1252 struct fresize_struct fresizes;
1253 int nxout, nyout, xoff, yoff;
1254 int nlead = nx;
1255
1256 nxout = nx / nbin; nyout = ny / nbin;
1257 if (gauss && nbin != 1)
1258 init_fresize_gaussian(&fresizes, (nbin / 2), (nbin / 2 * 2), nbin); // for nbin=3, sigma=1, half truncate width=2
1259 else
1260 init_fresize_bin(&fresizes, nbin);
1261 xoff = nbin / 2 + nbin / 2;
1262 yoff = nbin / 2 + nbin / 2;
1263 fresize(&fresizes, image_in, image_out, nx, ny, nlead, nxout, nyout, nxout, xoff, yoff, DRMS_MISSING_FLOAT);
1264
1265 }
|
Return to
smarp.c
CVS log
Up to [Development] / JSOC / proj / sharp / apps