util/apps/jsoc_rebin.c

/**
   @defgroup jsoc_rebin jsoc_rebin reduce/increase image size by integer multiples
   @ingroup su_util

   @brief Reduce (increase) the dimension of the input data by an integer factor.

   @par Synopsis:
   @code
   jsoc_rebin  in=input_data out=output_data  scale=<scale> method={boxcar,gaussian}
   where scale is <1 for size reduction and >1 for size increase.  The output image scale
   will be the nearest integer to "scale" or its reciprocal for scale < 1.0. 
   @endcode

   This is a general purpose module that takes a series of input 
   data and modifies its spatial size by a factor 
   specified by "scale".
   The method for avaraging (interpolation) can be specified 
   through the input "method". The current version handles a simple 
   boxcar average and Gaussian filtered sampling. If 'scale' < 1 then the input is reduced in size.

   The image is not registered to solar center by this module.  The image will be rotated by
   a flip-flip procedure is the CROTA2 parameter is near +-180.0 unless the -u flag (unchanged) is
   present.  If the -c flag is present the image will be cropped at the solar limb before scaling.
   If the -h flag or requestid parameter is present the output segments will have full headers.

   If gaussian smoothing is specified (via method) then the FWHM and nvector parameters should also
   be provided.  These define the Gaussian smoothing vector.  nvector is the full width of the
   smoothing function.  nvector will be adjusted such that nvector is odd if 1/scale rounds to an odd integer.

   @par Flags:
   @c
   -A  Process all segments
   -c  Crop before scaling.  Use rsun_obs/cdelt1 for limb radius.
   -h  Write full FITS headers.
   -u  Leave unchanged, do NOT rotate by 180 degrees if CROTA2 is near 180.  default is to do flip-flip method so
       image is norths up and no pixel values are changed.

   @par GEN_FLAGS:
   Ubiquitous flags present in every module.
   @ref jsoc_main

   @param in  The input data series.
   @param out The output series.

   @par Exit_Status:

   @par Example:
   Takes a series of 1024 X 1024 MDI full disk magnetogram and
   produces images with the resolution of HMI, 4096 X 4096.

   @code
   jsoc_rebin in='mdi.fd_M_96m_lev18[2003.10.20/1d]' out='su_phil.mdi_M_4k' scale=4 method='boxcar'
   @endcode

   @par Example:
   Reduces the resolution of HMI images 4096 X 4096 to that of MDI,
   1024 X 1024. Here the input is the HMI images and the output
   is the lower resolution HMI images.  Crop and rotate before rescaling.
   @code
   jsoc_rebin  -c in=hmi.M_45s[2011.10.20/1d]' out='su_phil.hmi_M_1k_45s' scale=0.25 method='boxcar'
   @endcode

   @bug
   None known so far.

*/

#include "jsoc.h"
#include "jsoc_main.h"
#include "fstats.h"

char *module_name = "jsoc_rebin";

#define DIE(msg) {fflush(stdout);fprintf(stderr,"%s, status=%d\n",msg,status); return(status);}

ModuleArgs_t module_args[] =
{
     {ARG_STRING, "in", "NOT SPECIFIED",  "Input data series."},
     {ARG_STRING, "out", "NOT SPECIFIED",  "Output data series."},
     {ARG_FLAG, "A", "0", "Process all segments in the input series"},
     {ARG_FLAG, "c", "0", "Crop at rsun_obs."},
     {ARG_FLAG, "h", "0", "Include full FITS header in output segment."},
     {ARG_FLAG, "u", "0", "do not rotate by 180 if needed."},
     {ARG_FLOAT, "scale", "1.0", "Scale factor."},
     {ARG_FLOAT, "FWHM", "-1.0", "Smoothing Gaussian FWHM for method=gaussian."},
     {ARG_INT, "nvector", "-1.0", "Smoothing Gaussian vector length for method=gaussian."},
     {ARG_INT, "inseg", "-1", "Input segment number"},
     {ARG_INT, "outseg", "-1", "Output segment number"},
     {ARG_STRING, "method", "boxcar", "conversion type, one of: boxcar, gaussian."},
     {ARG_STRING, "requestid", "NA", "RequestID if called as an export processing step."},
     {ARG_END}
};

#define     Deg2Rad    (M_PI/180.0)
#define     Rad2arcsec (3600.0/Deg2Rad)
#define     arcsec2Rad (Deg2Rad/3600.0)
#define     Rad2Deg    (180.0/M_PI)

