JSOC/proj/sharp/apps/sharp.c - diff

Return to sharp.c CVS log

Up to [Development] / JSOC / proj / sharp / apps

Diff for /JSOC/proj/sharp/apps/sharp.c between version 1.14 and 1.18

version 1.14, 2013/05/21 17:06:32

version 1.18, 2013/11/02 20:31:08

Line 22

* v0.2 Sep 04 2012

* v0.3 Dec 18 2012

* v0.4 Jan 02 2013

* v0.5 Jan 23 2013

* v0.6 Aug 12 2013

* Notes:

* v0.0

Line 39

Line 41

* Fixed memory leakage of 0.15G per rec; denoted with "Dec 18"

* v0.4

* Took out convert_inplace(). Was causing all the images to be int

* v0.5

* Corrected ephemeris keywords, added argument mInfo for setKeys()

* v0.6

* Changes in remapping of bitmap and conf_disambig, now near neighbor without anti-aliasing

* Example:

* sharp "mharp=hmi.Mharp_720s[1404][2012.02.20_10:00]" \

Line 90

Line 96

// Nyqvist rate at disk center is 0.03 degree. Oversample above 0.015 degree

#define NYQVIST (0.015)

// Some other things

#ifndef MIN

#define MIN(a,b) (((a)<(b)) ? (a) : (b))

#endif

Line 117

Line 124

#define INTERP 0

#define dpath "/home/jsoc/cvs/Development/JSOC"

/* ========================================================================================================== */

// Space weather keywords

