shpgeo.c

/******************************************************************************
 * Copyright (c) 1999, Carl Anderson
 *
 * This code is based in part on the earlier work of Frank Warmerdam
 * 
 * Permission is hereby granted, free of charge, to any person obtaining a
 * copy of this software and associated documentation files (the "Software"),
 * to deal in the Software without restriction, including without limitation
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 * and/or sell copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included
 * in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 * DEALINGS IN THE SOFTWARE.
 ******************************************************************************
 *
 * requires shapelib 1.2
 *   gcc shpproj shpopen.o dbfopen.o -lm -lproj -o shpproj
 * 
 * this may require linking with the PROJ4 projection library available from
 *
 * http://www.remotesensing.org/proj
 *
 * use -DPROJ4 to compile in Projection support
 *
 * $Log: shpgeo.c,v $
 * Revision 1.16  2017-07-10 18:01:35  erouault
 * * contrib/shpgeo.c: fix compilation on _MSC_VER < 1800 regarding lack
 * of NAN macro.
 *
 * Revision 1.15  2016-12-06 21:13:33  erouault
 * * configure.ac: change soname to 2:1:0 to be in sync with Debian soname.
 * http://bugzilla.maptools.org/show_bug.cgi?id=2628
 * Patch by Bas Couwenberg
 *
 * * contrib/doc/Shape_PointInPoly_README.txt, contrib/shpgeo.c: typo fixes.
 * http://bugzilla.maptools.org/show_bug.cgi?id=2629
 * Patch by Bas Couwenberg
 *
 * * web/*: use a local .css file to avoid a privacy breach issue reported
 * by the lintian QA tool.
 * http://bugzilla.maptools.org/show_bug.cgi?id=2630
 * Patch by Bas Couwenberg
 *
 *
 * Contributed by Sandro Mani: https://github.com/manisandro/shapelib/tree/autotools
 *
 * Revision 1.14  2016-12-05 12:44:07  erouault
 * * Major overhaul of Makefile build system to use autoconf/automake.
 *
 * * Warning fixes in contrib/
 *
 * Revision 1.13  2011-07-24 03:17:46  fwarmerdam
 * include string.h and stdlib.h where needed in contrib (#2146)
 *
 * Revision 1.12  2007-09-03 23:17:46  fwarmerdam
 * fix SHPDimension() function
 *
 * Revision 1.11  2006/11/06 20:45:58  fwarmerdam
 * Fixed SHPProject.
 *
 * Revision 1.10  2006/11/06 20:44:58  fwarmerdam
 * SHPProject() uses pj_transform now
 *
 * Revision 1.9  2006/01/25 15:33:50  fwarmerdam
 * fixed ppsC assignment maptools bug 1263
 *
 * Revision 1.8  2002/01/15 14:36:56  warmerda
 * upgrade to use proj_api.h
 *
 * Revision 1.7  2002/01/11 15:22:04  warmerda
 * fix many warnings.  Lots of this code is cruft.
 *
 * Revision 1.6  2001/08/30 13:42:31  warmerda
 * avoid use of auto initialization of PT for VC++
 *
 * Revision 1.5  2000/04/26 13:24:06  warmerda
 * made projUV handling safer
 *
 * Revision 1.4  2000/04/26 13:17:15  warmerda
 * check if projUV or UV
 *
 * Revision 1.3  2000/03/17 14:15:16  warmerda
 * Don't try to use system nan.h ... doesn't always exist.
 *
 * Revision 1.2  1999/05/26 02:56:31  candrsn
 * updates to shpdxf, dbfinfo, port from Shapelib 1.1.5 of dbfcat and shpinfo
 *
 */

#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "shapefil.h"

#include "shpgeo.h"

#if defined(_MSC_VER) && _MSC_VER < 1800
#include <float.h>
#define INFINITY (DBL_MAX + DBL_MAX)
#define NAN (INFINITY - INFINITY)
#endif


 /* I'm using some shorthand throughout this file
 *      R+ is a Clockwise Ring and is the positive portion of an object
 *      R- is a CounterClockwise Ring and is a hole in a R+
 *      A complex object is one having at least one R-
 *      A compound object is one having more than one R+
 *	A simple object has one and only one element (R+ or R-)
 *
 *	The closed ring constraint is for polygons and assumed here
 *	Arcs or LineStrings I am calling Rings (generically open or closed)
 *	Point types are vertices or lists of vertices but not Rings
 *
 *   SHPT_POLYGON, SHPT_POLYGONZ, SHPT_POLYGONM and SHPT_MULTIPATCH
 *   can have SHPObjects that are compound as well as complex
 *  
 *   SHP_POINT and its Z and M derivatives are strictly simple
 *   MULTI_POINT, SHPT_ARC and their derivatives may be simple or compound
 *
 */


