JSOC_Documentation: proj/mag/polcorr/apps/update

00001 /*
00002  * Module name:     update_poldb.c
00003  *
00004  * Description:
00005  *          Collect 360 deg of polar field data from Synoptic maps
00006  *          centered at Mar 06/Sep 08, save orig and smoothed version
00007  *          maybe overridden by T_CEN if not available
00008  *          Fit down to 60 latitude, using Stereographic and bilinear regardless
00009  *
00010  * Written by:      Xudong Sun (xudongs@stanford.edu)
00011  *                  Image fitting code by Rasmus Larsen (rmlarsen@gmail.com)
00012  *
00013  * External function:
00014  *          heliographic_coord.c    polar_remap.c   interp.c    fitimage.c
00015  *
00016  * Version:
00017  *          v1.0        Sep 10 2013
00018  *
00019  * Issues:
00020  *          v1.0
00021  *
00022  *
00023  * Example:
00024  *  update_poldb year=2012 ns=0 in=hmi.synoptic_mr_720s out=su_xudong.hmi_poldb scale=1 -o
00025  *
00026  */
00027 
00028 #include <jsoc_main.h>
00029 #include <stdio.h>
00030 #include <stdlib.h>
00031 #include <math.h>
00032 #include <sys/time.h>
00033 #include <sys/resource.h>
00034 
00035 // Heliographic coord, lib function
00036 #include "heliographic_coords.c"
00037 
00038 // Polar projection
00039 #include "polar_remap.c"
00040 
00041 // Image fitting
00042 #include "fitimage.c"
00043 
00044 /* ################## Preset ################## */
00045 
00046 // Image size
00047 #define IMGDIM  (600)
00048 // Lower limit of latitude for fitting
00049 #define LAT0    (60.)
00050 
00051 #define LAT1    (62.)
00052 #define LAT2    (75.)
00053 
00054 #define EPSL (1.e-5)
00055 #define SECSINDAY (86400.0)
00056 #define DTOR    (M_PI / 180.)
00057 
00058 #define S_DATE ".03.06_02:00:00_UTC"
00059 #define N_DATE ".09.08_02:00:00_UTC"
00060 
00061 #define EPOCH_DIFF ((double)220924792.000) /* 1970.01.01_00:00:00_UTC */
00062 
00063 #define INSEGNAME   "synopMr"
00064 
00065 /* ############# Macros ############# */
00066 
00067 #ifndef MIN
00068 #define MIN(a,b) (((a)<(b)) ? (a) : (b))
00069 #endif
00070 #ifndef MAX
00071 #define MAX(a,b) (((a)>(b)) ? (a) : (b))
00072 #endif
00073 
00074 #define ARRLENGTH(ARR) (sizeof(ARR) / sizeof(ARR[0]))
00075 #define SHOW(msg) {printf("%s", msg); fflush(stdout);}
00076 #define DIE(msg) {fflush(stdout); fprintf(stderr,"%s (status=%d)\n", msg, status); return(status);}
00077 
00078 #define FREE_ARR(arr) {if (arr) {free(arr);}}
00079 
00080 /* ############# Keywords ############# */
00081 
00082 char *synop_keys[] = {
00083     "CTYPE1", "CTYPE2", "CDELT1", "CDELT2", "CUNIT1", "CUNIT2", "WCSNAME"};
00084 
00085 /* ################## Aux function ################## */
00086 
00087 /* Get current year */
00088 int currentYear(void)
00089 {
00090     time_t timeval;
00091     struct tm *tp;
00092     time (&timeval);
00093     tp = gmtime(&timeval);
00094     return (tp->tm_year + 1900);
00095 }
00096 
00097 /* Convert TIME to current year */
00098 int time2Year(TIME t)
00099 {
00100     time_t timeval = (time_t)(t + EPOCH_DIFF);
00101     struct tm *tp = gmtime(&timeval);
00102     return (tp->tm_year + 1900);
00103 }
00104 
00105 /* ################## Checking functions ################## */
00106 
00107 /* Check parameters, set center time */
00108 int checkParams(int year, int ns, int nmax, float scale, TIME *t_cen_ptr);
00109 
00110 /* Check output availability */
00111 int checkOutput(char *outRSName, int year, int ns, int ovwt);
00112 
00113 /* Check input availability, determine crot, etc. */
00114 int checkInput(char *inRSName, char *insegName, int ns, TIME t_cen,
00115                  int *crot, int *crot1, double *lonc);
00116 
00117 /* Check input record, keys, called by checkInput */
00118 int checkInputRec(char *inRSName, char *insegName, int carr,
00119                   double *cdelt1, double *cdelt2, char **ctype1, char **ctype2, int *nx, int *ny);
00120 
00121 /* ################## Working functions ################## */
00122 
00123 /* Top level working wrapper
00124  * Get polar data, remap, smooth, remap back
00125  * Image allocated inside the function
00126  */
00127 int createPole(float **data_raw, float **data, float **img_raw, float **img,
00128                int *dataDims, int *imgDims,
00129                char *inRSName, char *insegName, int ns, int nmax, int crot, int crot1, double *lonc);
00130 
00131 /* Get maps from two CRs, obtain lat/lon and
00132  * dimensions of the polar data
00133  */
00134 int patchMaps(double **map, int *mapDims, double **lons, double **lats, int *dataDims,
00135               char *inRSName, char *insegName, int crot, int crot1, double *lonc);
00136 
00137 /* ################## Output functions ################## */
00138 
00139 /* Output data, free arrays inside and reset to NULL */
00140 int writeOutput(float *data_raw, float *data, float *img_raw, float *img, int *dataDims, int *imgDims,
00141                 char *inRSName, char *outRSName, char *insegName, int year, int ns, int nmax, int crot, double lonc, TIME t_cen, float scale);
00142 
00143 /* Compute mean field between LAT1 and LAT2 and above LAT2 */
00144 void getMean(float *mean1, float *mean2, float *data, int *dims,
00145              int ns, double crpix2, double crval2, double cdelt2);
00146 
00147 /* ################################################# */
00148 /* ################## Main Module ################## */
00149 /* ################################################# */
00150 
00151 // Manually correct keywords
00152 //#define CORRECT_EPHEM 1
00153 
00154 char *module_name = "update_poldb";
00155 
00156 ModuleArgs_t module_args[] =
00157 {
00158     {ARG_INT, "year", NULL, "Year of the data"},
00159     {ARG_INT, "ns", NULL, "North or South, 0 for N, 1 for S"},
00160     {ARG_INT, "nmax", "5", "Max order of Chebychev for smoothing"},
00161     {ARG_STRING, "in", NULL, "Input data series, name only."},
00162     {ARG_STRING, "out", NULL, "Output data series, name only."},
00163     {ARG_FLOAT, "scale", "1.", "Multiplier to scale the polar field"},
00164     {ARG_TIME, "t_cen", "1980.01.01_UTC", "Manual-set center time"},
00165     {ARG_STRING, "inseg", INSEGNAME, "segment name"},
00166     {ARG_FLAG, "o", "", "Flag to overwrite existing data records"},
00167     {ARG_END}
00168 };
00169 
00170 int DoIt(void)
00171 {
00172 
00173     int status = DRMS_SUCCESS;
00174     
00175     /* ============== */
00176     /* Get parameters */
00177     /* ============== */
00178     
00179     char *inRSName = (char *)params_get_str(&cmdparams, "in");
00180     char *outRSName = (char *)params_get_str(&cmdparams, "out");
00181     char *insegName = (char *)params_get_str(&cmdparams, "inseg");
00182     int year = params_get_int(&cmdparams, "year");
00183     int ns = params_get_int(&cmdparams, "ns");
00184     int nmax = params_get_int(&cmdparams, "nmax");
00185     float scale = params_get_float(&cmdparams, "scale");
00186     TIME t_cen = params_get_time(&cmdparams, "t_cen");
00187     int ovwt = params_isflagset(&cmdparams, "o");       // overwrite existing
00188     
00189     /* ================ */
00190     /* Check parameters */
00191     /* ================ */
00192     
00193     SHOW("Check parameters...\n");
00194     
00195     // t_cen will be updated if use default
00196     status = checkParams(year, ns, nmax, scale, &t_cen);
00197     
00198     if (status) DIE("Stop.");
00199     
00200     /* ================== */
00201     /* Open output record */
00202     /* ================== */
00203     
00204     SHOW("Check output data series...\n");
00205     
00206     status = checkOutput(outRSName, year, ns, ovwt);
00207     
00208     if (status) DIE("Stop.");
00209     
00210     /* =================== */
00211     /* Check input records */
00212     /* =================== */
00213     
00214     int crot, crot1;     // map center carr
00215     double lonc;      // center pixel longitude
00216     int dataDims[2], imgDims[2] = {IMGDIM, IMGDIM};
00217     
00218     SHOW("Check input data series...\n");
00219     
00220     // Determine maps to use, and array dimensions
00221     status = checkInput(inRSName, insegName, ns, t_cen,
00222                         &crot, &crot1, &lonc);
00223     
00224     if (status) DIE("Stop.");
00225 
00226     /* ===================== */
00227     /* Get data, then smooth */
00228     /* ===================== */
00229     
00230     // data_raw: cutout of polar data
00231     // data: smoothed polar data, used for filling
00232     // img: intermediate, projected polar view
00233     float *data_raw = NULL, *data = NULL, *img_raw = NULL, *img = NULL;
00234     
00235     SHOW("Processing...\n");
00236     
00237     // All the work!
00238     status = createPole(&data_raw, &data, &img_raw, &img,
00239                         dataDims, imgDims,
00240                         inRSName, insegName, ns, nmax, crot, crot1, &lonc);
00241     
00242     if (status) {
00243         FREE_ARR(data_raw); FREE_ARR(data); FREE_ARR(img); FREE_ARR(img_raw);
00244         DIE("Stop\n");
00245     }
00246     
00247     /* ============= */
00248     /* Output record */
00249     /* ============= */
00250     
00251     SHOW("Writing output...\n");
00252     
00253     // Arrays freed inside if success
00254     status = writeOutput(data_raw, data, img_raw, img, dataDims, imgDims,
00255                          inRSName, outRSName, insegName, year, ns, nmax, crot, lonc, t_cen, scale);
00256     
00257     if (status) {
00258         FREE_ARR(data_raw); FREE_ARR(data);
00259         FREE_ARR(img_raw); FREE_ARR(img);
00260         DIE("Stop\n");
00261     }
00262     
00263     /* ============= */
00264     
00265     return DRMS_SUCCESS;
00266     
00267 }   // DoIt
00268 
00269 
00270 // ========================================
00271 
00272 /* Check parameters, set center time */
00273 
00274 int checkParams(int year, int ns, int nmax, float scale, TIME *t_cen_ptr)
00275 {
00276     
00277     if ((year > currentYear()) || (year < 1996)) {
00278         SHOW("Year is out of range. ");  // HMI/MDI only
00279         return 1;
00280     }
00281     if ((ns < 0) || (ns > 1)) {
00282         SHOW("Please set ns to 0 for north, 1 for south. ");
00283         return 1;
00284     }
00285     if (nmax <= 0) {
00286         SHOW("Please set nmax to a positive integer. ");
00287         return 1;
00288     }
00289     if (scale <= 0) {
00290         SHOW("Please set scale to be positive. ");
00291         return 1;
00292     }
00293     printf("  Update %s polar field database for year %4d.\n",
00294            ((ns) ? "southern" : "northern"), year);
00295     
00296     // Check time, reject time that differs from the default by 30 days
00297     
00298     char t_cen_str[100];
00299     TIME t_cen_tmp;
00300 
00301     snprintf(t_cen_str, 100, "%4d%s", year, ((ns) ? S_DATE : N_DATE));
00302     t_cen_tmp = sscan_time(t_cen_str);      // Default, 09.08 for N, 03.06 for S
00303     
00304     int useTCen = (year == time2Year(*t_cen_ptr)) &&
00305                   ((t_cen_tmp - *t_cen_ptr) < (30. * SECSINDAY));   // Compare with default
00306     
00307     // Final t_cen
00308     
00309     if (!useTCen) {                 // Default center time
00310         *t_cen_ptr = t_cen_tmp;     // Update t_cen in main()
00311         printf("  Use default center time: %s\n", t_cen_str); fflush(stdout);
00312     } else {                        // User specified center time
00313         sprint_time(t_cen_str, *t_cen_ptr, "UTC", 0);
00314         printf("  Use user-specified center time: %s\n", t_cen_str); fflush(stdout);
00315     }
00316     
00317     return 0;
00318     
00319 }
00320 
00321 
00322 // ========================================
00323 
00324 /* Check output availability */
00325 
00326 int checkOutput(char *outRSName, int year, int ns, int ovwt)
00327 {
00328 
00329     int status = 0;
00330     
00331     char outQuery[100];
00332     snprintf(outQuery, 100, "%s[%4d][%1d]", outRSName, year, ns);
00333     
00334     DRMS_RecordSet_t *outRS = NULL;
00335     outRS = drms_open_records(drms_env, outQuery, &status);
00336     if (status) {
00337         if (outRS) drms_close_records(outRS, DRMS_FREE_RECORD);
00338         SHOW("Output record set error. ");
00339         return status;
00340     }
00341      
00342     if (outRS->n != 0 && !ovwt) {
00343         if (outRS) drms_close_records(outRS, DRMS_FREE_RECORD);
00344         SHOW("Output record exists, overwriting not requested. ");
00345         return 1;
00346     }
00347     
00348     return 0;
00349     
00350 }
00351 
00352 
00353 // ========================================
00354 
00355 /* Check input availability, determine crot, etc */
00356 
00357 int checkInput(char *inRSName, char *insegName, int ns, TIME t_cen,
00358                   int *crot, int *crot1, double *lonc)
00359 {
00360  
00361     int status = 0;
00362     int carr, carr1;
00363     
00364     // Center lon, carr_num for t_cen
00365     
00366     double lon_cen, lat_cen;
00367     HeliographicLocation(t_cen, &carr, &lon_cen, &lat_cen);
00368     printf("  Center pixel: crot=%4d, lon_cen=%f\n", carr, lon_cen); fflush(stdout);
00369     
00370     // Input query for crot, check keywords
00371 
00372     int nx = 0, ny = 0;     // output array size
00373     double cdelt1 = DRMS_MISSING_DOUBLE, cdelt2 = DRMS_MISSING_DOUBLE;
00374     char *ctype1 = NULL, *ctype2 = NULL;
00375     
00376     status = checkInputRec(inRSName, insegName, carr, &cdelt1, &cdelt2, &ctype1, &ctype2, &nx, &ny);
00377     if (status) return status;
00378     
00379     // Generally requires two frames, since lon_cen is not at map center
00380     // One full map only when the center longitude is within half of cdelt1 of 180-cdelt1/2.
00381     
00382     int x_off = round((lon_cen - (180. - fabs(cdelt1) / 2.)) / fabs(cdelt1));
00383     carr1 = (x_off == 0) ? carr : ((x_off < 0) ? (carr + 1) : (carr - 1));
00384     
00385     // Input query for a second map, crot+/-1
00386     
00387     int nx_1 = 0, ny_1 = 0;     // output array size
00388     double cdelt1_1 = DRMS_MISSING_DOUBLE, cdelt2_1 = DRMS_MISSING_DOUBLE;
00389     char *ctype1_1 = NULL, *ctype2_1 = NULL;
00390     
00391     status = checkInputRec(inRSName, insegName, carr1, &cdelt1_1, &cdelt2_1, &ctype1_1, &ctype2_1, &nx_1, &ny_1);
00392     if (status) return status;
00393     
00394     // Map dimensions and types must match
00395     
00396     if (carr != carr1) {
00397         if (isnan(cdelt1) || isnan (cdelt2) || isnan(cdelt1_1) || isnan(cdelt2_1) ||
00398             !ctype1 || !ctype2 || !ctype1_1 || !ctype2_1 ||
00399             fabs(cdelt1 - cdelt1_1) > EPSL ||
00400             fabs(cdelt2 - cdelt2_1) > EPSL ||
00401             strcmp(ctype1, ctype1_1) != 0 ||
00402             strcmp(ctype2, ctype2_1) != 0 ||
00403             (strcmp(ctype2, "CRLT-CEA") != 0 && strcmp(ctype2, "Sine Latitude") != 0) ||  // has to be sine-lat
00404             !nx || !ny || !nx_1 || !ny_1 ||
00405             nx != nx_1 || ny != ny_1)
00406         {
00407             SHOW("  Keywords of two input records don't match.\n");
00408             FREE_ARR(ctype1); FREE_ARR(ctype2);
00409             FREE_ARR(ctype1_1); FREE_ARR(ctype2_1);
00410             return 1;
00411         }
00412     }
00413     
00414     // Output
00415     
00416     *crot = carr;
00417     *crot1 = carr1;
00418     *lonc = lon_cen;
00419     
00420     // Clean up
00421     
00422     FREE_ARR(ctype1); FREE_ARR(ctype2);
00423     FREE_ARR(ctype1_1); FREE_ARR(ctype2_1);
00424     
00425     return 0;
00426 }
00427 
00428 
00429 // ========================================
00430 
00431 /* Check input record, keys, called by checkInput */
00432 
00433 int checkInputRec(char *inRSName, char *insegName, int carr,
00434                   double *cdelt1, double *cdelt2, char **ctype1, char **ctype2, int *nx, int *ny)
00435 {
00436     
00437     int status = 0;
00438     
00439     char inQuery[100];
00440     snprintf(inQuery, 100, "%s[%4d]", inRSName, carr);
00441     printf("  Check %s\n", inQuery);
00442     
00443     DRMS_RecordSet_t *inRS = NULL;
00444     inRS = drms_open_records(drms_env, inQuery, &status);
00445     if (status || inRS->n == 0) {
00446         if (inRS) drms_close_records(inRS, DRMS_FREE_RECORD);
00447         SHOW("  Input record error.\n");
00448         return ((status) ? status : 1);
00449     }
00450     
00451     DRMS_Record_t *inRec = inRS->records[0];
00452     *cdelt1 = drms_getkey_double(inRec, "CDELT1", &status);
00453     *cdelt2 = drms_getkey_double(inRec, "CDELT2", &status);
00454     *ctype1 = drms_getkey_string(inRec, "CTYPE1", &status);
00455     *ctype2 = drms_getkey_string(inRec, "CTYPE2", &status);
00456     
00457     // Check segment & array
00458     
00459     DRMS_Segment_t *inSeg = drms_segment_lookup(inRec, insegName);
00460     DRMS_Array_t *inArray = drms_segment_read(inSeg, DRMS_TYPE_FLOAT, &status);
00461     if (status) {
00462         drms_close_records(inRS, DRMS_FREE_RECORD);
00463         SHOW(" Array read error.\n");
00464         return status;
00465     }
00466     
00467     // Determine map size
00468     
00469     *nx = inArray->axis[0]; *ny = inArray->axis[1];
00470     
00471     // Clean up
00472     
00473     drms_free_array(inArray);
00474     drms_close_records(inRS, DRMS_FREE_RECORD);
00475     
00476     printf("    CDELT1=%f, CDELT2=%f, CTYPE1=%s, CTYPE2=%s\n",
00477            *cdelt1, *cdelt2, *ctype1, *ctype2);
00478     
00479     return 0;
00480     
00481 }
00482 
00483 
00484 // ========================================
00485 
00486 /* Top level working wrapper
00487  * Get polar data, remap, smooth, remap back
00488  * Image allocated inside the function
00489  */
00490 
00491 int createPole(float **data_raw, float **data, float **img_raw, float **img,
00492                int *dataDims, int *imgDims,
00493                char *inRSName, char *insegName, int ns, int nmax, int crot, int crot1, double *lonc)
00494 {
00495  
00496     int status = 0;
00497     
00498     // Open maps
00499     
00500     double *map = NULL;      // full patched map
00501     double *lons = NULL, *lats = NULL;     // longitude, latitude
00502     int mapDims[2];         // full map size
00503     
00504     status = patchMaps(&map, mapDims, &lons, &lats, dataDims,       // All these are updated
00505                        inRSName, insegName, crot, crot1, lonc);
00506     if (status) return status;
00507     
00508     // Allocate arrays
00509     
00510     int nx = dataDims[0], ny = dataDims[1], nxny = nx * ny;
00511     int imgDim2 = imgDims[0] * imgDims[1];
00512     
00513     printf("  nx=%d, ny=%d\n", nx, ny);
00514     
00515     *data_raw = (float *) (malloc(nxny * sizeof(float)));
00516     *data = (float *) (malloc(nxny * sizeof(float)));
00517     *img_raw = (float *) (malloc(imgDim2 * sizeof(float)));
00518     *img = (float *) (malloc(imgDim2 * sizeof(float)));
00519     
00520     double *image = (double *) (malloc(imgDim2 * sizeof(double)));  // operate in double
00521     
00522     // data_raw: direct cutout of the polar region
00523     
00524     int offset = (ns) ? 0 : (mapDims[0] * (mapDims[1] - ny));      // south start from 0
00525     
00526     for (int count = 0; count < nxny; count++) {
00527         (*data_raw)[count] = map[count+offset];
00528     }
00529 
00530     // Stereographic projection to polar view
00531     
00532     synop2pole(ns, 1, image, imgDims, map, mapDims, lons, lats, (LAT0*DTOR));   // in polar_remap.c
00533     
00534     // Smoothing
00535     
00536     for (int count = 0; count < imgDim2; count++) (*img_raw)[count] = image[count];     // copy out raw
00537     fitimage(image, imgDims[0], imgDims[1], nmax, 1);
00538     for (int count = 0; count < imgDim2; count++) (*img)[count] = image[count];     // copy out
00539     
00540     // Stereographic projection back to Carrington map, saved back to map
00541     
00542     pole2synop(ns, 1, image, imgDims, map, mapDims, lons, lats, (LAT0*DTOR));   // in polar_remap.c
00543     
00544     // data: smoothed polar region
00545     
00546     for (int count = 0; count < nxny; count++) {
00547         (*data)[count] = map[count+offset];
00548     }
00549     
00550     // Clean up
00551     
00552     FREE_ARR(image);
00553     FREE_ARR(map);
00554     FREE_ARR(lons); FREE_ARR(lats);
00555     
00556     return 0;
00557     
00558 }
00559 
00560 // ========================================
00561 
00562 /* Get maps from two CRs, obtain lat/lon and 
00563  * dimensions of the polar data
00564  * lonc will be updated according to resolution
00565  */
00566 
00567 int patchMaps(double **map, int *mapDims, double **lons, double **lats, int *dataDims,
00568               char *inRSName, char *insegName, int crot, int crot1, double *lonc)
00569 {
00570     
00571     int status = 0;
00572     char inQuery[100], inQuery1[100];
00573     DRMS_RecordSet_t *inRS = NULL, *inRS1 = NULL;
00574     DRMS_Record_t *inRec = NULL, *inRec1 = NULL;
00575     DRMS_Segment_t *inSeg = NULL, *inSeg1 = NULL;
00576     DRMS_Array_t *inArray = NULL, *inArray1 = NULL;
00577     
00578     // Open CROT
00579 
00580     snprintf(inQuery, 100, "%s[%4d]", inRSName, crot);
00581     inRS = drms_open_records(drms_env, inQuery, &status);
00582     if (status) {
00583         SHOW(" Getting data error.\n");
00584         return status;
00585     }
00586     inRec = inRS->records[0];
00587     inSeg = drms_segment_lookup(inRec, insegName);
00588     inArray = drms_segment_read(inSeg, DRMS_TYPE_FLOAT, &status);
00589     if (status) {
00590         SHOW(" Getting data error.\n");
00591         if (inArray) drms_free_array(inArray);
00592         drms_close_records(inRS, DRMS_FREE_RECORD);
00593         return status;
00594     }
00595     
00596     // Open CROT1 uf necessary
00597     
00598     if (crot != crot1) {
00599         snprintf(inQuery1, 100, "%s[%4d]", inRSName, crot1);
00600         inRS1 = drms_open_records(drms_env, inQuery1, &status);
00601         if (status) {
00602             SHOW(" Getting data error.\n");
00603             drms_free_array(inArray);
00604             drms_close_records(inRS, DRMS_FREE_RECORD);
00605             return status;
00606         }
00607         inRec1 = inRS1->records[0];
00608         inSeg1 = drms_segment_lookup(inRec1, insegName);
00609         inArray1 = drms_segment_read(inSeg1, DRMS_TYPE_FLOAT, &status);
00610         if (status) {
00611             SHOW(" Getting data error.\n");
00612             drms_free_array(inArray);
00613             if (inArray1) drms_free_array(inArray1);
00614             drms_close_records(inRS, DRMS_FREE_RECORD);
00615             drms_close_records(inRS1, DRMS_FREE_RECORD);
00616             return status;
00617         }
00618     }
00619     
00620     // Lon/lat of full map
00621     
00622     int xsz = inArray->axis[0], ysz = inArray->axis[1];
00623     mapDims[0] = xsz; mapDims[1] = ysz;
00624     
00625     *lons = (double *) (malloc(xsz * sizeof(double)));
00626     *lats = (double *) (malloc(ysz * sizeof(double)));
00627     
00628     double cdelt1 = drms_getkey_double(inRec, "CDELT1", &status);
00629     double cdelt2 = drms_getkey_double(inRec, "CDELT2", &status);
00630     double crpix1 = drms_getkey_double(inRec, "CRPIX1", &status);
00631     double crpix2 = drms_getkey_double(inRec, "CRPIX2", &status);
00632     double crval1 = drms_getkey_double(inRec, "CRVAL1", &status);   // Carrington time
00633 #ifdef CORRECT_EPHEM
00634     crpix1 += 1;        // crpix2 += 0.5; crval1 += fabs(cdelt1) / 2.0;     // no need now, longitude start from 0.1
00635 #endif
00636     double crval2 = drms_getkey_double(inRec, "CRVAL2", &status);
00637     printf("cdelt1=%f, crpix1=%f, crval1=%f\n", cdelt1, crpix1, crval1);
00638     
00639     // Outmap relative coords
00640     
00641     for (int col = 0; col < xsz; col++) {
00642         (*lons)[col] = (col + 1) * fabs(cdelt1) * DTOR - M_PI;      // -179.9, 180; changed to multiple of 0.1 on Apr 30 XS
00643     }
00644     for (int row = 0; row < ysz; row++) {
00645         (*lats)[row] = asin((row - (crpix2 - 1.)) * cdelt2 + crval2);
00646     }
00647     
00648     // Size of polar data
00649     
00650     int row = 0;
00651     while (row < (ysz / 2)) {
00652         if (fabs((*lats)[row++] / DTOR) < LAT0) break;
00653     }
00654     
00655     dataDims[0] = xsz; dataDims[1] = row;
00656     
00657     printf("  Orig map size: %d x %d\n", xsz, ysz);
00658     printf("  Polar data size: %d x %d\n", dataDims[0], dataDims[1]);
00659     
00660     // Map
00661     
00662     int xysz = xsz * ysz;
00663     *map = (double *) (malloc(xysz * sizeof(double)));      // 3600x1440
00664     
00665     float *synop = (float *) inArray->data;
00666     
00667     if (crot1 == crot) {
00668 
00669         // Need only one map
00670         for (int count = 0; count < xysz; count++) {(*map)[count] = synop[count];}      // copy everything
00671 
00672     } else {
00673         
00674         float *synop1 = (float *) inArray1->data;
00675         
00676         // longitude of map center
00677         
00678         double crlon = crot * (double)360. - crval1;            // generally 180
00679         float colc = (*lonc - crlon) / fabs(cdelt1) + crpix1 - 1;   // column for requested lonc, start from 0
00680         colc = (xsz % 2) ? round(colc) : (floor(colc) + 0.5);
00681         *lonc = (colc + 1 - crpix1) * fabs(cdelt1) + crlon;     // actual center longitude, modified
00682         printf("  colc=%f, lonc=%f\n", colc, *lonc);
00683         
00684         // Copy & patch data
00685         
00686         int col0;
00687         if (crot1 > crot) {     // crot on the right
00688             col0 = round(colc + (xsz - 1) / 2.);        // right edge, col0 included
00689             for (int row = 0; row < ysz; row++) {
00690                 int offset_s = row * xsz;
00691                 int offset = row * xsz + xsz - 1 - col0;
00692                 for (int col = 0; col <= col0; col++) {(*map)[offset + col] = synop[offset_s + col];}        // 0, col0
00693                 int offset1 = row * xsz - col0 - 1;
00694                 for (int col = col0 + 1; col < xsz; col++) {(*map)[offset1 + col] = synop1[offset_s + col];}    // col0 + 1, xsz-1
00695             }
00696         } else {                // crot on the left
00697             col0 = round(colc - (xsz - 1) / 2.);        // left edge, col0 included
00698             for (int row = 0; row < ysz; row++) {
00699                 int offset_s = row * xsz;
00700                 int offset = row * xsz - col0;
00701                 for (int col = col0; col < xsz; col++) {(*map)[offset + col] = synop[offset_s + col];}      // col0, xsz-1
00702                 int offset1 = row * xsz + xsz - col0;
00703                 for (int col = 0; col < col0; col++) {(*map)[offset1 + col] = synop1[offset_s + col];}      // 0, col0-1
00704             }
00705         }
00706         printf("  col0=%d, lonc=%f, map[0,100]=%f, map[0,101]=%f\n", col0, *lonc, (*map)[100l*xsz], (*map)[101l*xsz]);
00707     }
00708     
00709     // Clean up
00710     
00711     drms_free_array(inArray);
00712     drms_close_records(inRS, DRMS_FREE_RECORD);
00713     if (crot1 != crot) {
00714         drms_free_array(inArray1);
00715         drms_close_records(inRS1, DRMS_FREE_RECORD);
00716     }
00717     
00718     return 0;
00719     
00720 }
00721 
00722 // ========================================
00723 
00724 // ========================================
00725 
00726 /* Output data, free arrays inside and reset to NULL */
00727 
00728 int writeOutput(float *data_raw, float *data, float *img_raw, float *img, int *dataDims, int *imgDims,
00729                 char *inRSName, char *outRSName, char *insegName, int year, int ns, int nmax, int crot, double lonc, TIME t_cen, float scale)
00730 {
00731     
00732     int status = 0;
00733     
00734     // Create output
00735     
00736     DRMS_Record_t *outRec = drms_create_record(drms_env, outRSName, DRMS_PERMANENT, &status);
00737     if (status) {
00738         SHOW("  Creating output record failed.\n");
00739         return status;
00740     }
00741     
00742     // Segment
00743     
00744     DRMS_Segment_t *outSeg_data = NULL, *outSeg_data_raw = NULL;
00745     DRMS_Segment_t *outSeg_img = NULL, *outSeg_img_raw = NULL;
00746     outSeg_data = drms_segment_lookup(outRec, "data");          // smoothed
00747     outSeg_data_raw = drms_segment_lookup(outRec, "data_raw");      // original
00748     outSeg_img = drms_segment_lookup(outRec, "image");          // polar view
00749     outSeg_img_raw = drms_segment_lookup(outRec, "image_raw");          // polar view original
00750     if (!outSeg_data || !outSeg_data_raw || !outSeg_img || !outSeg_img_raw) {
00751         drms_close_record(outRec, DRMS_FREE_RECORD);
00752         SHOW("  Creating output segment failed.\n");
00753         return 1;
00754     }
00755     
00756     // Arrays
00757     
00758     DRMS_Array_t *outArray_data = drms_array_create(DRMS_TYPE_FLOAT, 2, dataDims, data, &status);
00759     if (status) {SHOW("  Creating data array failed.\n"); return 1;}
00760     DRMS_Array_t *outArray_data_raw = drms_array_create(DRMS_TYPE_FLOAT, 2, dataDims, data_raw, &status);
00761     if (status) {SHOW("  Creating data_raw array failed.\n"); return 1;}
00762     DRMS_Array_t *outArray_img = drms_array_create(DRMS_TYPE_FLOAT, 2, imgDims, img, &status);
00763     if (status) {SHOW("  Creating img array failed.\n"); return 1;}
00764     DRMS_Array_t *outArray_img_raw = drms_array_create(DRMS_TYPE_FLOAT, 2, imgDims, img_raw, &status);
00765     if (status) {SHOW("  Creating img array failed.\n"); return 1;}
00766     
00767     for (int i = 0; i < 2; i++) {
00768         outSeg_data->axis[i] = outArray_data->axis[i];
00769         outSeg_data_raw->axis[i] = outArray_data_raw->axis[i];
00770         outSeg_img->axis[i] = outArray_img->axis[i];
00771         outSeg_img_raw->axis[i] = outArray_img_raw->axis[i];
00772     }
00773     
00774     outArray_data->israw = 0;       // always compressed
00775     outArray_data_raw->israw = 0;
00776     outArray_img->israw = 0;
00777     outArray_img_raw->israw = 0;
00778     
00779     outArray_data->bzero = outSeg_data->bzero;
00780     outArray_data_raw->bzero = outSeg_data_raw->bzero;
00781     outArray_img->bzero = outSeg_img->bzero;
00782     outArray_img_raw->bzero = outSeg_img_raw->bzero;
00783     
00784     outArray_data->bscale = outSeg_data->bscale;
00785     outArray_data_raw->bscale = outSeg_data_raw->bscale;
00786     outArray_img->bscale = outSeg_img->bscale;
00787     outArray_img_raw->bscale = outSeg_img_raw->bscale;
00788     
00789     // Write
00790     
00791     status = drms_segment_write(outSeg_data, outArray_data, 0);
00792     if (status) {SHOW("  Writing data array failed.\n"); return 1;}
00793     status = drms_segment_write(outSeg_data_raw, outArray_data_raw, 0);
00794     if (status) {SHOW("  Writing data_raw array failed.\n"); return 1;}
00795     status = drms_segment_write(outSeg_img, outArray_img, 0);
00796     if (status) {SHOW("  Writing img array failed.\n"); return 1;}
00797     status = drms_segment_write(outSeg_img_raw, outArray_img_raw, 0);
00798     if (status) {SHOW("  Writing img_raw array failed.\n"); return 1;}
00799     
00800     // Keywords
00801     
00802     char inQuery[100];
00803     snprintf(inQuery, 100, "%s[%4d]", inRSName, crot);
00804     DRMS_RecordSet_t *inRS = drms_open_records(drms_env, inQuery, &status);
00805     DRMS_Record_t *inRec = inRS->records[0];
00806     DRMS_Segment_t *inSeg = drms_segment_lookup(inRec, insegName);
00807     int nKeys = ARRLENGTH(synop_keys);
00808     for (int iKey = 0; iKey < nKeys; iKey++) {
00809         drms_copykey(outRec, inRec, synop_keys[iKey]);
00810     }
00811     int ysz = inSeg->axis[1];
00812     double yc = drms_getkey_double(inRec, "CRPIX2", &status);
00813 #ifdef CORRECT_EPHEM
00814     // yc += 0.5;       // Apr 30
00815 #endif
00816     double yref = drms_getkey_double(inRec, "CRVAL2", &status);
00817     
00818     drms_setkey_int(outRec, "YEAR", year);
00819     drms_setkey_int(outRec, "NS", ns);
00820     drms_setkey_int(outRec, "NMAX", nmax);
00821     drms_setkey_int(outRec, "CAR_ROT", crot);
00822     drms_setkey_double(outRec, "CARRTIME", crot * (double)360. - lonc);
00823     drms_setkey_time(outRec, "T_OBS", t_cen);
00824     drms_setkey_float(outRec, "SCALE", scale);
00825     
00826     drms_setkey_double(outRec, "CRPIX1", (1. + dataDims[0]) / 2.);
00827     double crpix2 = (ns) ? yc : (yc - ysz + dataDims[1]);
00828     drms_setkey_double(outRec, "CRPIX2", crpix2);
00829     drms_setkey_double(outRec, "CRVAL1", crot * (double)360. - lonc);
00830     drms_setkey_double(outRec, "CRVAL2", yref);
00831     
00832     char timebuf[1024];
00833     double UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */
00834     sprint_time(timebuf, (double)time(NULL) + UNIX_epoch, "ISO", 0);
00835     drms_setkey_string(outRec, "DATE", timebuf);
00836     
00837     // Mean field between LAT1 and LAT2 and above LAT2
00838     
00839     double cdelt2 = drms_getkey_double(inRec, "CDELT2", &status);
00840     float mean1, mean2;
00841     getMean(&mean1, &mean2, data_raw, dataDims, ns, crpix2, yref, cdelt2);
00842     drms_setkey_float(outRec, "MEAN62", mean1);
00843     drms_setkey_float(outRec, "MEAN75", mean2);
00844     
00845     //
00846     
00847     drms_free_array(outArray_data); data = NULL;
00848     drms_free_array(outArray_data_raw); data_raw = NULL;
00849     drms_free_array(outArray_img); img = NULL;
00850     drms_free_array(outArray_img_raw); img_raw = NULL;
00851     
00852     drms_close_records(inRS, DRMS_FREE_RECORD);
00853     drms_close_record(outRec, DRMS_INSERT_RECORD);
00854 
00855     return 0;
00856     
00857 }
00858 
00859 // =====================
00860 
00861 /* Compute mean field between LAT1 and LAT2 and above LAT2 */
00862 
00863 void getMean(float *mean1, float *mean2, float *data, int *dims,
00864              int ns, double crpix2, double crval2, double cdelt2)
00865 {
00866     
00867     double lat;
00868     int row1, row2;
00869 
00870     for (int row = 0; row < dims[1]; row++) {
00871         lat = fabs(asin((row - crpix2) * cdelt2 + crval2) / DTOR);
00872         if (ns) {       // south
00873             if (lat >= LAT2) row2 = row;
00874             if (lat >= LAT1) row1 = row;
00875         } else {
00876             if (lat <= LAT2) row2 = row;
00877             if (lat <= LAT1) row1 = row;
00878         }
00879     }
00880     printf("  row1=%d, row2=%d\n", row1, row2);
00881     
00882     //
00883     
00884     double sum1 = 0.0, sum2 = 0.0;
00885     int count1 = 0, count2 = 0;
00886     double val;
00887     if (ns) {
00888         for (int row = 0; row < row2; row++) {
00889             for (int col = 0; col < dims[0]; col++) {
00890                 val = data[row * dims[0] + col];
00891                 if (isnan(val)) continue;
00892                 sum2 += val; count2++;
00893             }
00894         }
00895         for (int row = row2; row < row1; row++) {
00896             for (int col = 0; col < dims[0]; col++) {
00897                 val = data[row * dims[0] + col];
00898                 if (isnan(val)) continue;
00899                 sum1 += val; count1++;
00900             }
00901         }
00902     } else {
00903         for (int row = row1; row < row2; row++) {
00904             for (int col = 0; col < dims[0]; col++) {
00905                 val = data[row * dims[0] + col];
00906                 if (isnan(val)) continue;
00907                 sum1 += val; count1++;
00908             }
00909         }
00910         for (int row = row2; row < dims[1] - 1; row++) {
00911             for (int col = 0; col < dims[0]; col++) {
00912                 val = data[row * dims[0] + col];
00913                 if (isnan(val)) continue;
00914                 sum2 += val; count2++;
00915             }
00916         }
00917     }
00918     
00919     //
00920     
00921     *mean1 = (float) (sum1 / count1);
00922     *mean2 = (float) (sum2 / count2);
00923     
00924     printf("  mean1=%f, mean2=%f\n", *mean1, *mean2);
00925     
00926 }
proj/mag/polcorr/apps/update_poldb.c
