/*
 *  GVDR library for reading GVDR data files
 *  Copyright (C) 1994 Michael J. Maurer
 *
 *  This library is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU Library General Public
 *  License as published by the Free Software Foundation; either
 *  version 2 of the License, or (at your option) any later version.
 *
 *  This library is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *  Library General Public License for more details.
 *
 *  You should have received a copy of the GNU Library General Public
 *  License along with this library; if not, write to the Free
 *  Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 *
 *  Michael Maurer <maurer@nova.stanford.edu>
 *  Durand Bldg - Room 232
 *  Stanford, CA 94305-4055
 *  (415) 723-1024
 */
static char rcsid[]="$Id: gvp2mat.c,v 1.2 1994/05/13 01:51:39 maurer Exp $";
/******************************************************************************
  gvp2mat.c

     Function: Extract values from PDS-format GVDR files and write to
  Matlab format MAT files.

This file is part of the STARLab Magellan Altimeter Data Processing Software.
Michael Maurer, March 1994.
******************************************************************************/
/* $Log: gvp2mat.c,v $
 * Revision 1.2  1994/05/13  01:51:39  maurer
 * Replaced calls to fopen() with calls to gvfopen().
 *
 * Revision 1.1  1994/05/11  02:55:18  maurer
 * If'ed out file reads with new match variable.
 *
 * Revision 1.0  1994/04/29  02:04:54  maurer
 * Initial revision
 * */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <limits.h>
#include "libmisc.h"
#include "gvdr.h"
#include "pdstab.h"
#include "gvpdstab.h"	/* for pt_gvhdr and gvpdslib.p */
#include "gcvt.h"
#include "proj.h"	/* for gv_proj.p */
#include "tgvdrmat.h"

#define GVPDS2MAT_C
#include "gvp2mat.p"
#include "gvpdslib.p"
#include "pdstab.p"
#include "gvdr.p"
#include "gcvt.p"
#include "gv_proj.p"
#include "tgvdrmat.p"

static FILE	*fout	= NULL;
static int	bbxlo	= INT_MIN;
static int	bbxhi	= INT_MAX;
static int	bbylo	= INT_MIN;
static int	bbyhi	= INT_MAX;
static int	geog	= 0;				/* cartesian/geographic */
static char	*regexp	= ".*";
static int	nre_lo	= 0;
static int	nre_hi	= 0;
#define NRE 100
static re_limit_t	re_lo[NRE];
static re_limit_t	re_hi[NRE];

int
main(int argc,
     char *argv[])
{
   init(argc,argv);
   gvpds2mat();
   return 0;
}