struct ObsInfo_struct
  {
  // from observation info
  TIME  t_obs;
  double rsun_obs, obs_vr, obs_vw, obs_vn;
  double crpix1, crpix2, cdelt1, cdelt2, crota2;
  double crval1, crval2;
  double cosa, sina;
  double obs_b0;
  // observed point
  int i,j;
  // parameters for observed point
  double x,y,r;
  double rho;
  double lon;
    double lat;
  double sinlat, coslat;
  double sig;
  double mu;
  double chi;
  double obs_v;
  };

typedef struct ObsInfo_struct ObsInfo_t;

void rebinArraySF(DRMS_Array_t *out, DRMS_Array_t *in);
int upNcenter(DRMS_Array_t *arr, ObsInfo_t *ObsLoc);
int crop_image(DRMS_Array_t *arr, ObsInfo_t *ObsLoc);
const char *get_input_recset(DRMS_Env_t *drms_env, const char *in);

ObsInfo_t *GetObsInfo(DRMS_Segment_t *seg, ObsInfo_t *pObsLoc, int *rstatus);
ObsInfo_t *GetMinObsInfo(DRMS_Segment_t *seg, ObsInfo_t *pObsLoc, int *rstatus);

int DoIt(void)
  {
  int status = DRMS_SUCCESS;
  DRMS_RecordSet_t *inRS, *outRS;
  int irec, nrecs;
  const char *inQuery = params_get_str(&cmdparams, "in");
  char *inStr;
  const char *outSeries = params_get_str(&cmdparams, "out");
  const char *method = params_get_str(&cmdparams, "method");
  const char *requestid = params_get_str(&cmdparams, "requestid");
  int nvector = params_get_int(&cmdparams, "nvector");
  float fscale = params_get_float(&cmdparams, "scale");
  float fwhm = params_get_float(&cmdparams, "FWHM");
  int crop = params_get_int(&cmdparams, "c");
  int as_is = params_get_int(&cmdparams, "u");
  int wantAllSegs =  params_get_int(&cmdparams, "A");
  int inseg = params_get_int(&cmdparams, "inseg");
  int outseg = params_get_int(&cmdparams, "outseg");
  int full_header = params_get_int(&cmdparams, "h") || strcmp(requestid, "NA");
  char *in_filename = NULL;
  DRMS_Record_t *template = NULL;
  int firstNx = 0, firstNy = 0, firstNaxis = 0;
  int wantSeg = 1;

  int iscale, ivec, vec_half;
  double *vector;
  char history[4096];

  if (fscale < 1.0) // shrinking
    iscale = 1.0/fscale + 0.5;
  else  // enlarging
    iscale = fscale + 0.5;
  if (nvector < 0)
    nvector = iscale;
  // Both 1/scale and nvector must be odd or both even so add 1 to nvector if needed       
  if (((iscale & 1) && !(nvector & 1)) || ((!(iscale & 1) && (nvector & 1) )))
    nvector += 1;
  vector = (double *)malloc(nvector * sizeof(double));
  vec_half = nvector/2; // counts on truncate to int if nvector is odd.

  if (strcasecmp(method, "boxcar")==0 && fscale < 1)
    {
    for (ivec = 0; ivec < nvector; ivec++)
      vector[ivec] = 1.0;
    sprintf(history, "Boxcar bin by %d%s%s",
      iscale,
      (crop ? ", Cropped at rsun_obs" : ""),
      (!as_is ? ", North is up" : "") );
    }
  else if (strcasecmp(method, "boxcar")==0 && fscale >= 1)
    {
    if (nvector != iscale)
      DIE("For fscale>=1 nvector must be fscale");
    for (ivec = 0; ivec < nvector; ivec++)
      vector[ivec] = 1.0;
    sprintf(history, "Replicate to expand by %d%s%s",
      iscale,
      (crop ? ", Cropped at rsun_obs" : ""),
      (!as_is ? ", North is up" : "") );
    }
  else if (strcasecmp(method, "gaussian")==0) // do 2-D vector weights calculated as Gaussian
    {
    if (fwhm < 0)
      DIE("Need FWHM parameter");
    for (ivec = 0; ivec < nvector; ivec++)
      {
      double arg = (ivec - (nvector-1)/2.0) * (ivec - (nvector-1)/2.0);
      vector[ivec] = exp(-arg/(fwhm*fwhm*0.52034));
      }
    sprintf(history, "Scale by %f with Gasussian smoothing FWHM=%f, nvector=%d%s%s",
      fscale, fwhm, nvector,
      (crop ? ", Cropped at rsun_obs" : ""),
      (!as_is ? ", North is up" : "") );
    }
  else
    DIE("invalid conversion method");

  printf(history);

  inStr = strdup(get_input_recset(drms_env, (char *)inQuery));
  if (!inStr || *inStr=='\0') DIE("Cant make special cadence recordset list file");
  inRS = drms_open_records(drms_env, inStr, &status);
  if (strcmp(inStr, inQuery) && *inStr == '@')
    unlink(inStr+1);
  if (status || inRS->n == 0)
    DIE("No input data found");

  nrecs = inRS->n;
  if (nrecs == 0)
    DIE("No records found");
  drms_stage_records(inRS, 1, 1);

  outRS = drms_create_records(drms_env, nrecs, (char *)outSeries, DRMS_PERMANENT, &status);
  if (status)
    DIE("Output recordset not created");
      
  int iSet;
  int iSubRec;
  int nSubRecs = 0;

  /* We need the record-set specfication for each record, so loop by recordset subset. */
  for (iSet = 0, irec = 0; iSet < inRS->ss_n; iSet++)
  {
      int hasSeglist = wantAllSegs ? 1 : 0;
      
      /* For each record, we need to know if the segments it contains came from a record-set specification that had a segment list.
       * To do that, we have to search for '{' in the specification that resolved to this record. So we need the spec for this record,
       * which is in inRS->ss_queries[iSet]. */
      if (strchr(inRS->ss_queries[iSet], '{'))
      {
          hasSeglist = 1;
      }
      
      nSubRecs = drms_recordset_getssnrecs(inRS, iSet, &status);
      
      for (iSubRec = 0; iSubRec < nSubRecs; iSubRec++, irec++)
      {
          ObsInfo_t *ObsLoc;
          DRMS_Record_t *outRec, *inRec;
          DRMS_Segment_t *outSeg, *inSeg;
          DRMS_Segment_t *orig = NULL;
          DRMS_Array_t *inArray, *outArray;
          float *inData, *outData;
          float val;
          int quality=0;
          
          inRec = inRS->records[(inRS->ss_starts)[iSet] + iSubRec];
          
          quality = drms_getkey_int(inRec, "QUALITY", &status);
          if (status || (!status && quality >= 0))
          {
              /* Segment loop. */
              HIterator_t *segIter = NULL;
              
              while ((inSeg = drms_record_nextseg2(inRec, &segIter, 1, &orig)) != NULL)
              {
                  /* A segment can only be processed if it is dimension 2 and if not first segment has same dimension
                   * as first segment.
                   */
                  if (inseg >= 0)
                  {
                      /* This module was called with the inseg argument (which is the input series' segment number).
                       * Process only that segment. But inSeg could be a linked segment, so we have to obtain the
                       * segment struct in the original series. There is no pointer from the linked segment to
                       * the original segment. And we do not have the name or number of the original segment. Or we
                       * didn't, so I made drms_record_nextseg2(). */
                      if (orig->info->segnum != inseg)
                      {
                          /* Skip this segment. It isn't the one the caller has requested. */
                          continue;
                      }
                  }
                  wantSeg = 1; // start by assuming OK
                  
                  inArray = drms_segment_read(inSeg, DRMS_TYPE_FLOAT, &status);
                  if (status)
                      {
                      printf(" No data file found but QUALITY not bad, status=%d\n", status);
                      drms_free_array(inArray);
                      continue;
                      }
                  int inx, iny, outx, outy, i, j;
                  int in_nx = inArray->axis[0];
                  int in_ny = inArray->axis[1];
                  int out_nx = in_nx * fscale + 0.5;
                  int out_ny = in_ny * fscale + 0.5;
                  int outDims[2] = {out_nx, out_ny};
                  if (firstNaxis == 0)
                      {
                      if (inArray->naxis != 2)
                          wantSeg = 0; // Skip segment if not 2 axes. mark for skip below
                      else
                          {
                          firstNaxis = inArray->naxis;
                          firstNx = in_nx;
                          firstNy = in_ny;
                          }
                      }
                  if (inArray->naxis != firstNaxis || in_nx != firstNx || in_ny != firstNy)
                      {  // this segment now wanted, if no segment list provided, copy this segment as-is
                      wantSeg = 0;
                      }
 
                  if (wantSeg && (crop || !as_is))
                      {
                      ObsLoc = GetObsInfo(inSeg, NULL, &status);
                      if (!as_is) upNcenter(inArray, ObsLoc);
                      if (crop) crop_image(inArray, ObsLoc);
                      }
                  else if (wantSeg)
                      {
                      ObsLoc = GetMinObsInfo(inSeg, NULL, &status);
                      }
                  
                  if (wantSeg)
                      {
                      inData = (float *)inArray->data;
                      outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, outDims, NULL, &status);
                      outData = (float *)outArray->data;
                      }
                  
                  if (wantSeg && fscale > 1.0)
                  {
                      int out_go = (iscale-1)/2.0 + 0.5;
                      for (iny = 0; iny < in_ny; iny += 1)
                          for (inx = 0; inx < in_nx; inx += 1)
                          {
                              val = inData[in_nx*iny + inx];
                              for (j = 0; j < nvector; j += 1)
                              {
                                  outy = iny*iscale + out_go + j - vec_half;
                                  for (i = 0; i < nvector; i += 1)
                                  {
                                      outx = inx*iscale + out_go + i - vec_half;
                                      if (outx >= 0 && outx < out_nx && outy >= 0 && outy < out_ny)
                                          outData[out_nx*outy + outx] = val;
                                  }
                              }
                          }
                  }
                  else if (wantSeg)
                  {
                      int in_go = (iscale-1)/2.0 + 0.5;
                      for (outy = 0; outy < out_ny; outy += 1)
                          for (outx = 0; outx < out_nx; outx += 1)
                          {
                              double total = 0.0;
                              double weight = 0.0;
                              int nn = 0;
                              for (j = 0; j < nvector; j += 1)
                              {
                                  iny = outy*iscale + in_go + j - vec_half;
                                  for (i = 0; i < nvector; i += 1)
                                  {
                                      inx = outx*iscale + in_go + i - vec_half;
                                      if (inx >= 0 && inx < in_nx && iny >=0 && iny < in_ny)
                                      {
                                          val = inData[in_nx*(iny) + inx];
                                          if (!drms_ismissing_float(val))
                                          {
                                              double w = vector[i]*vector[j];
                                              total += w*val;
                                              weight += w;
                                              nn++;
                                          }
                                      }
                                  }
                              }
                              outData[out_nx*outy + outx] = (nn > 0 ? total/weight : DRMS_MISSING_FLOAT);
                          }
                  }
                  
                  
                  // write data file
                  outRec = outRS->records[irec];
                  
                  /* The original intent was for the output series' segments to parallel the input series' segments.
                   * However, the names were allowed to vary. So the segment descriptions match, and they are in the
                   * same order. However, the names given to these descriptions might vary. Therefore, we need to
                   * look-up output segments by segment number, not name.
                   *
                   * If outseg was provided in the arguments, use that. Otherwise, use the segment number of the input
                   * segment (the original input segment, in case a link was followed).
                   */
                  outSeg = drms_segment_lookupnum(outRec, (outseg >= 0) ? outseg : orig->info->segnum);
                  
                  if (!outSeg)
                  {
                      DIE("Unable to look-up output segment.");
                  }
                  
                  // copy all keywords
                  drms_copykeys(outRec, inRec, 0, kDRMS_KeyClass_Explicit);
                  
                  // Now fixup coordinate keywords
                  // Only CRPIX1,2 and CDELT1,2 and CROTA2 should need repair.
                  drms_setkey_double(outRec, "CDELT1", ObsLoc->cdelt1/fscale);
                  drms_setkey_double(outRec, "CDELT2", ObsLoc->cdelt2/fscale);
                  drms_setkey_double(outRec, "CRPIX1", (ObsLoc->crpix1-0.5) * fscale + 0.5);
                  drms_setkey_double(outRec, "CRPIX2", (ObsLoc->crpix2-0.5) * fscale + 0.5);
                  drms_setkey_double(outRec, "CROTA2", ObsLoc->crota2);
                  
                  if (strcasecmp(method, "gaussian")==0)
                      drms_setkey_double(outRec, "FWHM", fwhm);
                  drms_appendhistory(outRec, history, 1);
                  drms_setkey_time(outRec, "DATE", CURRENT_SYSTEM_TIME);
                  if (strcmp(requestid, "NA") != 0)
                      drms_setkey_string(outRec, "RequestID", requestid);
                  

                  if (!wantSeg)
                      {
                      drms_free_array(inArray);
                      inArray = drms_segment_read(inSeg,DRMS_TYPE_RAW,&status);
                      if (!inArray || status)
                          DIE("Problem reading ancillary array");
                      status = drms_segment_write(outSeg, inArray, 0);
                      if (status)
                          DIE("Problem writing ancillary array");
                      drms_free_array(inArray);
                      }
                  else
                      {
                      // get info for array from input segment
                      outArray->parent_segment = outSeg;
                      
                      drms_setkey_int(outRec, "TOTVALS", out_nx*out_ny);
                      set_statistics(outSeg, outArray, 1);
                      
                      // Use the input array as the best guess for scale and zero
                      outArray->bzero = inArray->bzero;
                      outArray->bscale = inArray->bscale;

                      drms_free_array(inArray);
                      if (full_header)
                          status = drms_segment_writewithkeys(outSeg, outArray, 0);
                      else
                          status = drms_segment_write(outSeg, outArray, 0);
                      if (status)
                          DIE("problem writing file");
                      drms_free_array(outArray);
                      }
                  
                  if (inseg >= 0)
                  {
                      /* The caller requested the processing of a single segment. Exit loop. */
                      break;
                  }
                  
                  /* If the user did not provide a segment specifier (which means that all segment structs are in the rec->segments container,
                   * and there was no seglist), then process only the first segment. */
                  if (inRec->seriesinfo && *inRec->seriesinfo->seriesname)
                  {
                     template = drms_template_record(drms_env, inRec->seriesinfo->seriesname, &status);
                  }
                  else
                  {
                     template = NULL;
                  }

                  if (status || !template)
                  {
                      DIE("Cannot obtain DRMS-record template.");
                  }

                  if (drms_record_numsegments(inRec) == drms_record_numsegments(template) && !hasSeglist)
                  {
                      break;
                  }
              } // end of segment loop
              
              if (segIter)
              {
                  hiter_destroy(&segIter);
              }
          } // record quality check
      }
  } //end of "irec" loop

  drms_close_records(inRS, DRMS_FREE_RECORD);
  drms_close_records(outRS, DRMS_INSERT_RECORD);
  return (DRMS_SUCCESS);
  } // end of DoIt