/* **************************************************************************
 * asFileName
 *
 * utility function, toss part of filename after last dot
 *
 * **************************************************************************/ 
char * asFileName ( const char *fil, char *ext ) {
  char	pszBasename[120];
  static char	pszFullname[120];
  int	i;
/* -------------------------------------------------------------------- */
/*	Compute the base (layer) name.  If there is any extension	*/
/*	on the passed in filename we will strip it off.			*/
/* -------------------------------------------------------------------- */
//    pszFullname = (char*) malloc(( strlen(fil)+5 ));
//    pszBasename = (char *) malloc(strlen(fil)+5);
    strcpy( pszBasename, fil );
    for( i = strlen(pszBasename)-1; 
	 i > 0 && pszBasename[i] != '.' && pszBasename[i] != '/'
	       && pszBasename[i] != '\\';
	 i-- ) {}

    if( pszBasename[i] == '.' )
        pszBasename[i] = '\0';

/* -------------------------------------------------------------------- */
/*	Note that files pulled from										*/
/*	a PC to Unix with upper case filenames won't work!				*/
/* -------------------------------------------------------------------- */
//    pszFullname = (char *) malloc(strlen(pszBasename) + 5);
    sprintf( pszFullname, "%s.%s", pszBasename, ext );

    return ( pszFullname );
}


/************************************************************************/
/*                             SfRealloc()                              */
/*                                                                      */
/*      A realloc cover function that will access a NULL pointer as     */
/*      a valid input.                                                  */
/************************************************************************/
/* copied directly from shpopen.c -- maybe expose this in shapefil.h	*/  
static void * SfRealloc( void * pMem, int nNewSize )

{
    if( pMem == NULL )
        return( (void *) malloc(nNewSize) );
    else
        return( (void *) realloc(pMem,nNewSize) );
}


/* **************************************************************************
 * SHPPRoject
 *
 * Project points using projection handles, for use with PROJ4.3
 *
 * act as a wrapper to protect against library changes in PROJ
 * 
 * **************************************************************************/ 
int SHPProject ( SHPObject *psCShape, projPJ inproj, projPJ outproj ) {
#ifdef	PROJ4

    int    j;

    if ( pj_is_latlong(inproj) ) {
        for(j=0; j < psCShape->nVertices; j++) {
            psCShape->padfX[j] *= DEG_TO_RAD;
            psCShape->padfY[j] *= DEG_TO_RAD;
        }
    }   

    pj_transform(inproj, outproj, psCShape->nVertices, 0, psCShape->padfX,
                 psCShape->padfY, NULL);

    if ( pj_is_latlong(outproj) ) {
        for(j=0; j < psCShape->nVertices; j++) {
            psCShape->padfX[j] *= RAD_TO_DEG;
            psCShape->padfY[j] *= RAD_TO_DEG;
        }
    }   

    /* Recompute new Extents of projected Object								*/
    SHPComputeExtents ( psCShape );
#endif  

    return ( 1 );
}


/* **************************************************************************
 * SHPSetProjection
 *
 * establish a projection handle for use with PROJ4.3
 *
 * act as a wrapper to protect against library changes in PROJ
 *
 * **************************************************************************/
projPJ SHPSetProjection ( int param_cnt, char **params ) {
#ifdef PROJ4
  projPJ	*p = NULL;

  if ( param_cnt > 0 && params[0] )
  {
      p = pj_init ( param_cnt, params );
  }
  else
  {
      char* params_local[] = { "+proj=longlat", NULL };
      p = pj_init ( 1, params_local );
  }

  return ( p );
#else
  return ( NULL );
#endif
}


/* **************************************************************************
 * SHPFreeProjection
 *
 * release a projection handle for use with PROJ4.3
 *
 * act as a wrapper to protect against library changes in PROJ
 * 
 * **************************************************************************/
int SHPFreeProjection ( projPJ p) {
#ifdef PROJ4
  if ( p )
    pj_free ( p );
#endif
  return ( 1 );
}


/* **************************************************************************
 * SHPOGisType
 *
 * Convert Both ways from and to OGIS Geometry Types
 *
 * **************************************************************************/