static void
gvpds2mat(void)
{
   FILE *fp;
   char *p;
   gvhdr_t gvh;
   gvdrhdr_t hdr;
   gvdrcell_t gv;
   gpds_tidx_t *tidx = NULL;
   gpds_pidx_t *pidx = NULL;
   gvdrxif_rec *xif = NULL;
   gvdredf_rec *edf = NULL;
   gvdradf_rec *adf = NULL;
   void *anf = NULL;
   proj_t P;
   int err,npx,npy,nx,ny,txlo,tylo,pxlo,pylo,txhi,tyhi,pxhi,pyhi;
   int ty,ylo,yhi,tx,xlo,xhi,py,px,xr,yr,x,y;
   gvmatch_t match;

   printf("reading header in GVHDR\n");
   if ((fp=gvfopen(p="GVHDR.TAB","r"))==NULL)
      error(-1,errno,"[gvpds2mat] gvfopen %s failed",p);
   if (gvhdr_read(fp,&gvh))
      error(-1,0,"[gvpds2mat] gvhdr_read failed");
   fclose(fp);
   gcvt_setmode(GCVT_READ,gvh.gvh_fmt_float,gvh.gvh_fmt_byte);
   gvdrhdr_fill(&gvh,&hdr);
   if (gvdr_rproj(&hdr.hdr3,&P))
      error(-1,0,"[gvpds2mat] gvdr_rproj failed");

   get_bbox(&hdr);
   nx=bbxhi-bbxlo+1;
   ny=bbyhi-bbylo+1;
   npx=hdr.hdr3.gh_tile_npx;
   npy=hdr.hdr3.gh_tile_npy;
   tgvdrmat_init(&hdr.hdr3,nx,ny,regexp,re_lo,re_hi,nre_lo,nre_hi,&match);
   if (!match.any)
      error(-1,0,"no GVDR variables match pattern %s",regexp);

   gvpds_p2t(&gvh,bbxlo,bbyhi,&txlo,&tylo,&pxlo,&pylo);
   gvpds_p2t(&gvh,bbxhi,bbylo,&txhi,&tyhi,&pxhi,&pyhi);
   if (txhi<txlo || tyhi<tylo ||
       (txhi==txlo && pxhi<pxlo) ||
       (tyhi==tylo && pyhi<pylo))
      error(-1,0,"[gvpds2mat] invalid bounding box");
   if (0) {
      printf("bbxlo=%d bbxhi=%d   bbylo=%d bbyhi=%d\n",
	     bbxlo,bbxhi,bbylo,bbyhi);
      printf("txlo=%4d pxlo=%4d tylo=%4d pylo=%4d\n",
	     txlo,pxlo,tylo,pylo);
      printf("txhi=%4d pxhi=%4d tyhi=%4d pyhi=%4d\n",
	     txhi,pxhi,tyhi,pyhi);
   }

   printf("reading index in GVTIDX\n");
   tidx=gvtidx_alloc(&gvh);
   if ((fp=gvfopen(p="GVTIDX.TAB","r"))==NULL)
      error(-1,errno,"[gvpds2mat] gvfopen %s failed",p);
   if (gvtidx_read(fp,&gvh,tidx))
      error(-1,0,"[gvpds2mat] gvtidx_read failed");
   fclose(fp);
	       
   printf("allocating arrays\n");
   if (match.dat) pidx=gvpidx_alloc(&gvh);
   if (match.gx) xif=gvxif_alloc(&gvh);
   if (match.ge) edf=gvedf_alloc(&gvh);
   if (match.gd) adf=gvadf_alloc(&gvh);
   if (match.ga) anf=gvanf_alloc(&gvh);
   gvpds_alloc(&gvh,&gv);

   for (ty=tylo; ty<=tyhi; ty++) {
      ylo = (ty==tylo) ? pylo : 0;
      yhi = (ty==tyhi) ? pyhi : npy-1;
      for (tx=txlo; tx<=txhi; tx++) {
	 printf("\nProcessing tile %d,%d\n",tx,ty);
	 if (match.dat) {
	    printf("   reading tile from PIDX\n");
	    if ((fp=gvfopen(p="GVPIDX.TAB","r"))==NULL)
	       error(-1,errno,"[gvpds2mat] gvfopen %s failed",p);
	    if (err=gvpidx_rtile(fp,&gvh,tidx,tx,ty,pidx))
	       error(-1,0,"[gvpds2mat] gvpidx_rtile failed %d",err);
	    fclose(fp);

	    if (match.gx) {
	       printf("   reading tile from XIF\n");
	       if ((fp=gvfopen(p="GVXIF.TAB","r"))==NULL)
		  error(-1,errno,"[gvpds2mat] gvfopen %s failed",p);
	       if (err=gvxif_rtile(fp,&gvh,tidx,tx,ty,xif))
		  error(-1,0,"[gvpds2mat] gvxif_rtile failed %d",err);
	       fclose(fp);
	    }

	    if (match.ge) {
	       printf("   reading tile from EDF\n");
	       if ((fp=gvfopen(p="GVRDF.TAB","r"))==NULL)
		  error(-1,errno,"[gvpds2mat] gvfopen %s failed",p);
	       if (err=gvedf_rtile(fp,&gvh,tidx,tx,ty,edf))
		  error(-1,0,"[gvpds2mat] gvedf_rtile failed %d",err);
	       fclose(fp);
	    }

	    if (match.gd) {
	       printf("   reading tile from ADF\n");
	       if ((fp=gvfopen(p="GVADF.TAB","r"))==NULL)
		  error(-1,errno,"[gvpds2mat] gvfopen %s failed",p);
	       if (err=gvadf_rtile(fp,&gvh,tidx,tx,ty,adf))
		  error(-1,0,"[gvpds2mat] gvadf_rtile failed %d",err);
	       fclose(fp);
	    }

	    if (match.ga) {
	       printf("   reading tile from ANF\n");
	       if ((fp=gvfopen(p="GVANF.TAB","r"))==NULL)
		  error(-1,errno,"[gvpds2mat] gvfopen %s failed",p);
	       if (err=gvanf_rtile(fp,&gvh,tidx,tx,ty,anf))
		  error(-1,0,"[gvpds2mat] gvanf_rtile failed %d",err);
	       fclose(fp);
	    }
	 }
	 printf("   extracting desired variables\n");
	 xlo = (tx==txlo) ? pxlo : 0;
	 xhi = (tx==txhi) ? pxhi : npx-1;
	 for (py=ylo; py<=yhi; py++) {
	    for (px=xlo; px<=xhi; px++) {
	       if (gvpds_rpixel(&gvh,pidx,xif,edf,adf,anf,px,py,&gv))
		  error(-1,0,"[gvpds2mat] gvpds_rpixel failed");
	       gvpds_t2p(&gvh,tx,ty,px,py,&x,&y);
	       gv.gv_X=x;
	       gv.gv_Y=y;
	       xr=px-pxlo+(tx-txlo)*npx;
	       yr=py-pylo+(ty-tylo)*npy;
	       tgvdrmat_pixel(xr,yr,&gv,&match);
	    }
	 }
      }
   }
   printf("\nwriting output\n");
   tgvdrmat_finish(fout);
   fclose(fout);
}