// ----------------------------------------------------------------------

/* center whith whole pixel shifts and rotate by 180 if needed */
/* Only apply center if it will not result in an image crop.  I.e. not ever
   for AIA, and not for HMI or MDI or other if a shift of more than 20 arcsec
   is implied  */
int upNcenter(DRMS_Array_t *arr, ObsInfo_t *ObsLoc)
  {
  int nx, ny, ix, iy, i, j, xoff, yoff, max_off;
  double rot, x0, y0, midx, midy;
  float *data;
  float *data2;
  if (!arr || !ObsLoc)
    return(1);
  data = arr->data;
  nx = arr->axis[0];
  ny = arr->axis[1];
  x0 = ObsLoc->crpix1 - 1;
  y0 = ObsLoc->crpix2 - 1;
  midx = (nx-1.0)/2.0;
  midy = (ny-1.0)/2.0;
  if ((rot = fabs(ObsLoc->crota2)) > 179 && rot < 181)
    {
    // rotate image by 180 degrees by a flip flip
    float val;
    int half = ny / 2;
    int odd = ny & 1;
    for (iy=0; iy<half; iy++)
      {
      for (ix=0; ix<nx; ix++)
        {
        i = iy*nx + ix;
        j = (ny - 1 - iy)*nx + (nx - 1 - ix);
        val = data[i];
        data[i] = data[j];
        data[j] = val;
        }
      }
    // reverse middle row if ny is odd
    if (odd) {
        for (ix=0; ix<nx/2; ++ix) {
            i = nx*half + ix;
            j = nx*half + nx - ix - 1;
            val = data[i];
            data[i] = data[j];
            data[j] = val;
        }
    }
    x0 = nx - 1 - x0;
    y0 = ny - 1 - y0;
    rot = ObsLoc->crota2 - 180.0;
    if (rot < -90.0) rot += 360.0;
    ObsLoc->crota2 = rot;
    }
  // Center to nearest pixel - if OK to do so
  xoff = round(x0 - midx);
  yoff = round(y0 - midy);
  max_off = 20.0 / ObsLoc->cdelt1;
  if (arr->parent_segment &&
      arr->parent_segment->record &&
      arr->parent_segment->record->seriesinfo && 
      arr->parent_segment->record->seriesinfo->seriesname && 
      strncasecmp(arr->parent_segment->record->seriesinfo->seriesname, "aia", 3) &&
      abs(xoff) < max_off && abs(yoff) < max_off) 
    {
    if (abs(xoff) >= 1 || abs(yoff) >= 1)
      {
      data2 = malloc(4*nx*ny);
      for (iy=0; iy<ny; iy++)
        {
        int jy = iy + yoff;
        for (ix=0; ix<nx; ix++)
          {
          int jx = ix + xoff;
          int idx = jy*nx + jx;
          int idx2 = iy*nx + ix;
          if (jx<0 || jx>=nx || jy<0 || jy>=ny)
            data2[idx2] = DRMS_MISSING_FLOAT;
          else
            data2[idx2] = data[idx];
          }
        }
      x0 -= xoff;
      y0 -= yoff;
      free(data);
      arr->data = data2;
      }
    }
  // update center location
  ObsLoc->crpix1 = x0 + 1;
  ObsLoc->crpix2 = y0 + 1;
  return(0);
  }

