shpgeo.c

Source code, and some other odds and ends can be downloaded from http://shapelib.maptools.org/dl.

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

}


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


}


/* **************************************************************************
 * 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 ;
    	case  SHPT_ARC         :	dimension = SHPD_LINE;
    	case  SHPT_POLYGON     :	dimension = SHPD_AREA;    	
    	case  SHPT_MULTIPOINT  :	dimension = SHPD_POINT;    	
    	case  SHPT_POINTZ      :	dimension = SHPD_POINT | SHPD_Z;    	
    	case  SHPT_ARCZ        :	dimension = SHPD_LINE | SHPD_Z;
    	case  SHPT_POLYGONZ    :	dimension = SHPD_AREA | SHPD_Z;
    	case  SHPT_MULTIPOINTZ :	dimension = SHPD_POINT | SHPD_Z;    	
    	case  SHPT_POINTM      :	dimension = SHPD_POINT | SHPD_MEASURE;    	
    	case  SHPT_ARCM        :	dimension = SHPD_LINE | SHPD_MEASURE;
    	case  SHPT_POLYGONM    :	dimension = SHPD_AREA | SHPD_MEASURE;
    	case  SHPT_MULTIPOINTM :	dimension = SHPD_POINT | SHPD_MEASURE;   	 	
    	case  SHPT_MULTIPATCH  :	dimension = SHPD_AREA;  
    }

   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) )  
       return ( (PT) {NAN,NAN} );

   sPT = SHPPointsinPoly_2d ( psCShape );
   
   if ( sPT ) {
     rPT.x = sPT[0].x;
     rPT.y = sPT[0].y; 
    } else {
     rPT =  (( PT ){ NAN, 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;
   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 );
     
     /* does SHPCreateobject use preallocated memory or does it copy the 	*/
     /* contents.  To be safe conditionally release CLx, CLy, CLst, CLstt	*/
     if ( CLx )   free ( CLx );
     if ( CLy )   free ( CLy );
     if ( CLst )  free ( CLst );
     if ( CLstt ) free ( CLstt );
    }
    
  return ( PIP );
}


/* **************************************************************************
 * SHPCentrd_2d
 *
 * Return the single mathematical / geometric centroid of a potentially 
 * complex/compound SHPObject
 *
 * reject non area SHP Types
 * 
 * **************************************************************************/
PT SHPCentrd_2d ( SHPObject *psCShape ) {
    int		ring, ringPrev, ring_nVertices, rStart;
    double	Area, ringArea;
    PT		ringCentrd, C;
    
  
   if ( !(SHPDimension (psCShape->nSHPType) & SHPD_AREA) )  
       return ( (PT) {NAN,NAN} );

#ifdef DEBUG
	printf ("for Object with %d vtx, %d parts [ %d, %d] \n",
		psCShape->nVertices, psCShape->nParts,
		psCShape->panPartStart[0],psCShape->panPartStart[1]);
#endif
   
   Area = 0;
   C.x = 0.0;
   C.y = 0.0;
   
   /* for each ring in compound / complex object calc the ring cntrd		*/
   
   ringPrev = psCShape->nVertices;
   for ( ring = (psCShape->nParts - 1); ring >= 0; ring-- ) {
     rStart = psCShape->panPartStart[ring];
     ring_nVertices = ringPrev - rStart;

     RingCentroid_2d ( ring_nVertices, (double*) &(psCShape->padfX [rStart]), 
     	(double*) &(psCShape->padfY [rStart]), &ringCentrd, &ringArea);  

#ifdef DEBUG
	printf ("(SHPCentrd_2d)  Ring %d, vtxs %d, area: %f, ring centrd %f, %f \n",
		ring, ring_nVertices, ringArea, ringCentrd.x, ringCentrd.y);
#endif

     /* use Superposition of these rings to build a composite Centroid		*/
     /* sum the ring centrds * ringAreas,  at the end divide by total area	*/
     C.x +=  ringCentrd.x * ringArea;
     C.y +=  ringCentrd.y * ringArea; 
     Area += ringArea; 
     ringPrev = rStart;
    }    

     /* hold on the division by AREA until were at the end					*/
     C.x = C.x / Area;
     C.y = C.y / Area;
#ifdef DEBUG
	printf ("SHPCentrd_2d) Overall Area: %f, Centrd %f, %f \n",
		Area, C.x, C.y);
#endif   
   return ( C );
}


/* **************************************************************************
 * RingCentroid_2d
 *
 * Return the mathematical / geometric centroid of a single closed ring
 *
 * **************************************************************************/