static void
get_bbox(gvdrhdr_t *hdr)
{
   if (!geog) {						/* specified in cartesian coords */
      if (bbxlo==INT_MIN || bbxlo<hdr->hdr3.gh_proj_bbxlo)
	 bbxlo=hdr->hdr3.gh_proj_bbxlo;
      if (bbxhi==INT_MAX || bbxhi>hdr->hdr3.gh_proj_bbxhi)
	 bbxhi=hdr->hdr3.gh_proj_bbxhi;
      if (bbylo==INT_MIN || bbylo<hdr->hdr3.gh_proj_bbylo)
	 bbylo=hdr->hdr3.gh_proj_bbylo;
      if (bbyhi==INT_MAX || bbyhi>hdr->hdr3.gh_proj_bbyhi)
	 bbyhi=hdr->hdr3.gh_proj_bbyhi;
   } else {						/* specified in geographic coords */
      double latlo=bbylo;
      double lathi=bbyhi;
      double lonlo=bbxlo;
      double lonhi=bbxhi;
      int xul,yul,xur,yur,xll,yll,xlr,ylr,err;
      gvdrcell_t gv;
      proj_t P;

      if (lonlo==INT_MIN || lonlo<hdr->hdr3.gh_proj_minlon)
	 lonlo=hdr->hdr3.gh_proj_minlon;
      if (lonhi==INT_MAX || lonhi>hdr->hdr3.gh_proj_maxlon)
	 lonhi=hdr->hdr3.gh_proj_maxlon;
      if (latlo==INT_MIN || latlo<hdr->hdr3.gh_proj_minlat)
	 latlo=hdr->hdr3.gh_proj_minlat;
      if (lathi==INT_MAX || lathi>hdr->hdr3.gh_proj_maxlat)
	 lathi=hdr->hdr3.gh_proj_maxlat;

      if (err=gvdr_rproj(&hdr->hdr3,&P))
	 error(-1,0,"[get_bbox] gvdr_rproj failed %d",err);

      /* upper left */
      gv.gv_lon=lonlo;
      gv.gv_lat=lathi;
      if (err=gvdr_fromlat(&P,&gv))
	 error(-1,0,"[get_bbox] gvdr_fromlat failed %d",err);
      xul=gv.gv_X;
      yul=gv.gv_Y;

      /* upper right */
      gv.gv_lon=lonhi;
      gv.gv_lat=lathi;
      if (err=gvdr_fromlat(&P,&gv))
	 error(-1,0,"[get_bbox] gvdr_fromlat failed %d",err);
      xur=gv.gv_X;
      yur=gv.gv_Y;

      /* lower left */
      gv.gv_lon=lonlo;
      gv.gv_lat=latlo;
      if (err=gvdr_fromlat(&P,&gv))
	 error(-1,0,"[get_bbox] gvdr_fromlat failed %d",err);
      xll=gv.gv_X;
      yll=gv.gv_Y;

      /* lower right */
      gv.gv_lon=lonhi;
      gv.gv_lat=latlo;
      if (err=gvdr_fromlat(&P,&gv))
	 error(-1,0,"[get_bbox] gvdr_fromlat failed %d",err);
      xlr=gv.gv_X;
      ylr=gv.gv_Y;

      bbxlo=xul;
      if (xur<bbxlo) bbxlo=xur;
      if (xll<bbxlo) bbxlo=xll;
      if (xlr<bbxlo) bbxlo=xlr;

      bbylo=yul;
      if (yur<bbylo) bbylo=yur;
      if (yll<bbylo) bbylo=yll;
      if (ylr<bbylo) bbylo=ylr;

      bbxhi=xul;
      if (xur>bbxhi) bbxhi=xur;
      if (xll>bbxhi) bbxhi=xll;
      if (xlr>bbxhi) bbxhi=xlr;

      bbyhi=yul;
      if (yur>bbyhi) bbyhi=yur;
      if (yll>bbyhi) bbyhi=yll;
      if (ylr>bbyhi) bbyhi=ylr;

      printf("Translated geographic limits to X=(%d:%d) Y=(%d:%d)\n",
	     bbxlo,bbxhi,bbylo,bbyhi);
   }
}