// ----------------------------------------------------------------------

int crop_image(DRMS_Array_t *arr, ObsInfo_t *ObsLoc)
  {
  int nx, ny, ix, iy, i, j, xoff, yoff;
  double x0, y0;
  double rsun = ObsLoc->rsun_obs/ObsLoc->cdelt1;
  double scale, crop_limit2;
  float *data;
  if (!arr || !ObsLoc)
    return(1);
  data = arr->data;
  nx = arr->axis[0];
  ny = arr->axis[1];
  x0 = ObsLoc->crpix1 - 1;
  y0 = ObsLoc->crpix2 - 1;
  scale = 1.0/rsun;
  // crop_limit = 0.99975; // 1 - 1/4000, 1/2 HMI pixel.
  crop_limit2 = 0.99950; // square of 1 - 1/4000, 1/2 HMI pixel.
  for (iy=0; iy<ny; iy++)
    for (ix=0; ix<nx; ix++)
      {
      double x, y, R2;
      float *Ip = data + iy*nx + ix;
      if (drms_ismissing_float(*Ip))
        continue;
      x = ((double)ix - x0) * scale; /* x,y in pixel coords */
      y = ((double)iy - y0) * scale;
      R2 = x*x + y*y;
      if (R2 > crop_limit2)
        *Ip = DRMS_MISSING_FLOAT;
      }
  return(0);
  }