int SHPOGisType ( int GeomType, int toOGis) {

    if ( toOGis == 0 )  /* connect OGis -> SHP types  					*/
	switch (GeomType) {
            case (OGIST_POINT):		return ( SHPT_POINT );      break;
            case (OGIST_LINESTRING):	return ( SHPT_ARC );        break;
            case (OGIST_POLYGON):		return ( SHPT_POLYGON );    break;
            case (OGIST_MULTIPOINT):	return ( SHPT_MULTIPOINT ); break;
            case (OGIST_MULTILINE):	return ( SHPT_ARC );	    break;
            case (OGIST_MULTIPOLYGON):	return ( SHPT_POLYGON );    break;
        }
    else  /* ok so its SHP->OGis types 									*/
        switch (GeomType) {
	    case (SHPT_POINT):		return ( OGIST_POINT );	    break;
	    case (SHPT_POINTM):		return ( OGIST_POINT );	    break;
	    case (SHPT_POINTZ):		return ( OGIST_POINT );	    break;
	    case (SHPT_ARC):		return ( OGIST_LINESTRING );break;
	    case (SHPT_ARCZ):		return ( OGIST_LINESTRING );break;
	    case (SHPT_ARCM):		return ( OGIST_LINESTRING );break;
	    case (SHPT_POLYGON):	return ( OGIST_MULTIPOLYGON );break;
	    case (SHPT_POLYGONZ):	return ( OGIST_MULTIPOLYGON );break;
	    case (SHPT_POLYGONM):	return ( OGIST_MULTIPOLYGON );break;
	    case (SHPT_MULTIPOINT):	return ( OGIST_MULTIPOINT );break;
	    case (SHPT_MULTIPOINTZ):	return ( OGIST_MULTIPOINT );break;
	    case (SHPT_MULTIPOINTM):	return ( OGIST_MULTIPOINT );break;
	    case (SHPT_MULTIPATCH):	return ( OGIST_GEOMCOLL );  break;
        } 	  	  	  

    return 0;
}


/* **************************************************************************
 * SHPReadSHPStream
 *
 * Encapsulate entire SHPObject for use with Postgresql
 *
 * **************************************************************************/
int SHPReadSHPStream ( SHPObject *psCShape, char *stream_obj) {

int	obj_storage;
int	my_order, need_swap =0, GeoType ;
int	use_Z = 0;
int	use_M = 0;

  need_swap = stream_obj[0];
  my_order = 1;
  my_order = ((char*) (&my_order))[0];
  need_swap = need_swap & my_order;

  if ( need_swap )
     swapW (stream_obj, (void*) &GeoType, sizeof (GeoType) );
   else
     memcpy (stream_obj, &GeoType, sizeof (GeoType) );

           
  if ( need_swap ) {
  
  } else {
    memcpy (stream_obj, &(psCShape->nSHPType),  sizeof (psCShape->nSHPType) );
    memcpy (stream_obj, &(psCShape->nShapeId),  sizeof (psCShape->nShapeId) );
    memcpy (stream_obj, &(psCShape->nVertices), sizeof (psCShape->nVertices) );
    memcpy (stream_obj, &(psCShape->nParts), 	 sizeof (psCShape->nParts) );
    memcpy (stream_obj, &(psCShape->dfXMin), 	 sizeof (psCShape->dfXMin) );
    memcpy (stream_obj, &(psCShape->dfYMin), 	 sizeof (psCShape->dfYMin) );
    memcpy (stream_obj, &(psCShape->dfXMax), 	 sizeof (psCShape->dfXMax) );
    memcpy (stream_obj, &(psCShape->dfYMax), 	 sizeof (psCShape->dfYMax) );
    if ( use_Z ) {
      memcpy (stream_obj, &(psCShape->dfZMin), 	 sizeof (psCShape->dfZMin) );
      memcpy (stream_obj, &(psCShape->dfZMax), 	 sizeof (psCShape->dfZMax) );
    }
   
    memcpy (stream_obj, psCShape->panPartStart, psCShape->nParts * sizeof (int) );
    memcpy (stream_obj, psCShape->panPartType,  psCShape->nParts * sizeof (int) );
    
/* get X and Y coordinate arrarys */    
    memcpy (stream_obj, psCShape->padfX, psCShape->nVertices * 2 * sizeof (double) );

/* get Z coordinate array if used */
    if ( use_Z ) 
      memcpy (stream_obj, psCShape->padfZ, psCShape->nVertices * 2 * sizeof (double) );
/* get Measure coordinate array if used */
    if ( use_M ) 
      memcpy (stream_obj, psCShape->padfM, psCShape->nVertices * 2 * sizeof (double) );
   }    /* end put data without swap */
   
   return (0);
}


/* **************************************************************************
 * SHPWriteSHPStream
 *
 * Encapsulate entire SHPObject for use with Postgresql
 *
 * **************************************************************************/
