JSOC_Documentation: proj/cgem/prep/apps/cgem

00001 /*
00002  *  cgem_map.c
00003  *
00004  *  This module maps a scaler map of one projection to another
00005  *  Written as prototype for functions needed for cgem_prep
00006  *  One application is to map a Plate-Carree Br map into Mercater and back
00007  *
00008  *  Author:
00009  *      Xudong Sun
00010  *
00011  *  Version:
00012  *              v0.0 Jul 18 2016
00013  *  Notes:
00014  *      v0.0
00015  *
00016  *  Example Calls:
00017  *      > cgem_map "in=su_xudong.B_720s_carree[377][2011.02.15_00:00]" "seg=Bz" "out=su_xudong.B_720s_mercator" "map=Mercator" "cols=600" "rows=428" -a
00018  *
00019  */
00020 
00021 
00022 #include <stdio.h>
00023 #include <stdlib.h>
00024 #include <time.h>
00025 #include <sys/time.h>
00026 #include <math.h>
00027 #include <string.h>
00028 #include "jsoc_main.h"
00029 #include "astro.h"
00030 #include "cartography_cgem.c"
00031 #include "finterpolate.h"
00032 
00033 #define PI              (M_PI)
00034 #define RADSINDEG       (PI/180.)
00035 #define RAD2ARCSEC      (648000./M_PI)
00036 #define SECINDAY        (86400.)
00037 #define FOURK           (4096)
00038 #define FOURK2          (16777216)
00039 
00040 #define ARRLENGTH(ARR) (sizeof(ARR) / sizeof(ARR[0]))
00041 
00042 // Some other things
00043 #ifndef MIN
00044 #define MIN(a,b) (((a)<(b)) ? (a) : (b))
00045 #endif
00046 #ifndef MAX
00047 #define MAX(a,b) (((a)>(b)) ? (a) : (b))
00048 #endif
00049 
00050 #define DIE(msg) {fflush(stdout); fprintf(stderr,"%s, status=%d\n", msg, status); return(status);}
00051 #define SHOW(msg) {printf("%s", msg); fflush(stdout);}
00052 
00053 #define kNotSpecified "Not Specified"
00054 #define dpath    "/home/jsoc/cvs/Development/JSOC"
00055 
00056 // Macros for WCS transformations.  assume crpix1, crpix2 = CRPIX1, CRPIX2, sina,cosa = sin and cos of CROTA2 resp.
00057 // and crvalx and crvaly are CRVAL1 and CRVAL2, cdelt = CDELT1 == CDELT2, then
00058 // PIX_X and PIX_Y are CCD pixel addresses, WX and WY are arc-sec W and N on the Sun from disk center.
00059 #define PIX_X(wx,wy) ((((wx-crvalx)*cosa + (wy-crvaly)*sina)/cdelt)+crpix1)
00060 #define PIX_Y(wx,wy) ((((wy-crvaly)*cosa - (wx-crvalx)*sina)/cdelt)+crpix2)
00061 #define WX(pix_x,pix_y) (((pix_x-crpix1)*cosa - (pix_y-crpix2)*sina)*cdelt+crvalx)
00062 #define WY(pix_x,pix_y) (((pix_y-crpix2)*cosa + (pix_x-crpix1)*sina)*cdelt+crvaly)
00063 
00064 // Ephemeris information
00065 struct ephemeris {
00066     double disk_lonc, disk_latc;
00067     double disk_xc, disk_yc;
00068     double rSun, asd, pa;
00069 };
00070 
00071 // Mapping information
00072 struct reqInfo {
00073     double xc, yc;
00074     int ncol, nrow;
00075     double dx, dy;
00076     int proj;
00077 };
00078 
00079 /* ========================================================================================================== */
00080 
00081 char *mapName[] = {"PlateCarree", "Cassini-Soldner", "Mercator",
00082     "LambertCylindrical", "Sanson-Flamsteed", "gnomonic", "Postel",
00083     "stereographic", "orthographic", "LambertAzimuthal"};
00084 char *wcsCode[] = {"CAR", "CAS", "MER", "CEA", "GLS", "TAN", "ARC", "STG",
00085     "SIN", "ZEA"};
00086 
00087 /* ========================================================================================================== */
00088 
00089 /* ========================================================================================================== */
00090 
00091 char *module_name = "cgem_map";
00092 
00093 ModuleArgs_t module_args[] =
00094 {
00095     {ARG_STRING, "in", kNotSpecified, "Input series."},
00096     {ARG_STRING, "out", kNotSpecified, "Output series."},
00097     {ARG_STRING, "seg", kNotSpecified, "Iutput segment."},
00098     {ARG_NUME, "map", "carree", "Projetion method, carree by default.",
00099         "carree, Cassini, Mercator, cyleqa, sineqa, gnomonic, Postel, stereographic, orthographic, Lambert"},
00100     {ARG_FLOAT, "xref", "0", "Reference output map x center, ususally Carrington lon."},
00101     {ARG_FLOAT, "yref", "0", "Reference output map y center, ususally Carrington lon."},
00102     {ARG_FLOAT, "dx", "0", "X pixel size, unit depending on projection (default deg)."},
00103     {ARG_FLOAT, "dy", "0", "Y pixel size, unit depending on projection (default deg)."},
00104     {ARG_FLAG, "a", "", "Automatic, i.e. adopting original center and scale, override xyref, dxy."},
00105     {ARG_INT, "cols", "500", "Columns of output."},
00106     {ARG_INT, "rows", "500", "Rows of output."},
00107     {ARG_END}
00108 };
00109 
00110 int DoIt(void)
00111 {
00112     
00113     int status = DRMS_SUCCESS;
00114     
00115     /* Get data series */
00116     
00117     char *inQuery = NULL, *segQuery = NULL;      // Input query
00118     char *outQuery = NULL;                      // Output query
00119     
00120     inQuery = (char *) params_get_str(&cmdparams, "in");
00121     segQuery = (char *) params_get_str(&cmdparams, "seg");
00122     outQuery = (char *) params_get_str(&cmdparams, "out");
00123     
00124     if (!strcmp(inQuery, kNotSpecified) ||
00125         !strcmp(segQuery, kNotSpecified) ||
00126         !strcmp(outQuery, kNotSpecified)) {
00127         DIE("Input/output not available");
00128     }
00129     
00130     /* Get arguments */
00131     
00132     int automap = params_isflagset(&cmdparams, "a");
00133     
00134     struct reqInfo req;
00135     req.xc = params_get_float(&cmdparams, "xref") * RADSINDEG;
00136     req.yc = params_get_float(&cmdparams, "yref") * RADSINDEG;
00137     req.ncol = params_get_int(&cmdparams, "cols");
00138     req.nrow = params_get_int(&cmdparams, "rows");
00139     req.dx = params_get_float(&cmdparams, "dx") * RADSINDEG;        // deg to rad, for now
00140     req.dy = params_get_float(&cmdparams, "dy") * RADSINDEG;
00141     req.proj = params_get_int(&cmdparams, "map");
00142     
00143     if (req.ncol <= 2 || req.ncol > 4000 ||
00144         req.nrow <= 2 || req.nrow > 4000) {
00145         DIE("requested map too small/too big");
00146     }
00147 
00148     /* Input Data */
00149     
00150     DRMS_RecordSet_t *inRS = drms_open_records(drms_env, inQuery, &status);
00151     int nrecs = inRS->n;
00152     if (status || nrecs == 0 || !inRS) DIE("Input data series error");
00153     DRMS_Segment_t *inSeg = drms_segment_lookup(inRS->records[0], segQuery);        // check segment
00154     if (!inSeg) {
00155         drms_close_records(inRS, DRMS_FREE_RECORD);
00156         DIE("Requested segment not available");
00157     }
00158     
00159     /* Output Data */
00160     
00161     DRMS_RecordSet_t *outRS = drms_create_records(drms_env, nrecs, outQuery, DRMS_PERMANENT, &status);
00162     if (status || !outRS) {
00163         drms_close_records(outRS, DRMS_FREE_RECORD);
00164         DIE("Error creating output series");
00165     }
00166     
00167     /* Loop */
00168     
00169     for (int irec = 0; irec < nrecs; irec++) {
00170         
00171         DRMS_Record_t *inRec = inRS->records[irec];
00172         DRMS_Record_t *outRec = outRS->records[irec];
00173         printf("Processing rec #%d of %d:\n", irec+1, nrecs);
00174         
00175         // Determine center, scale, xysize, etc.
00176         
00177         double xc_o, yc_o, dx_o, dy_o, colc_o, rowc_o;  // original
00178         double crln_obs;
00179         int ncol_o, nrow_o, proj_o;
00180         char proj_str[4], *ctype1;
00181         double xc, yc, dx, dy;          // requested
00182         int ncol, nrow, proj;
00183         
00184         // input projection, no error checking for now
00185         xc_o = drms_getkey_double(inRec, "CRVAL1", &status) * RADSINDEG;        // lonc
00186         yc_o = drms_getkey_double(inRec, "CRVAL2", &status) * RADSINDEG;        // latc
00187         crln_obs = drms_getkey_double(inRec, "CRLN_OBS", &status) * RADSINDEG;
00188         xc_o -= crln_obs;       // Stonyhurst
00189         colc_o = drms_getkey_double(inRec, "CRPIX1", &status) - 1;
00190         rowc_o = drms_getkey_double(inRec, "CRPIX2", &status) - 1;
00191         dx_o = drms_getkey_double(inRec, "CDELT1", &status) * RADSINDEG;
00192         dy_o = drms_getkey_double(inRec, "CDELT2", &status) * RADSINDEG;
00193         ctype1 = drms_getkey_string(inRec, "CTYPE1", &status);
00194         strncpy(proj_str, ctype1+5, 3);
00195         for (proj_o = 0; proj_o++; proj_o < ARRLENGTH(wcsCode)) {
00196             if (strcmp(proj_str, wcsCode[proj_o])== 0) break;
00197         }
00198         if (proj_o >= ARRLENGTH(wcsCode) || status) {
00199             printf("Erroneous geometry, image skipped.\n");
00200             continue;
00201         } else {
00202             proj_o -= 1;
00203         }
00204         if (automap) {
00205             xc = xc_o; yc = yc_o;
00206             dx = dx_o; dy = dy_o;
00207             // determin maximum size
00208             ;
00209         } else {
00210             xc = req.xc - crln_obs; yc = req.yc;
00211             dx = req.dx; dy = req.dy;
00212         }
00213         ncol = req.ncol; nrow = req.nrow;
00214         proj = req.proj;
00215         printf("xc_o=%f, yc_o=%f, dx_o=%f, dy_o=%f, proj_o=%s\n",
00216                xc_o/RADSINDEG, yc_o/RADSINDEG, dx_o/RADSINDEG, dy_o/RADSINDEG, mapName[proj_o]);
00217         printf("xc=%f, yc=%f, dx=%f, dy=%f, proj=%s\n",
00218                xc/RADSINDEG, yc/RADSINDEG, dx/RADSINDEG, dy/RADSINDEG, mapName[proj]);
00219         
00220         // Input image
00221         
00222         inSeg = drms_segment_lookup(inRec, segQuery);
00223         DRMS_Array_t *inArray = drms_segment_read(inSeg, DRMS_TYPE_FLOAT, &status);
00224         if (status) return 1;
00225         float *map_in = (float *)inArray->data;
00226         ncol_o = inSeg->axis[0]; nrow_o = inSeg->axis[1];
00227         
00228         // Coordinates
00229         
00230         int npix = ncol * nrow;
00231         float *col_o = (float *) (malloc(npix * sizeof(float)));
00232         float *row_o = (float *) (malloc(npix * sizeof(float)));
00233         double rowc = (nrow - 1.) / 2., colc = (ncol - 1.) / 2.;
00234         printf("colc=%f, rowc=%f, colc_o=%f, rowc_o=%f\n", colc, rowc, colc_o, rowc_o);
00235         
00236         // Convert
00237         
00238         printf("Mapping..."); fflush(stdout);
00239         double x, y, x_o, y_o, lon, lat;
00240         int idx;
00241         for (int row = 0; row < nrow; row++) {
00242             y = (row - rowc) * dy;
00243             for (int col = 0; col < ncol; col++) {
00244                 x = (col - colc) * dx;
00245                 idx = row * ncol + col;
00246                 if (plane2sphere(x, y, yc, xc, &lat, &lon, proj)) {     // requested
00247                     col_o[idx] = DRMS_MISSING_FLOAT;
00248                     row_o[idx] = DRMS_MISSING_FLOAT;
00249                     continue;
00250                 }
00251                 if (sphere2plane(lat, lon, yc_o, xc_o, &x_o, &y_o, proj_o)) {       // orig
00252                     col_o[idx] = DRMS_MISSING_FLOAT;
00253                     row_o[idx] = DRMS_MISSING_FLOAT;
00254                     continue;
00255                 }
00256                 col_o[idx] = x_o / dx_o + colc_o;      // index in input images
00257                 row_o[idx] = y_o / dy_o + rowc_o;
00258             }
00259         }
00260         printf(" done\n"); fflush(stdout);
00261         
00262         // Interpolate
00263         
00264         printf("Interpolate..."); fflush(stdout);
00265         float *map_out = (float *) (malloc(npix * sizeof(float)));      // output array
00266         struct fint_struct pars;
00267         init_finterpolate_linear(&pars, 0.);        // bi-linear
00268         
00269         finterpolate(&pars, map_in, col_o, row_o, map_out,
00270                      ncol_o, nrow_o, ncol_o, ncol, nrow, ncol, DRMS_MISSING_FLOAT);
00271         
00272         free(col_o); free(row_o);       // Clean up
00273         drms_free_array(inArray);
00274         printf(" done\n"); fflush(stdout);
00275         
00276         // Output
00277         
00278          printf("Write data..."); fflush(stdout);
00279         int outDims[2] = {ncol, nrow};
00280         DRMS_Segment_t *outSeg = drms_segment_lookupnum(outRec, 0);
00281         DRMS_Array_t *outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, outDims, map_out, &status);
00282         if (status) {
00283             SHOW("Output array error");
00284             free(map_out);
00285             continue;
00286         }
00287         outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1];
00288         outArray->israw = 0;        // always compressed
00289         outArray->bzero = outSeg->bzero;
00290         outArray->bscale = outSeg->bscale;
00291         
00292         status = drms_segment_write(outSeg, outArray, 0);
00293         drms_free_array(outArray);
00294         if (status) {
00295             SHOW("Output array write error");
00296             continue;
00297         }
00298         
00299         // Keywords
00300         
00301         drms_copykeys(outRec, inRec, 0, kDRMS_KeyClass_Explicit);     // copy all keys
00302         
00303         TIME val, trec, tnow, UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */
00304         tnow = (double)time(NULL);
00305         tnow += UNIX_epoch;
00306         drms_setkey_time(outRec, "DATE", tnow);
00307 
00308         drms_setkey_float(outRec, "CRPIX1", ncol/2. + 0.5);
00309         drms_setkey_float(outRec, "CRPIX2", nrow/2. + 0.5);
00310         
00311         drms_setkey_float(outRec, "CRVAL1", (xc + crln_obs) / RADSINDEG);
00312         drms_setkey_float(outRec, "CRVAL2", yc / RADSINDEG);
00313         drms_setkey_float(outRec, "CDELT1", dx / RADSINDEG);
00314         drms_setkey_float(outRec, "CDELT2", dy / RADSINDEG);
00315         drms_setkey_string(outRec, "CUNIT1", "degree");
00316         drms_setkey_string(outRec, "CUNIT2", "degree");
00317         
00318         char key[64];
00319         snprintf (key, 64, "CRLN-%s", wcsCode[proj]);
00320         drms_setkey_string(outRec, "CTYPE1", key);
00321         snprintf (key, 64, "CRLT-%s", wcsCode[proj]);
00322         drms_setkey_string(outRec, "CTYPE2", key);
00323         drms_setkey_float(outRec, "CROTA2", 0.0);
00324         
00325         printf(" done\n"); fflush(stdout);
00326         
00327     }       // irec
00328     
00329     drms_close_records(inRS, DRMS_FREE_RECORD);
00330     drms_close_records(outRS, DRMS_INSERT_RECORD);
00331     
00332     return DRMS_SUCCESS;
00333     
00334 }   // DoIt
00335
proj/cgem/prep/apps/cgem_map.c