// ----------------------------------------------------------------------

#define CHECK(keyname) {if (status) {fprintf(stderr,"Keyword failure to find: %s, status=%d\n",keyname,status); *rstatus=status; return(NULL);}}

ObsInfo_t *GetObsInfo(DRMS_Segment_t *seg, ObsInfo_t *pObsLoc, int *rstatus)
  {
  TIME t_prev;
  DRMS_Record_t *rec;
  TIME t_obs;
  double dv;
  ObsInfo_t *ObsLoc;
  int status;

  if (!seg || !(rec = seg->record))
    { *rstatus = 1; return(NULL); }

  ObsLoc = (pObsLoc ? pObsLoc : (ObsInfo_t *)malloc(sizeof(ObsInfo_t)));
  if (!pObsLoc)
    memset(ObsLoc, 0, sizeof(ObsInfo_t));

  t_prev = ObsLoc->t_obs;
  t_obs = drms_getkey_time(rec, "T_OBS", &status); CHECK("T_OBS");

  if (t_obs <= 0.0)
    { *rstatus = 2; return(NULL); }

  if (t_obs != t_prev)
    {
    ObsLoc->crpix1 = drms_getkey_double(rec, "CRPIX1", &status); CHECK("CRPIX1");
    ObsLoc->crpix2 = drms_getkey_double(rec, "CRPIX2", &status); CHECK("CRPIX2");
    ObsLoc->crval1 = drms_getkey_double(rec, "CRVAL1", &status); CHECK("CRVAL1");
    ObsLoc->crval2 = drms_getkey_double(rec, "CRVAL2", &status); CHECK("CRVAL2");
    ObsLoc->cdelt1 = drms_getkey_double(rec, "CDELT1", &status); CHECK("CDELT1");
    ObsLoc->cdelt2 = drms_getkey_double(rec, "CDELT2", &status); CHECK("CDELT1");
    ObsLoc->crota2 = drms_getkey_double(rec, "CROTA2", &status); if (status) ObsLoc->crota2 = 0.0; // WCS default
    ObsLoc->sina = sin(ObsLoc->crota2*Deg2Rad);
    ObsLoc->cosa = sqrt (1.0 - ObsLoc->sina*ObsLoc->sina);
    ObsLoc->rsun_obs = drms_getkey_double(rec, "RSUN_OBS", &status);
    if (status)
      {
      double dsun_obs = drms_getkey_double(rec, "DSUN_OBS", &status); CHECK("DSUN_OBS");
      ObsLoc->rsun_obs = asin(696000000.0/dsun_obs)/arcsec2Rad;
      }
    ObsLoc->obs_vr = drms_getkey_double(rec, "OBS_VR", &status); CHECK("OBS_VR");
    ObsLoc->obs_vw = drms_getkey_double(rec, "OBS_VW", &status); CHECK("OBS_VW");
    ObsLoc->obs_vn = drms_getkey_double(rec, "OBS_VN", &status); CHECK("OBS_VN");
    ObsLoc->obs_b0 = drms_getkey_double(rec, "CRLT_OBS", &status); CHECK("CRLT_OBS");
    ObsLoc->t_obs = t_obs;
    }
  *rstatus = 0;
  return(ObsLoc);
  }