int SHPWriteSHPStream ( WKBStreamObj *stream_obj, SHPObject *psCShape ) {

/*int	obj_storage = 0;*/
int	need_swap = 0, my_order, GeoType;
int	use_Z = 0;
int	use_M = 0;

  need_swap = 1;
  need_swap = ((char*) (&need_swap))[0];
  
  /*realloc (stream_obj, obj_storage );*/
  
  if ( need_swap ) {
  
  } else {
    memcpy (stream_obj, psCShape, 4 * sizeof (int) );
    memcpy (stream_obj, psCShape, 4 * sizeof (double) );
    if ( use_Z ) 
      memcpy (stream_obj, psCShape, 2 * sizeof (double) );
    if ( use_M ) 
      memcpy (stream_obj, psCShape, 2 * sizeof (double) );
          
    memcpy (stream_obj, psCShape, psCShape->nParts * 2 * sizeof (int) );
    memcpy (stream_obj, psCShape, psCShape->nVertices * 2 * sizeof (double) );
    if ( use_Z ) 
      memcpy (stream_obj, psCShape, psCShape->nVertices * 2 * sizeof (double) );
    if ( use_M ) 
      memcpy (stream_obj, psCShape, psCShape->nVertices * 2 * sizeof (double) );
   }
   
  return (0);
}


/* **************************************************************************
 * WKBStreamWrite
 *
 * Encapsulate entire SHPObject for use with Postgresql
 *
 * **************************************************************************/
int WKBStreamWrite ( WKBStreamObj* wso, void* this, int tcount, int tsize ) {

   if ( wso->NeedSwap )
     SwapG ( &(wso->wStream[wso->StreamPos]), this, tcount, tsize );
   else
     memcpy ( &(wso->wStream[wso->StreamPos]), this, tsize * tcount );
     
   wso->StreamPos += tsize;

   return 0;
}



/* **************************************************************************
 * WKBStreamRead
 *
 * Encapsulate entire SHPObject for use with Postgresql
 *
 * **************************************************************************/
int WKBStreamRead ( WKBStreamObj* wso, void* this, int tcount, int tsize ) {

   if ( wso->NeedSwap )
     SwapG ( this, &(wso->wStream[wso->StreamPos]), tcount, tsize );
   else
     memcpy ( this, &(wso->wStream[wso->StreamPos]), tsize * tcount );
     
   wso->StreamPos += tsize;

   return 0;
}



/* **************************************************************************
 * SHPReadOGisWKB
 *
 * Encapsulate entire SHPObject for use with Postgresql
 *
 * **************************************************************************/
SHPObject* SHPReadOGisWKB ( WKBStreamObj *stream_obj) {
  SHPObject	*psCShape;
  char		WKB_order;
  int		need_swap = 0, my_order, GeoType = 0;
  int		use_Z = 0, use_M = 0;
  int		nSHPType, thisDim;

  WKBStreamRead ( stream_obj, &WKB_order, 1, sizeof(char));
  my_order = 1;
  my_order = ((char*) (&my_order))[0];
  stream_obj->NeedSwap = !(WKB_order & my_order);

  /* convert OGis Types to SHP types  */
  nSHPType = SHPOGisType ( GeoType, 0 );

  WKBStreamRead ( stream_obj, &GeoType, 1, sizeof(int));
     
  thisDim = SHPDimension ( nSHPType );
  
  if ( thisDim && SHPD_AREA ) 
    { psCShape = SHPReadOGisPolygon ( stream_obj ); } 
   else {
    if ( thisDim && SHPD_LINE )
      { psCShape = SHPReadOGisLine ( stream_obj ); }
    else {
      if ( thisDim && SHPD_POINT )
        { psCShape = SHPReadOGisPoint ( stream_obj ); }
      }
    }
   

   return (0);
}



/* **************************************************************************
 * SHPWriteOGisWKB
 *
 * Encapsulate entire SHPObject for use with Postgresql
 *
 * **************************************************************************/
