/*-*- mode: C; c-basic-offset:2 ; tab-width: 2; -*-
Simple hack to loop through a shapefile and extract only those shapes
  that lie within a specified bounding box

  Coordinate system of the bounding box specification must match the 
  coordinate system of the shapefile.  No checking of this requirement is done.
*/
#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include <string.h>
#include <math.h>
#include <shapefil.h>

int main(int argc, char **argv)
{

  char inshape[120] = "";
  char outshape[120] = "";
  double xmin,xmax,ymin,ymax;

  SHPHandle InShape=NULL, OutShape=NULL;
  DBFHandle InDBF=NULL, OutDBF=NULL;
  int nShapeType, nEntities, nVertices, *panParts, i, iPart,j;
  SHPObject	*psCShape;
  char	*DBFRow = NULL;

  if (argc < 7)
    {
      fprintf(stderr,"Usage: %s inshape outshape xmin xmax ymin ymax\n",
              argv[0]);
      exit(0);
    }

    strncpy( inshape, *++argv, sizeof(inshape));
    --argc;
    strncpy( outshape, *++argv, sizeof(outshape));
    --argc; 
    printf("In: %s, out: %s\n",inshape,outshape);

    xmin=atof(*++argv);
    --argc;
    xmax=atof(*++argv);
    --argc;
    ymin=atof(*++argv);
    --argc;
    ymax=atof(*++argv);
    --argc;

    printf("Selecting shapes inside of bounding box X=(%lf,%lf), Y=(%lf,%lf)\n",xmin,xmax,ymin,ymax);
    if (argc != 1)
      {
        fprintf(stderr,"Usage: %s inshape outshape xmin xmax ymin ymax\n",
              argv[0]);
        fprintf(stderr,"Ignoring %d extra arguments.\n",argc-1);
      }

  InShape = SHPOpen(inshape, "rb");
  InDBF = DBFOpen( inshape, "rb");

  if (InShape == NULL || InDBF == NULL)
    {
      fprintf (stderr,"Unable to open source shapefiles: %s\n", inshape);
      exit(1);
    }
  
  SHPGetInfo(InShape, &nEntities, &nShapeType, NULL, NULL);
  OutShape = SHPCreate( outshape, nShapeType);
  OutDBF = DBFCloneEmpty(InDBF, outshape);
  if (OutShape == NULL || OutDBF == NULL)
    {
      fprintf (stderr,"Unable to open destination shapefiles: %s\n", outshape);
      exit(1);
    }
  DBFRow = (char *) malloc ( (InDBF->nRecordLength) + 15 );

  for( i = 0; i < nEntities; i++ )
    {
      psCShape = SHPReadObject ( InShape, i );
      
      // This is meant to match the way xastir selects bounding boxes
      // as compared to its view region
      if (ymin <= psCShape->dfYMax
          && ymax >= psCShape->dfYMin
          && xmin <= psCShape->dfXMax
          && xmax >= psCShape->dfXMin)
        {
          SHPWriteObject ( OutShape, -1, psCShape );   
          SHPDestroyObject ( psCShape );
          
          memcpy ( DBFRow, DBFReadTuple ( InDBF, i ), InDBF->nRecordLength );
          DBFWriteTuple ( OutDBF, OutDBF->nRecords, DBFRow );
        }
    }
  SHPClose(InShape);
  SHPClose(OutShape);
  DBFClose(InDBF);
  DBFClose(OutDBF);

}