/* GetMinObsInfo - gets minimum standard WCS keywords for e.g. heliographic mapped data */
ObsInfo_t *GetMinObsInfo(DRMS_Segment_t *seg, ObsInfo_t *pObsLoc, int *rstatus)
  { 
  TIME t_prev;
  DRMS_Record_t *rec;
  TIME t_obs;
  double dv;
  ObsInfo_t *ObsLoc;
  int status;

  if (!seg || !(rec = seg->record))
    { *rstatus = 1; return(NULL); }

  ObsLoc = (pObsLoc ? pObsLoc : (ObsInfo_t *)malloc(sizeof(ObsInfo_t)));
  if (!pObsLoc)
    memset(ObsLoc, 0, sizeof(ObsInfo_t));

  t_prev = ObsLoc->t_obs;
  t_obs = drms_getkey_time(rec, "T_OBS", &status); CHECK("T_OBS");

  if (t_obs <= 0.0)
    { *rstatus = 2; return(NULL); }

  if (t_obs != t_prev)
    {
    ObsLoc->crpix1 = drms_getkey_double(rec, "CRPIX1", &status); CHECK("CRPIX1");
    ObsLoc->crpix2 = drms_getkey_double(rec, "CRPIX2", &status); CHECK("CRPIX2");
    ObsLoc->crval1 = drms_getkey_double(rec, "CRVAL1", &status); CHECK("CRVAL1");
    ObsLoc->crval2 = drms_getkey_double(rec, "CRVAL2", &status); CHECK("CRVAL2");
    ObsLoc->cdelt1 = drms_getkey_double(rec, "CDELT1", &status); CHECK("CDELT1");
    ObsLoc->cdelt2 = drms_getkey_double(rec, "CDELT2", &status); CHECK("CDELT1");
    ObsLoc->crota2 = drms_getkey_double(rec, "CROTA2", &status); if (status) ObsLoc->crota2 = 0.0; // WCS default
    ObsLoc->sina = sin(ObsLoc->crota2*Deg2Rad);
    ObsLoc->cosa = sqrt (1.0 - ObsLoc->sina*ObsLoc->sina);
    ObsLoc->rsun_obs = drms_getkey_double(rec, "RSUN_OBS", &status);
    ObsLoc->obs_vr = drms_getkey_double(rec, "OBS_VR", &status);
    ObsLoc->obs_vw = drms_getkey_double(rec, "OBS_VW", &status);
    ObsLoc->obs_vn = drms_getkey_double(rec, "OBS_VN", &status);
    ObsLoc->obs_b0 = drms_getkey_double(rec, "CRLT_OBS", &status);
    ObsLoc->t_obs = t_obs;
    }
  *rstatus = 0;
  return(ObsLoc);
  }