int SHPWriteOGisWKB ( WKBStreamObj* stream_obj, SHPObject *psCShape ) {

  int		need_swap = 0, my_order, GeoType, thisDim;
  int		use_Z = 0, use_M = 0;  
  char		LSB = 1;
  /* indicate that this WKB is in LSB Order	*/

  /* OGis WKB can handle either byte order,  but if I get to choose I'd
  /* rather have it predicatable system-to-system							*/

  if ( stream_obj ) {
    if ( stream_obj->wStream )
      free ( stream_obj->wStream );
   } else 
    { stream_obj = calloc ( 3, sizeof (int ) ); }

  /* object size needs to be 9 bytes for the wrapper, and for each polygon	*/
  /* another 9 bytes all plus twice the total number of vertices 		*/
  /* times the sizeof (double) and just pad with 10 more chars for fun 		*/
  stream_obj->wStream = calloc (1, (9 * (psCShape->nParts + 1)) +
  	( sizeof(double) * 2 * psCShape->nVertices ) + 10 );
  
 #ifdef DEBUG2
   printf (" I just allocated %d bytes to wkbObj \n", 
  	(int)(sizeof (int) + sizeof (int) + sizeof(int) +
        ( sizeof(int) * psCShape->nParts + 1 ) +
  	( sizeof(double) * 2 * psCShape->nVertices ) + 10) );
 #endif 
 
  my_order = 1;
  my_order = ((char*) (&my_order))[0];
  /* Need to swap if this system is not  LSB (Intel Order)					*/
  stream_obj->NeedSwap = ( my_order != LSB );

  stream_obj->StreamPos = 0;

  
  #ifdef DEBUG2
    printf ("this system is (%d) LSB recorded as needSwap %d\n",my_order, stream_obj->NeedSwap);
  #endif

  WKBStreamWrite ( stream_obj, & LSB, 1, sizeof(char) );
  
  #ifdef DEBUG2
    printf ("this system in LSB \n");
  #endif

  
  /* convert SHP Types to OGis types  */
  GeoType = SHPOGisType ( psCShape->nSHPType, 1 );
  WKBStreamWrite ( stream_obj, &GeoType, 1, sizeof(int) );
     
  thisDim = SHPDimension ( psCShape->nSHPType );
  
  if ( thisDim && SHPD_AREA ) 
    { SHPWriteOGisPolygon ( stream_obj, psCShape ); } 
   else {
    if ( thisDim && SHPD_LINE )
      { SHPWriteOGisLine ( stream_obj, psCShape ); }
    else {
      if ( thisDim && SHPD_POINT )
        { SHPWriteOGisPoint ( stream_obj, psCShape ); }
      }
    }

#ifdef DEBUG2
  printf("(SHPWriteOGisWKB) outta here when stream pos is %d \n", stream_obj->StreamPos);
#endif   
  return (0);
}


/* **************************************************************************
 * SHPWriteOGisPolygon
 *
 * for this pass code to more generic OGis MultiPolygon Type
 * later add support for OGis Polygon Type
 *
 * Encapsulate entire SHPObject for use with Postgresql
 *
 * **************************************************************************/
int SHPWriteOGisPolygon ( WKBStreamObj *stream_obj, SHPObject *psCShape ) {
   SHPObject		**ppsC;
   SHPObject		*psC;
   int			rPart, ring, rVertices, cpart, cParts, nextring, i, j;
   char			Flag = 1;
   int			GeoType = OGIST_POLYGON;
   
   /* cant have more than nParts complex objects in this object */
   ppsC = calloc ( psCShape->nParts, sizeof(int) );

   
   nextring = 0;
   cParts=0;
   while ( nextring >= 0 ) {
       ppsC[cParts] = SHPUnCompound ( psCShape, &nextring ); 
       cParts++;
    }
   
#ifdef DEBUG2
     printf ("(SHPWriteOGisPolygon) Uncompounded into %d parts \n", cParts);
#endif
      
   WKBStreamWrite ( stream_obj, &cParts, 1, sizeof(int) );

   for ( cpart = 0; cpart < cParts; cpart++) {

     WKBStreamWrite ( stream_obj, & Flag, 1, sizeof(char) );
     WKBStreamWrite ( stream_obj, & GeoType, 1, sizeof(int) );
  
     psC = (SHPObject*) ppsC[cpart];
     WKBStreamWrite ( stream_obj, &(psC->nParts), 1, sizeof(int) );    
     
     for ( ring = 0; (ring < (psC->nParts)) && (psC->nParts > 0); ring ++) {
       if ( ring < (psC->nParts-2) )
         { rVertices = psC->panPartStart[ring+1] - psC->panPartStart[ring]; }
       else
         { rVertices = psC->nVertices - psC->panPartStart[ring]; }
#ifdef DEBUG2
     printf ("(SHPWriteOGisPolygon) scanning part %d, ring %d %d vtxs \n", 
     		cpart, ring, rVertices);
#endif             
       rPart = psC->panPartStart[ring];
       WKBStreamWrite ( stream_obj, &rVertices, 1, sizeof(int) );       
       for ( j=rPart; j < (rPart + rVertices); j++ ) {
         WKBStreamWrite ( stream_obj, &(psC->padfX[j]), 1, sizeof(double) );
         WKBStreamWrite ( stream_obj, &(psC->padfY[j]), 1, sizeof(double) );
        } /* for each vertex */
      }  /* for each ring */
    }  /* for each complex part */

#ifdef DEBUG2
     printf ("(SHPWriteOGisPolygon) outta here \n");
#endif   
  return (1);
}