Line 173 enum projection {

Line 181 enum projection {

lambert

};

// WSC code

char *wcsCode[] = {"CAR", "CAS", "MER", "CEA", "GLS", "TAN", "ARC", "STG",

"SIN", "ZEA"};

// Ephemeris information

struct ephemeris {

double disk_lonc, disk_latc;

Line 191 struct mapInfo {

Line 203 struct mapInfo {

float *xi_out, *zeta_out; // coordinate on full disk image to sample at

};

/* ========================================================================================================== */

/* Get all input data series */

Line 232 int getEphemeris(DRMS_Record_t *inRec, s

Line 243 int getEphemeris(DRMS_Record_t *inRec, s

void findCoord(struct mapInfo *mInfo);

/* Mapping function */

int performSampling(float *outData, float *inData, struct mapInfo *mInfo);

int performSampling(float *outData, float *inData, struct mapInfo *mInfo, int interpOpt);

/* Performing local vector transformation */

void vectorTransform(float *bx_map, float *by_map, float *bz_map, struct mapInfo *mInfo);

Line 269 void computeSWIndex(struct swIndex *swKe

Line 280 void computeSWIndex(struct swIndex *swKe

void setSWIndex(DRMS_Record_t *outRec, struct swIndex *swKeys_ptr);

/* Set all keywords, no error checking for now */

void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *inRec);

void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *inRec, struct mapInfo *mInfo);

// ===================

Line 311 char *CutSegs[] = {"magnetogram", "bitma

Line 322 char *CutSegs[] = {"magnetogram", "bitma

"disambig", "conf_disambig"};

char *CEASegs[] = {"magnetogram", "bitmap", "Dopplergram", "continuum", "disambig",

BR_SEG_CEA, BT_SEG_CEA, BP_SEG_CEA, BR_ERR_SEG_CEA, BT_ERR_SEG_CEA, BP_ERR_SEG_CEA};

/* ========================================================================================================== */

char *module_name = "sharp";

char *version_id = "2012 Dec 18"; /* Version number */

int seed;

ModuleArgs_t module_args[] =

Line 352 int DoIt(void)

Line 363 int DoIt(void)

contQuery = (char *) params_get_str(&cmdparams, "cont");

sharpCeaQuery = (char *) params_get_str(&cmdparams, "sharp_cea");

sharpCutQuery = (char *) params_get_str(&cmdparams, "sharp_cut");

seed = params_get_int(&cmdparams, "seed");

Line 592 int createCeaRecord(DRMS_Record_t *mharp

Line 602 int createCeaRecord(DRMS_Record_t *mharp

mInfo.proj = (enum projection) cyleqa; // projection method

mInfo.xscale = XSCALE;

mInfo.yscale = YSCALE;

mInfo.nbin = NBIN;

int ncol0, nrow0; // oversampled map size

// Get ephemeris

Line 608 int createCeaRecord(DRMS_Record_t *mharp

Line 619 int createCeaRecord(DRMS_Record_t *mharp

return 1;

}

// ========================================

// Do this for all bitmaps, Aug 12 2013 XS

// ========================================

mInfo.nbin = 1; // for bitmaps. suppress anti-aliasing

ncol0 = mInfo.ncol;

nrow0 = mInfo.nrow;

mInfo.xi_out = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));

mInfo.zeta_out = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));

findCoord(&mInfo); // compute it here so it could be shared by the following 4 functions

if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "bitmap")) {

SHOW("CEA: mapping bitmap error\n");

return 1;

}

printf("Bitmap mapping done.\n");

if (mapScaler(sharpRec, bharpRec, mharpRec, &mInfo, "conf_disambig")) {

SHOW("CEA: mapping conf_disambig error\n");

return 1;

}

printf("Conf disambig mapping done.\n");

free(mInfo.xi_out);

free(mInfo.zeta_out);

// ========================================

// Do this again for floats, Aug 12 2013 XS

// ========================================

// Create xi_out, zeta_out array in mInfo:

// Coordinates to sample in original full disk image

int ncol0, nrow0; // oversampled map size

mInfo.nbin = NBIN;

ncol0 = mInfo.ncol * mInfo.nbin + (mInfo.nbin / 2) * 2; // pad with nbin/2 on edge to avoid NAN

nrow0 = mInfo.nrow * mInfo.nbin + (mInfo.nbin / 2) * 2;

Line 626 int createCeaRecord(DRMS_Record_t *mharp

Line 668 int createCeaRecord(DRMS_Record_t *mharp

SHOW("CEA: mapping magnetogram error\n");

return 1;

}

if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "bitmap")) {

SHOW("CEA: mapping bitmap error\n");

return 1;

}

printf("Magnetogram mapping done.\n");

if (mapScaler(sharpRec, dopRec, mharpRec, &mInfo, "Dopplergram")) {

SHOW("CEA: mapping dopplergram error\n");

return 1;

Line 643 int createCeaRecord(DRMS_Record_t *mharp

Line 682 int createCeaRecord(DRMS_Record_t *mharp

}

printf("Intensitygram mapping done.\n");

if (mapScaler(sharpRec, bharpRec, mharpRec, &mInfo, "conf_disambig")) {

SHOW("CEA: mapping conf_disambig error\n");

return 1;

}

printf("Conf disambig mapping done.\n");

// Mapping vector B

if (mapVectorB(sharpRec, bharpRec, &mInfo)) {

Line 675 int createCeaRecord(DRMS_Record_t *mharp

Line 708 int createCeaRecord(DRMS_Record_t *mharp

DRMS_Link_t *bHarpLink = hcon_lookup_lower(&sharpRec->links, "BHARP");

if (bHarpLink) drms_link_set("BHARP", sharpRec, bharpRec);

setKeys(sharpRec, bharpRec); // Set all other keywords

setKeys(sharpRec, bharpRec, &mInfo); // Set all other keywords

drms_copykey(sharpRec, mharpRec, "QUALITY"); // copied from los records

// Space weather

Line 717 int mapScaler(DRMS_Record_t *sharpRec, D

Line 750 int mapScaler(DRMS_Record_t *sharpRec, D

int status = 0;

int nx = mInfo->ncol, ny = mInfo->nrow, nxny = nx * ny;

int dims[2] = {nx, ny};

int interpOpt = INTERP; // Aug 12 XS, default, overridden below for bitmaps and conf_disambig

// Input full disk array

Line 731 int mapScaler(DRMS_Record_t *sharpRec, D

Line 765 int mapScaler(DRMS_Record_t *sharpRec, D

float *inData;

int xsz = inArray->axis[0], ysz = inArray->axis[1];

if ((xsz != FOURK) || (ysz != FOURK)) { // for bitmap, make tmp full disk

interpOpt = 3; // Aug 12 XS, near neighbor

float *inData0 = (float *) inArray->data;

inData = (float *) (calloc(FOURK2, sizeof(float)));

int x0 = (int) drms_getkey_float(harpRec, "CRPIX1", &status) - 1;

Line 750 int mapScaler(DRMS_Record_t *sharpRec, D

Line 785 int mapScaler(DRMS_Record_t *sharpRec, D

// Mapping

float *map = (float *) (malloc(nxny * sizeof(float)));

if (performSampling(map, inData, mInfo))

if (performSampling(map, inData, mInfo, interpOpt)) // Add interpOpt for different types, Aug 12 XS

{if (inArray) drms_free_array(inArray); free(map); return 1;}

// Write out

Line 813 int mapVectorB(DRMS_Record_t *sharpRec,

Line 848 int mapVectorB(DRMS_Record_t *sharpRec,

float *bx_map = NULL, *by_map = NULL, *bz_map = NULL; // intermediate maps, in CCD bxyz representation

bx_map = (float *) (malloc(nxny * sizeof(float)));

if (performSampling(bx_map, bx_img, mInfo))

if (performSampling(bx_map, bx_img, mInfo, INTERP))

{free(bx_img); free(by_img); free(bz_img); free(bx_map); return 1;}

by_map = (float *) (malloc(nxny * sizeof(float)));

if (performSampling(by_map, by_img, mInfo))

if (performSampling(by_map, by_img, mInfo, INTERP))

{free(bx_img); free(by_img); free(bz_img); free(bz_map); return 1;}

bz_map = (float *) (malloc(nxny * sizeof(float)));

if (performSampling(bz_map, bz_img, mInfo))

if (performSampling(bz_map, bz_img, mInfo, INTERP))

{free(bx_img); free(by_img); free(bz_img); free(bz_map); return 1;}

free(bx_img); free(by_img); free(bz_img);

Line 997 int getEphemeris(DRMS_Record_t *inRec, s

Line 1032 int getEphemeris(DRMS_Record_t *inRec, s

float crpix1 = drms_getkey_float(inRec, "IMCRPIX1", &status);

float crpix2 = drms_getkey_float(inRec, "IMCRPIX2", &status);

float cdelt = drms_getkey_float(inRec, "CDELT1", &status); // in arcsec, assumimg dx=dy

printf("cdelt=%f\n",cdelt);

ephem->disk_xc = PIX_X(0.0,0.0) - 1.0; // Center of disk in pixel, starting at 0

ephem->disk_yc = PIX_Y(0.0,0.0) - 1.0;

Line 1101 void findCoord(struct mapInfo *mInfo)

Line 1135 void findCoord(struct mapInfo *mInfo)

int performSampling(float *outData, float *inData, struct mapInfo *mInfo)

int performSampling(float *outData, float *inData, struct mapInfo *mInfo, int interpOpt)

{

int status = 0;

int ind_map;

int ncol0 = mInfo->ncol * mInfo->nbin + (mInfo->nbin / 2) * 2; // pad with nbin/2 on edge to avoid NAN

int nrow0 = mInfo->nrow * mInfo->nbin + (mInfo->nbin / 2) * 2;

float *outData0 = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));

// Changed Aug 12 2013, XS, for bitmaps

float *outData0;

if (interpOpt == 3 && mInfo->nbin == 1) {

outData0 = outData;

} else {

outData0 = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));

}

float *xi_out = mInfo->xi_out;

float *zeta_out = mInfo->zeta_out;

Line 1117 int performSampling(float *outData, floa

Line 1158 int performSampling(float *outData, floa

// Interpolation

struct fint_struct pars;

int interpOpt = INTERP; // Use Wiener by default, 6 order, 1 constraint

// Aug 12 2013, passed in as argument now

switch (interpOpt) {

case 0: // Wiener, 6 order, 1 constraint

Line 1129 int performSampling(float *outData, floa

Line 1170 int performSampling(float *outData, floa

case 2: // Bilinear

init_finterpolate_linear(&pars, 1.);

break;

case 3: // Near neighbor

break;

default:

return 1;

}

printf("interpOpt = %d, nbin = %d ", interpOpt, mInfo->nbin);

if (interpOpt == 3) { // Aug 6 2013, Xudong

for (int row0 = 0; row0 < nrow0; row0++) {

for (int col0 = 0; col0 < ncol0; col0++) {

ind_map = row0 * ncol0 + col0;

outData0[ind_map] = nnb(inData, FOURK, FOURK, xi_out[ind_map], zeta_out[ind_map]);

}

} else {

finterpolate(&pars, inData, xi_out, zeta_out, outData0,

FOURK, FOURK, FOURK, ncol0, nrow0, ncol0, DRMS_MISSING_FLOAT);

}

// Rebinning, smoothing

if (interpOpt == 3 && mInfo->nbin == 1) {

return 0;

} else {

frebin(outData0, outData, ncol0, nrow0, mInfo->nbin, 1); // Gaussian

free(outData0); // Dec 18 2012

}

Line 1561 int createCutRecord(DRMS_Record_t *mharp

Line 1618 int createCutRecord(DRMS_Record_t *mharp

if (bHarpLink) drms_link_set("BHARP", sharpRec, bharpRec);

setSWIndex(sharpRec, swKeys_ptr); // Set space weather indices

setKeys(sharpRec, bharpRec); // Set all other keywords

setKeys(sharpRec, bharpRec, NULL); // Set all other keywords, NULL specifies cutout

// Stats

Line 1716 void computeSWIndex(struct swIndex *swKe

Line 1773 void computeSWIndex(struct swIndex *swKe

float *bx_err = (float *) bx_errArray->data; // bx_err

// Get emphemeris

float cdelt1_orig = drms_getkey_float(inRec, "CDELT1", &status);

//float cdelt1_orig = drms_getkey_float(inRec, "CDELT1", &status);

float cdelt1 = drms_getkey_float(inRec, "CDELT1", &status);

float dsun_obs = drms_getkey_float(inRec, "DSUN_OBS", &status);

double rsun_ref = drms_getkey_double(inRec, "RSUN_REF", &status);

double rsun_obs = drms_getkey_double(inRec, "RSUN_OBS", &status);

Line 1727 void computeSWIndex(struct swIndex *swKe

Line 1782 void computeSWIndex(struct swIndex *swKe

float crpix1 = drms_getkey_float(inRec, "CRPIX1", &status);

float crpix2 = drms_getkey_float(inRec, "CRPIX2", &status);

// Temp arrays

// convert cdelt1_orig from degrees to arcsec

float cdelt1 = (atan((rsun_ref*cdelt1_orig*RADSINDEG)/(dsun_obs)))*(1/RADSINDEG)*(3600.);

// Temp arrays

float *bh = (float *) (malloc(nxny * sizeof(float)));

float *bt = (float *) (malloc(nxny * sizeof(float)));

float *jz = (float *) (malloc(nxny * sizeof(float)));

Line 1841 void computeSWIndex(struct swIndex *swKe

Line 1898 void computeSWIndex(struct swIndex *swKe

swKeys_ptr->totpot_err = DRMS_MISSING_FLOAT;

}

if (computeShearAngle(bx_err, by_err, bh_err, bx, by, bz, bpx, bpy, bpz, dims,

if (computeShearAngle(bx_err, by_err, bz_err, bx, by, bz, bpx, bpy, bpz, dims,

&(swKeys_ptr->meanshear_angle), &(swKeys_ptr->meanshear_angle_err), &(swKeys_ptr->area_w_shear_gt_45),

mask, bitmask)) {

swKeys_ptr->meanshear_angle = DRMS_MISSING_FLOAT; // If fail, fill in NaN

Line 1914 void setSWIndex(DRMS_Record_t *outRec, s

Line 1971 void setSWIndex(DRMS_Record_t *outRec, s

void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *inRec)

void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *inRec, struct mapInfo *mInfo)

{

copy_me_keys(inRec, outRec);

copy_patch_keys(inRec, outRec);

copy_geo_keys(inRec, outRec);

copy_ambig_keys(inRec, outRec);

int status = 0;

// Change a few geometry keywords for CEA records

if (mInfo != NULL) {

drms_setkey_float(outRec, "CRPIX1", mInfo->ncol/2. + 0.5);

drms_setkey_float(outRec, "CRPIX2", mInfo->nrow/2. + 0.5);

drms_setkey_float(outRec, "CRVAL1", mInfo->xc);

drms_setkey_float(outRec, "CRVAL2", mInfo->yc);

drms_setkey_float(outRec, "CDELT1", mInfo->xscale);

drms_setkey_float(outRec, "CDELT2", mInfo->yscale);

drms_setkey_string(outRec, "CUNIT1", "degree");

drms_setkey_string(outRec, "CUNIT2", "degree");

char key[64];

snprintf (key, 64, "CRLN-%s", wcsCode[(int) mInfo->proj]);

drms_setkey_string(outRec, "CTYPE1", key);

snprintf (key, 64, "CRLT-%s", wcsCode[(int) mInfo->proj]);

drms_setkey_string(outRec, "CTYPE2", key);

drms_setkey_float(outRec, "CROTA2", 0.0);

} else {

float disk_xc = drms_getkey_float(inRec, "IMCRPIX1", &status);

float disk_yc = drms_getkey_float(inRec, "IMCRPIX2", &status);

float x_ll = drms_getkey_float(inRec, "CRPIX1", &status);

float y_ll = drms_getkey_float(inRec, "CRPIX2", &status);

// Defined as disk center's pixel address wrt lower-left of cutout

drms_setkey_float(outRec, "CRPIX1", disk_xc - x_ll + 1.);

drms_setkey_float(outRec, "CRPIX2", disk_yc - y_ll + 1.);

// Always 0.

drms_setkey_float(outRec, "CRVAL1", 0);

drms_setkey_float(outRec, "CRVAL2", 0);

}

char timebuf[1024];

float UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */

double val;

int status = DRMS_SUCCESS;

val = drms_getkey_double(inRec, "DATE",&status);

drms_setkey_double(outRec, "DATE_B", val);

Line 1933 void setKeys(DRMS_Record_t *outRec, DRMS

Line 2026 void setKeys(DRMS_Record_t *outRec, DRMS

// set cvs commit version into keyword HEADER

char *cvsinfo = strdup("$Id$");

// status = drms_setkey_string(outRec, "HEADER", cvsinfo);

char *cvsinfo2 = sw_functions_version();

status = drms_setkey_string(outRec, "CODEVER7", cvsinfo);

char cvsinfoall[2048];

strcat(cvsinfoall,cvsinfo);

strcat(cvsinfoall,"\n");

strcat(cvsinfoall,cvsinfo2);

status = drms_setkey_string(outRec, "CODEVER7", cvsinfoall);

};

Legend:

Removed from v.1.14
changed lines
	Added in v.1.18

Karen Tian