int RingCentroid_2d ( int nVertices, double *a, double *b, PT *C, double *Area ) {
  int		iv,jv;
  int		sign_x, sign_y;
  double	dy_Area, dx_Area, Cx_accum, Cy_accum, ppx, ppy;
  double 	x_base, y_base, x, y;

/* the centroid of a closed Ring is defined as
 *
 *      Cx = sum (cx * dArea ) / Total Area
 *  and
 *      Cy = sum (cy * dArea ) / Total Area
 */      
   
  x_base = a[0];
  y_base = b[0];
  
  Cy_accum = 0.0;
  Cx_accum = 0.0;

  ppx = a[1] - x_base;
  ppy = b[1] - y_base;
  *Area = 0;

/* Skip the closing vector */
  for ( iv = 2; iv <= nVertices - 2; iv++ ) {
    x = a[iv] - x_base;
    y = b[iv] - y_base;

    /* calc the area and centroid of triangle built out of an arbitrary 	*/
    /* base_point on the ring and each successive pair on the ring			*/
    
    /* Area of a triangle is the cross product of its defining vectors		*/
    /* Centroid of a triangle is the average of its vertices				*/

    dx_Area =  ((x * ppy) - (y * ppx)) * 0.5;
    *Area += dx_Area;
    
    Cx_accum += ( ppx + x ) * dx_Area;       
    Cy_accum += ( ppy + y ) * dx_Area;
#ifdef DEBUG2
    printf("(ringcentrd_2d)  Pp( %f, %f), P(%f, %f)\n", ppx, ppy, x, y);
    printf("(ringcentrd_2d)    dA: %f, sA: %f, Cx: %f, Cy: %f \n", 
		dx_Area, *Area, Cx_accum, Cy_accum);
#endif    
    ppx = x;
    ppy = y;
  }

#ifdef DEBUG2
  printf("(ringcentrd_2d)  Cx: %f, Cy: %f \n", 
  	( Cx_accum / ( *Area * 3) ), ( Cy_accum / (*Area * 3) ));
#endif

  /* adjust back to world coords 											*/
  C->x = ( Cx_accum / ( *Area * 3)) + x_base;
  C->y = ( Cy_accum / ( *Area * 3)) + y_base;
      
  return ( 1 );
}




/* **************************************************************************
 * SHPRingDir_2d
 *
 * Test Polygon for CW / CCW  ( R+ / R- )
 *
 * return 1  for R+
 * return -1 for R-
 * return 0  for error
 * **************************************************************************/
int SHPRingDir_2d ( SHPObject *psCShape, int Ring ) {
   int		i, ti, last_vtx;
   double	tX;
   double 	*a, *b;
   double	dx0, dx1, dy0, dy1, v1, v2 ,v3;
   
   tX = 0.0;
   a = psCShape->padfX;
   b = psCShape->padfY;
   
   if ( Ring >= psCShape->nParts ) return ( 0 );
   
   if ( Ring >= psCShape->nParts -1 )
     { last_vtx = psCShape->nVertices; }
    else
     { last_vtx = psCShape->panPartStart[Ring + 1]; }
      
   /* All vertices at the corners of the extrema (rightmost lowest, leftmost lowest, 	*/
   /* topmost rightest, ...) must be less than pi wide.  If they werent they couldnt be	*/
   /* extrema.																			*/   
   /* of course the following will fail if the Extents are even a little wrong 			*/      
   
   for ( i = psCShape->panPartStart[Ring]; i < last_vtx; i++ ) {
     if ( b[i] == psCShape->dfYMax && a[i] > tX ) 
      { ti = i; }
   }

#ifdef DEBUG2
   printf ("(shpgeo:SHPRingDir) highest Rightmost Pt is vtx %d (%f, %f)\n", ti, a[ti], b[ti]);
#endif   
   
   /* cross product */
   /* the sign of the cross product of two vectors indicates the right or left half-plane	*/
   /* which we can use to indicate Ring Dir													*/ 
   if ( ti > psCShape->panPartStart[Ring] & ti < last_vtx ) 
    { dx0 = a[ti-1] - a[ti];
      dx1 = a[ti+1] - a[ti];
      dy0 = b[ti-1] - b[ti];
      dy1 = b[ti+1] - b[ti];
   }
   else
   /* if the tested vertex is at the origin then continue from 0 */ 
   {  dx1 = a[1] - a[0];
      dx0 = a[last_vtx] - a[0];
      dy1 = b[1] - b[0];
      dy0 = b[last_vtx] - b[0];
   }
   
//   v1 = ( (dy0 * 0) - (0 * dy1) );
//   v2 = ( (0 * dx1) - (dx0 * 0) );
/* these above are always zero so why do the math */
   v3 = ( (dx0 * dy1) - (dx1 * dy0) );

#ifdef DEBUG2   
   printf ("(shpgeo:SHPRingDir)  cross product for vtx %d was %f \n", ti, v3); 
#endif

   if ( v3 > 0 )
    { return (1); }
   else
    { return (-1); }
}



/* **************************************************************************
 * SHPArea_2d
 *
 * Calculate the XY Area of Polygon ( can be compound / complex )
 *
 * **************************************************************************/
double SHPArea_2d ( SHPObject *psCShape ) {
    double 	cArea;
    int		ring, ring_vtx, ringDir, ring_nVertices;
  
   cArea = 0;
   if ( !(SHPDimension (psCShape->nSHPType) & SHPD_AREA) )  
       return ( -1 );
 
 
   /* Walk each ring adding its signed Area,  R- will return a negative 	*/
   /* area, so we don't have to test for them								*/