// ----------------------------------------------------------------------

// In cases known to not have compact slotted series and cadence is specified
// generate explicit recordset list of closest good record to desired grid
// First get vector of times and quality
// Then if vector is not OK, quit.
// then: make temp file to hold recordset list
//       start with first time to define desired grid,
//       make array of desired times.
//       make empty array of recnums
//       search vector for good images nearest desired times
//       for each found time, write record query


const char *get_input_recset(DRMS_Env_t *drms_env, const char *inQuery)
  {
  static char newInQuery[102];
  TIME epoch = (cmdparams_exists(&cmdparams, "epoch")) ? params_get_time(&cmdparams, "epoch") : 0;
  DRMS_Array_t *data;
  TIME t_start, t_stop, t_now, t_want, t_diff, this_t_diff;
  int status = 1;
  int nrecs, irec;
  int nslots, islot;
  long long *recnums;
  TIME *t_this, half;
  TIME cadence;
  double *drecnum, *dquality;
  int quality;
  long long recnum;
  char keylist[DRMS_MAXQUERYLEN];
  static char filename[100];
  char *tmpdir;
  FILE *tmpfile;
  char newIn[DRMS_MAXQUERYLEN];
  char seriesname[DRMS_MAXQUERYLEN];
  char *lbracket;
  char *at = index(inQuery, '@');
      
  if (at && *at && (strncmp(inQuery,"aia.lev1[", 9)==0 ||
                    strncmp(inQuery,"hmi.lev1[", 9)==0 ||
                    strncmp(inQuery,"aia.lev1_nrt2[",14)==0 ||
                    strncmp(inQuery,"hmi.lev1_nrt[", 13)==0 ))
    {
    char *ip=(char *)inQuery, *op=newIn, *p;
    long n, mul;
        
    char *segSpec = NULL;
    const char *psl = NULL;
    char *pbracket = NULL;
        
    psl = strchr(inQuery, '{');
    if (psl)
    {
        segSpec = strdup(psl + 1);
        
        if (!segSpec)
        {
            fprintf(stderr, "No memory.\n");
            return NULL;
        }
        
        pbracket = strchr(segSpec, '}');
        
        if (!pbracket)
        {
            fprintf(stderr, "Invalid segment specification.\n");
            return NULL;
        }
        
        *pbracket = '\0';
        
        if (!*segSpec)
        {
            fprintf(stderr, "Invalid segment specification.\n");
            return NULL;
        }
    }
        
    while ( *ip && ip<at )
      *op++ = *ip++;
    ip++; // skip the '@'
    n = strtol(ip, &p, 10); // get digits only
    if (*p == 's') mul = 1;
    else if (*p == 'm') mul = 60;
    else if (*p == 'h') mul = 3600;
    else if (*p == 'd') mul = 86400;
    else 
      {
      fprintf(stderr,"cant make sense of @xx cadence spec for aia or hmi lev1 data");
      return(NULL);
      }
    cadence = n * mul;
    ip = ++p;  // skip cadence multiplier
    while ( *ip )
      *op++ = *ip++;
    *op = '\0';
    half = cadence/2.0;
    sprintf(keylist, "T_OBS,QUALITY,recnum");
    data = drms_record_getvector(drms_env, newIn, keylist, DRMS_TYPE_DOUBLE, 0, &status);
    if (!data || status)
      {
      fprintf(stderr,"getkey_vector failed status=%d\n", status);
      return(NULL);
      }
    nrecs = data->axis[1];
    irec = 0;
    t_this = (TIME *)data->data;
    dquality = (double *)data->data + 1*nrecs;
    drecnum = (double *)data->data + 2*nrecs;
    if (epoch > 0.0)
      {
      int s0 = (t_this[0] - epoch)/cadence;
      TIME t0 = s0*cadence + epoch;
      t_start = (t0 < t_this[0] ? t0 + cadence : t0);
      }
    else
      t_start = t_this[0];
    t_stop = t_this[nrecs-1];
    nslots = (t_stop - t_start + cadence/2)/cadence;
    recnums = (long long *)malloc(nslots*sizeof(long long));
    for (islot=0; islot<nslots; islot++)
      recnums[islot] = 0;
    islot = 0;
    t_want = t_start;
    t_diff = 1.0e9;
    for (irec = 0; irec<nrecs; irec++)
        {
        t_now = t_this[irec];
        quality = (int)dquality[irec] & 0xFFFFFFFF;
        recnum = (long long)drecnum[irec];
        this_t_diff = fabs(t_now - t_want);
        if (quality < 0)
          continue;
        if (t_now <= (t_want-half))
          continue;
        while (t_now > (t_want+half))
          {
          islot++;
          if (islot >= nslots)
             break;
          t_want = t_start + cadence * islot;
          this_t_diff = fabs(t_now - t_want);
          t_diff = 1.0e8;
          }
        if (islot < nslots && this_t_diff <= t_diff)
          recnums[islot] = recnum;
        t_diff = fabs(t_now - t_want);
        }
    if (islot+1 < nslots)
      nslots = islot+1;  // take what we got.
    strcpy(seriesname, inQuery);
    lbracket = index(seriesname,'[');
    if (lbracket) *lbracket = '\0';
    
    tmpdir = getenv("TMPDIR");
    if (!tmpdir) tmpdir = "/tmp";
    sprintf(filename, "%s/hg_patchXXXXXX", tmpdir);
    mkstemp(filename);
    tmpfile = fopen(filename,"w");
    for (islot=0; islot<nslots; islot++)
      if (recnums[islot])
      {
          if (!segSpec)
          {
              fprintf(tmpfile, "%s[:#%lld]\n", seriesname, recnums[islot]);
          }
          else
          {
              fprintf(tmpfile, "%s[:#%lld]{%s}\n", seriesname, recnums[islot], segSpec);
          }
      }
        
    if (segSpec)
        {
        free(segSpec);
        }
        
    fclose(tmpfile);
    free(recnums);
    drms_free_array(data);
    sprintf(newInQuery,"@%s", filename);
    return(newInQuery);
    }
  else
    return(inQuery);
  }
Revision:	1.21
Committed:	Mon Jul 17 16:45:23 2017 UTC (6 years, 2 months ago) by kehcheng
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Ver_9-1, Ver_LATEST, Ver_9-3, Ver_9-41, Ver_9-2, Ver_9-5, Ver_9-4
Changes since 1.20:	+10 -1 lines
Log Message:	fix upNcenter() for images with odd number of rows