/**
 * \file xy2ncvector.c
 *
 * $Id$
 *
 * Convert .xy files into NetCDF containing vectors with pairs
 * of coordinates to describe lines,
 * suitable to plot in ferret with plot/vs/line
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <netcdf.h>

#define ncerror(e) { if(e) {fprintf(stderr,"Error: %s\n", nc_strerror(e)); exit(1);}}

static char buf[1000];
static int dbrecnr = 0;

static
void usage(char *myname, char *err) {
  if (err) fprintf(stderr, "%s\n", err);
  fprintf(stderr, "usage: %s infilename outbasename\n", myname);
  exit(1);
}

static int ncid, linesxid, linesyid;
static size_t linesindex[2] = { 0, 0 };

static void create_cdf(const char *filename, const char *myname)
{
  float missfloat = NC_FILL_FLOAT;
  int rc,dim[2];
  char s[1000];
  time_t t;
  int linelen_dim_id, lines_dim_id; // record dimension is number of lines

  rc=nc_create(filename,NC_CLOBBER,&ncid);
  if(rc) {
    fprintf(stderr,"ERROR426: Cannot create file '%s': %s.\n",
            filename,nc_strerror(rc));
    return;
  }

  rc=nc_def_dim(ncid,"routelen",3,&linelen_dim_id);
  ncerror(rc);
  rc=nc_def_dim(ncid,"lines",NC_UNLIMITED,&lines_dim_id);
  ncerror(rc);

  snprintf(s,999,"%s %s",myname, filename);
  rc=nc_put_att_text(ncid,NC_GLOBAL,"source",strlen(s),s);
  ncerror(rc);
  time(&t);
  //snprintf(s,999,"Created for user %s on %s at %s",getuser(),gethost(),
  //       ctime(&t));
  //s[strlen(s)-1]='\0';
  //rc=nc_put_att_text(ncid,NC_GLOBAL,"history",strlen(s),s);
  //ncerror(rc);


  dim[1] = linelen_dim_id;
  dim[0] = lines_dim_id; // record dimension (if used) must be first
  rc=nc_def_var(ncid,"linesx",NC_FLOAT,2,dim,&linesxid);
  ncerror(rc);
  rc=nc_put_att_text(ncid, linesxid,"units",strlen("degrees_east"),"degrees_east");
  ncerror(rc);
  nc_put_att_float(ncid, linesxid,"missing_value",NC_FLOAT,1,&missfloat);
  rc=nc_put_att_float(ncid, linesxid,"_FillValue",NC_FLOAT,1,&missfloat);
  ncerror(rc);
  rc=nc_put_att_text(ncid, linesxid,"ferret_command", strlen("plot/vs/line/over linesx,linesy"), "plot/vs/line/over routex,linesy");
  ncerror(rc);

  rc=nc_def_var(ncid,"linesy",NC_FLOAT,2,dim,&linesyid);
  ncerror(rc);
  rc=nc_put_att_text(ncid, linesyid,"units",strlen("degrees_north"),"degrees_north");
  ncerror(rc);
  nc_put_att_float(ncid, linesyid,"missing_value",NC_FLOAT,1,&missfloat);
  rc=nc_put_att_float(ncid, linesyid,"_FillValue",NC_FLOAT,1,&missfloat);
  ncerror(rc);

  rc=nc_enddef(ncid);
  ncerror(rc);
  return;
}


static int readpair(FILE *infile, float *x, float *y) {
    int err;
    if ('>' == buf[0]) return 0;
    err = sscanf(buf, "%g %g", x, y);
    if (err != 2) {
      fprintf(stderr, "sscanf did not match 2 floats\n");
      return 0;
    }
    //if (x < -180) fprintf(stderr, "\nfound small X %g\n", *x);
    // avoid wrap over date line
    if (*x < 0) *x += 360;
    if (*x > 360) *x -= 360;
    if (*x < 0) *x += 360;
    fgets(buf, sizeof(buf), infile);
    return 1;
}

static
void readonepoly(FILE *infile) {
  /* upon entering, the first data line is already in buf */
  float x[3], y[3];
  int foundone;
  int rc;

  fprintf(stderr, "starting polygon nr %d ", dbrecnr);
  x[2] = y[2] = NC_FILL_FLOAT;
  readpair(infile, x+0, y+0);
  do {
    foundone = readpair(infile, x+1, y+1);
    if (foundone) {
      // write out one pair of points == two pairs of coordinates for one line
      linesindex[1] = 0;
      rc = nc_put_var1_float(ncid, linesxid, linesindex, x+0);
      ncerror(rc);
      rc = nc_put_var1_float(ncid, linesyid, linesindex, y+0);
      ncerror(rc);
      linesindex[1] = 1;
      rc = nc_put_var1_float(ncid, linesxid, linesindex, x+1);
      ncerror(rc);
      rc = nc_put_var1_float(ncid, linesyid, linesindex, y+1);
      ncerror(rc);
      // leave 3rd entry as missing_value
      //linesindex[1] = 2;
      //rc = nc_put_var1_float(ncid, linesxid, linesindex, &missfloat);
      //error(rc);
      //rc = nc_put_var1_float(ncid, linesyid, linesindex, &missfloat);
      //error(rc);
      linesindex[0]++;
      x[0] = x[1];
      y[0] = y[1];
    } else return;
  } while (!feof(infile));
  dbrecnr++;
}

int main(int argc, char **argv) {
  FILE *infile;
  char *p;
  char *myname = argv[0];

  if (argc != 3) { usage(myname, "need input and output file names"); }
  infile = fopen(argv[1], "r");
  if (!infile) { perror(argv[1]); exit(1); }
  create_cdf(argv[2], myname);

  while (!feof(infile)) {
    p = fgets(buf, sizeof(buf), infile);
    if (!p) break;
    if ('>' == *p) {
      //fprintf(stderr, "skipping %s\n", buf);
      continue;
    }
    readonepoly(infile);
  }
  fprintf(stderr, "found %d polygons\n", dbrecnr);
  nc_close(ncid);
  return 0;
}