proj/util/apps/jsoc_rebin.c

#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);
  }

---