static void
init(int argc, char *argv[])
{
   extern int optind, opterr;
   extern char *optarg;
   int c, err=0, getopt();
   char *var;

   /* Get options and arguments */
   progn = argv[0];
   while ((c = getopt(argc, argv, "r:L:R:T:B:gl:u:")) != EOF)
      switch (c) {
       case 'r':
	 regexp = optarg;
	 break;
       case 'L':
	 bbxlo = strtol(optarg,NULL,10);
	 break;
       case 'R':
	 bbxhi = strtol(optarg,NULL,10);
	 break;
       case 'T':
	 bbyhi = strtol(optarg,NULL,10);
	 break;
       case 'B':
	 bbylo = strtol(optarg,NULL,10);
	 break;
       case 'g':
	 geog = !geog;
	 break;
       case 'l':
	 if (nre_lo==NRE) {
	    error(0,0,"[init] too many -l options");
	    err++;
	    break;
	 }
	 re_lo[nre_lo].regexp = strtok(optarg,"=");
	 var=strtok(NULL,"=");
	 if (!var || sscanf(var,"%f",&re_lo[nre_lo].val)!=1) {
	    error(0,0,"[init] could not parse %s",optarg);
	    err++;
	    break;
	 }
	 nre_lo++;
	 break;
       case 'u':
	 if (nre_hi==NRE) {
	    error(0,0,"[init] too many -l options");
	    err++;
	    break;
	 }
	 re_hi[nre_hi].regexp = strtok(optarg,"=");
	 var=strtok(NULL,"=");
	 if (!var || sscanf(var,"%f",&re_hi[nre_hi].val)!=1) {
	    error(0,0,"[init] could not parse %s",optarg);
	    err++;	
	    break;
	 }
	 nre_hi++;
	 break;
       default:
	 err++;
      }
   if (err) {
      int i;
      fprintf(stderr, "Usage: %s [options] [outfile]\n", progn);
      fprintf(stderr, "Available options: [current setting]\n");
      fprintf(stderr, "	-r regexp     Variable name regular expression [%s]\n", regexp);
      fprintf(stderr, "	-L left       Left side bounding box [%d]\n", bbxlo);
      fprintf(stderr, "	-R right      Right side bounding box [%d]\n", bbxhi);
      fprintf(stderr, "	-T top        Top side bounding box [%d]\n", bbyhi);
      fprintf(stderr, "	-B bottom     Bottom side bounding box [%d]\n", bbylo);
      fprintf(stderr, "	-g            Bounding box in geographic coordinates (lon/lat) [%d]\n", geog);
      fprintf(stderr, "	-l var=lo     Set cohort selection lower limit for var to lo\n");
      for (i=0; i<nre_lo; i++)
	 fprintf(stderr, "	   %10s = %g\n",re_lo[i].regexp,re_lo[i].val);
      fprintf(stderr, "	-u var=hi     Set cohort selection upper limit for var to hi\n");
      for (i=0; i<nre_hi; i++)
	 fprintf(stderr, "	   %10s = %g\n",re_hi[i].regexp,re_hi[i].val);
      exit(-1);
   }
   if (optind<argc) {
      if (strcmp(argv[optind],"-") && strcmp(argv[optind],"stdout")) {
	 if ((fout=fopen(argv[optind],"w"))==NULL)
	    error(-1,errno,"[init] fopen %s failed",argv[optind]);
      } else
	 fout=stdout;
      optind++;
   }
}