/* **************************************************************************
 * SHPWriteOGisLine
 *
 * for this pass code to more generic OGis MultiXXXXXXX Type
 * later add support for OGis LineString Type
 *
 * Encapsulate entire SHPObject for use with Postgresql
 *
 * **************************************************************************/
int SHPWriteOGisLine ( WKBStreamObj *stream_obj, SHPObject *psCShape ) {

  return ( SHPWriteOGisPolygon( stream_obj, psCShape ));
}


/* **************************************************************************
 * SHPWriteOGisPoint
 *
 * for this pass code to more generic OGis MultiPoint Type
 * later add support for OGis Point Type
 *
 * Encapsulate entire SHPObject for use with Postgresql
 *
 * **************************************************************************/
int SHPWriteOGisPoint ( WKBStreamObj *stream_obj, SHPObject *psCShape ) {
   int			j;
   
   WKBStreamWrite ( stream_obj, &(psCShape->nVertices), 1, sizeof(int) );

   for ( j=0; j < psCShape->nVertices; j++ ) {
     WKBStreamWrite ( stream_obj, &(psCShape->padfX[j]), 1, sizeof(double) );
     WKBStreamWrite ( stream_obj, &(psCShape->padfY[j]), 1, sizeof(double) );
    } /* for each vertex */

  return (1);
}



/* **************************************************************************
 * SHPReadOGisPolygon
 *
 * for this pass code to more generic OGis MultiPolygon Type
 * later add support for OGis Polygon Type
 *
 * Encapsulate entire SHPObject for use with Postgresql
 *
 * **************************************************************************/
SHPObject* SHPReadOGisPolygon ( WKBStreamObj *stream_obj ) {
   SHPObject		**ppsC;
   SHPObject		*psC;
   int			rPart, ring, rVertices, cpart, cParts, nextring, i, j;
   int			totParts, totVertices, pRings, nParts;
   
   psC = SHPCreateObject ( SHPT_POLYGON, -1, 0, NULL, NULL, 0,
        	NULL, NULL, NULL, NULL ); 
   /* initialize a blank SHPObject 											*/        	

   WKBStreamRead ( stream_obj, &cParts, 1, sizeof(char) );

   totParts = cParts;
   totVertices = 0;
   
   SfRealloc ( psC->panPartStart, cParts * sizeof(int));
   SfRealloc ( psC->panPartType, cParts * sizeof(int));

   for ( cpart = 0; cpart < cParts; cpart++) {
     WKBStreamRead ( stream_obj, &nParts, 1, sizeof(int) );
     pRings = nParts;     
     /* pRings is the number of rings prior to the Ring loop below			*/
     
     if ( nParts > 1 ) {
       totParts += nParts - 1;
       SfRealloc ( psC->panPartStart, totParts * sizeof(int));
       SfRealloc ( psC->panPartType, totParts * sizeof(int));
      }

     rPart = 0;
     for ( ring = 0; ring < (nParts - 1); ring ++) {
       WKBStreamRead ( stream_obj, &rVertices, 1, sizeof(int) );
       totVertices += rVertices;
               
       psC->panPartStart[ring+pRings] = rPart;
       if ( ring == 0 )
         { psC->panPartType[ring + pRings] = SHPP_OUTERRING; }
        else
         { psC->panPartType[ring + pRings] = SHPP_INNERRING; }

       SfRealloc ( psC->padfX, totVertices * sizeof (double));
       SfRealloc ( psC->padfY, totVertices * sizeof (double));
              
       for ( j=rPart; j < (rPart + rVertices); j++ ) {
         WKBStreamRead ( stream_obj, &(psC->padfX[j]), 1, sizeof(double) );
         WKBStreamRead ( stream_obj, &(psC->padfY[j]), 1, sizeof(double) );        
        } /* for each vertex */
       rPart += rVertices;

      }  /* for each ring */
      
    }  /* for each complex part */
    
    return ( psC );
    
}


/* **************************************************************************
 * SHPReadOGisLine
 *
 * for this pass code to more generic OGis MultiLineString Type
 * later add support for OGis LineString Type
 * 
 * Encapsulate entire SHPObject for use with Postgresql
 *
 * **************************************************************************/
SHPObject* SHPReadOGisLine ( WKBStreamObj *stream_obj ) {
   SHPObject		**ppsC;
   SHPObject		*psC;
   int			rPart, ring, rVertices, cpart, cParts, nextring, i, j;
   int			totParts, totVertices, pRings, nParts;
   
   psC = SHPCreateObject ( SHPT_ARC, -1, 0, NULL, NULL, 0,
        	NULL, NULL, NULL, NULL ); 
   /* initialize a blank SHPObject 											*/        	

   WKBStreamRead ( stream_obj, &cParts, 1, sizeof(int) );

   totParts = cParts;
   totVertices = 0;
   
   SfRealloc ( psC->panPartStart, cParts * sizeof(int));
   SfRealloc ( psC->panPartType, cParts * sizeof(int));

   for ( cpart = 0; cpart < cParts; cpart++) {
     WKBStreamRead ( stream_obj, &nParts, 1, sizeof(int) );
     pRings = totParts;     
     /* pRings is the number of rings prior to the Ring loop below			*/
     
     if ( nParts > 1 ) {
       totParts += nParts - 1;
       SfRealloc ( psC->panPartStart, totParts * sizeof(int));
       SfRealloc ( psC->panPartType, totParts * sizeof(int));
      }

     rPart = 0;
     for ( ring = 0; ring < (nParts - 1); ring ++) {
       WKBStreamRead ( stream_obj, &rVertices, 1, sizeof(int) );
       totVertices += rVertices;
               
       psC->panPartStart[ring+pRings] = rPart;
       if ( ring == 0 )
         { psC->panPartType[ring + pRings] = SHPP_OUTERRING; }
        else
         { psC->panPartType[ring + pRings] = SHPP_INNERRING; }

       SfRealloc ( psC->padfX, totVertices * sizeof (double));
       SfRealloc ( psC->padfY, totVertices * sizeof (double));
              
       for ( j=rPart; j < (rPart + rVertices); j++ ) {
         WKBStreamRead ( stream_obj, &(psC->padfX[j]), 1, sizeof(double) );
         WKBStreamRead ( stream_obj, &(psC->padfY[j]), 1, sizeof(double) );        
        } /* for each vertex */
       rPart += rVertices;

      }  /* for each ring */
      
    }  /* for each complex part */
    
    return ( psC );
}


/* **************************************************************************
 * SHPReadOGisPoint
 *
 * Encapsulate entire SHPObject for use with Postgresql
 *
 * **************************************************************************/
SHPObject* SHPReadOGisPoint ( WKBStreamObj *stream_obj ) {
   SHPObject		*psC;
   int			nVertices, j;
   
   psC = SHPCreateObject ( SHPT_MULTIPOINT, -1, 0, NULL, NULL, 0,
        	NULL, NULL, NULL, NULL ); 
   /* initialize a blank SHPObject 											*/        	

   WKBStreamRead ( stream_obj, &nVertices, 1, sizeof(int) );

   SfRealloc ( psC->padfX, nVertices * sizeof (double));
   SfRealloc ( psC->padfY, nVertices * sizeof (double));
              
   for ( j=0; j < nVertices; j++ ) {
     WKBStreamRead ( stream_obj, &(psC->padfX[j]), 1, sizeof(double) );
     WKBStreamRead ( stream_obj, &(psC->padfY[j]), 1, sizeof(double) );        
    } /* for each vertex */
    
    return ( psC );
}




/* **************************************************************************
 * RingReadOGisWKB  
 *
 * this accepts OGisLineStrings which are basic building blocks
 *
 * Encapsulate entire SHPObject for use with Postgresql
 *
 * **************************************************************************/
int RingReadOgisWKB ( SHPObject *psCShape, char *stream_obj) {
    return 0;
}


/* **************************************************************************
 * RingWriteOGisWKB
 *
 * this emits OGisLineStrings which are basic building blocks
 * 
 * Encapsulate entire SHPObject for use with Postgresql
 *
 * **************************************************************************/
int RingWriteOgisWKB ( SHPObject *psCShape, char *stream_obj) {

    return 0;
}


/* **************************************************************************
 * SHPDimension
 *
 * Return the Dimensionality of the SHPObject 
 * a handy utility function
 * 
 * **************************************************************************/
int SHPDimension ( int SHPType ) {
    int 	dimension;
    
    dimension = 0;
      
    switch ( SHPType ) {
      case  SHPT_POINT       :	dimension = SHPD_POINT; break;
      case  SHPT_ARC         :	dimension = SHPD_LINE; break;  
      case  SHPT_POLYGON     :	dimension = SHPD_AREA; break;   	
      case  SHPT_MULTIPOINT  :	dimension = SHPD_POINT; break;
      case  SHPT_POINTZ      :	dimension = SHPD_POINT | SHPD_Z; break;
      case  SHPT_ARCZ        :	dimension = SHPD_LINE | SHPD_Z; break;
      case  SHPT_POLYGONZ    :	dimension = SHPD_AREA | SHPD_Z; break;
      case  SHPT_MULTIPOINTZ :	dimension = SHPD_POINT | SHPD_Z; break;
      case  SHPT_POINTM      :	dimension = SHPD_POINT | SHPD_MEASURE; break;
      case  SHPT_ARCM        :	dimension = SHPD_LINE | SHPD_MEASURE; break;
      case  SHPT_POLYGONM    :	dimension = SHPD_AREA | SHPD_MEASURE; break;
      case  SHPT_MULTIPOINTM :	dimension = SHPD_POINT | SHPD_MEASURE; break;
      case  SHPT_MULTIPATCH  :	dimension = SHPD_AREA; break;
    }

    return ( dimension );
}


/* **************************************************************************
 * SHPPointinPoly_2d
 *
 * Return a Point inside an R+ of a potentially 
 * complex/compound SHPObject suitable for labelling
 * return only one point even if if is a compound object
 *
 * reject non area SHP Types
 * 
 * **************************************************************************/
PT	SHPPointinPoly_2d ( SHPObject *psCShape ) {
   PT	*sPT, rPT;
   
   if ( !(SHPDimension (psCShape->nSHPType) & SHPD_AREA) )
   {
       rPT.x = NAN;
       rPT.y = NAN;
       return rPT;
   }

   sPT = SHPPointsinPoly_2d ( psCShape );
   
   if ( sPT ) {
     rPT.x = sPT[0].x;
     rPT.y = sPT[0].y; 
    } else {
     rPT.x = NAN;
     rPT.y = NAN;
    }
   return ( rPT );
}


/* **************************************************************************
 * SHPPointsinPoly_2d
 *
 * Return a Point inside each R+ of a potentially 
 * complex/compound SHPObject suitable for labelling
 * return one point for each R+ even if it is a compound object
 *
 * reject non area SHP Types
 * 
 * **************************************************************************/
PT*	SHPPointsinPoly_2d ( SHPObject *psCShape ) {
   PT		*PIP = NULL;
   int		cRing;
   SHPObject	*psO, *psInt, *CLine;
   double	*CLx, *CLy;
   int		*CLstt, *CLst, nPIP, ring, rMpart, ring_vtx, ring_nVertices;
   double	rLen, rLenMax;

   if ( !(SHPDimension (psCShape->nSHPType) & SHPD_AREA) )  
      return ( NULL );

   while (  psO = SHPUnCompound (psCShape, &cRing)) {
     CLx = calloc ( 4, sizeof(double));
     CLy = calloc ( 4, sizeof(double));
     CLst = calloc ( 2, sizeof(int));
     CLstt = calloc ( 2, sizeof(int));
     
     /* a horizontal & vertical compound line though the middle of the 		*/
     /* extents 															*/
     CLx [0] = psO->dfXMin;
     CLy [0] = (psO->dfYMin + psO->dfYMax ) * 0.5;
     CLx [1] = psO->dfXMax;   
     CLy [1] = (psO->dfYMin + psO->dfYMax ) * 0.5;
   
     CLx [2] = (psO->dfXMin + psO->dfXMax ) * 0.5;
     CLy [2] = psO->dfYMin;   
     CLx [3] = (psO->dfXMin + psO->dfXMax ) * 0.5;
     CLy [3] = psO->dfYMax;
     
     CLst[0] = 0;   CLst[1] = 2; 
     CLstt[0] = SHPP_RING; CLstt[1] = SHPP_RING;      
     
     CLine = SHPCreateObject ( SHPT_POINT, -1, 2, CLst, CLstt, 4,
        	CLx, CLy, NULL, NULL ); 

     /* with the H & V centrline compound object, intersect it with the OBJ	*/
     psInt = SHPIntersect_2d ( CLine, psO );
     /* return SHP type is lowest common dimensionality of the input types 	*/


     /* find the longest linestring returned by the intersection			*/
     ring_vtx = psInt->nVertices ;
     for ( ring = (psInt->nParts - 1); ring >= 0; ring-- ) {
       ring_nVertices = ring_vtx - psInt->panPartStart[ring];

       rLen += RingLength_2d ( ring_nVertices, 
     	(double*) &(psInt->padfX [psInt->panPartStart[ring]]), 
     	(double*) &(psInt->padfY [psInt->panPartStart[ring]]) );
 
       if ( rLen > rLenMax )  
         { rLenMax = rLen;
           rMpart = psInt->panPartStart[ring];
         }
       ring_vtx = psInt->panPartStart[ring];
      } 

     /* add the centerpoint of the longest ARC of the intersection to the	*/
     /*	PIP list															*/
     nPIP ++;
     SfRealloc ( PIP, sizeof(double) * 2 * nPIP);
     PIP[nPIP].x = (psInt ->padfX [rMpart] + psInt ->padfX [rMpart]) * 0.5;
     PIP[nPIP].y = (psInt ->padfY [rMpart] + psInt ->padfY [rMpart]) * 0.5;
     
     SHPDestroyObject ( psO );
     SHPDestroyObject ( CLine );