shpgeo.cpp

projapi是一个关于GIS行业投影转换的程序库

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

     /* 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;

	 return ( C );
}


/* **************************************************************************
 * RingCentroid_2d
 *
 * Return the mathematical / geometric centroid of a single closed ring
 *
 * **************************************************************************/
PROJAPI_CALL 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;

    ppx = x;
    ppy = y;
  }

  /* 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
 * **************************************************************************/
PROJAPI_CALL 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; }
   }

   /* 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) );

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



/* **************************************************************************
 * SHPArea_2d
 *
 * Calculate the XY Area of Polygon ( can be compound / complex )
 *
 * **************************************************************************/
PROJAPI_CALL 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								*/

   /* I just start at the last ring and work down to the first				*/
   ring_vtx = psCShape->nVertices ;
   for ( ring = (psCShape->nParts - 1); ring >= 0; ring-- ) 
   {
     ring_nVertices = ring_vtx - psCShape->panPartStart[ring];

     cArea += RingArea_2d ( ring_nVertices, 
     	(double*) &(psCShape->padfX [psCShape->panPartStart[ring]]), 
     	(double*) &(psCShape->padfY [psCShape->panPartStart[ring]]) );

     ring_vtx = psCShape->panPartStart[ring];
    } 

    /* Area is signed, negative Areas are R-									*/    
    return ( cArea );
    
}


/* **************************************************************************
 * SHPLength_2d
 *
 * Calculate the Planar ( XY ) Length of Polygon ( can be compound / complex )
 *    or Polyline ( can be compound ).  Length on Polygon is its Perimeter
 *
 * **************************************************************************/
PROJAPI_CALL double SHPLength_2d ( SHPObject *psCShape ) 
{
    double 	Length;
    int		i, j;
    double 	dx, dy;
    
   if ( !(SHPDimension (psCShape->nSHPType) & (SHPD_AREA || SHPD_LINE)) )  
       return ( (double) -1 );    
    
    Length = 0;
    j = 1;
    for ( i = 1; i < psCShape->nVertices; i++ ) {  
      if ( psCShape->panPartStart[j] == i ) 
       { j ++; }
    /* skip the moves with "pen up" from ring to ring */       
      else
       { 
        dx = psCShape->padfX[i] - psCShape->padfX[i-1];
        dy = psCShape->padfY[i] - psCShape->padfY[i-1];
        Length += sqrt ( ( dx * dx ) + ( dy * dy ) );
       }
     /* simplify this equation */
     }
     
   return ( Length );
}


/* **************************************************************************
 * RingLength_2d
 *
 * Calculate the Planar ( XY ) Length of Polygon ( can be compound / complex )
 *    or Polyline ( can be compound ).  Length of Polygon is its Perimeter
 *
 * **************************************************************************/
PROJAPI_CALL double RingLength_2d ( int nVertices, double *a, double *b ) 
{
    double 	Length;
    int		i, j;
    double 	dx, dy;
    
    Length = 0;
    j = 1;
    for ( i = 1; i < nVertices; i++ ) {  
      dx = a[i] - b[i-1];
      dy = b[i] - b[i-1];
      Length += sqrt ( ( dx * dx ) + ( dy * dy ) );
     /* simplify this equation */
     }
     
   return ( Length );
}


/* **************************************************************************
 * RingArea_2d
 *
 * Calculate the Planar Area of a single closed ring 
 *
 * **************************************************************************/
PROJAPI_CALL double RingArea_2d ( int nVertices, double *a, double *b ) 
{
  int		iv,jv;
  double	ppx, ppy;
  static 	double	Area;
  double	dx_Area;
  double 	x_base, y_base, x, y;
  
  x_base = a[0];
  y_base = b[0];
  
  ppx = a[1] - x_base;
  ppy = b[1] - y_base;
  Area = 0.0;

  for ( iv = 2; iv <= ( nVertices - 1 ); iv++ ) 
  {
    x = a[iv] - x_base;
    y = b[iv] - y_base;

    /* Area of a triangle is the cross product of its defining vectors		*/
    
    dx_Area = ((x * ppy) - (y * ppx)) * 0.5;

    Area += dx_Area;

    ppx = x;
    ppy = y;
  }

  return ( Area );
}



/* **************************************************************************
 * SHPUnCompound
 *
 * ESRI calls this function explode
 * Return a non compound ( possibly complex ) object
 *
 * ring_number is R+ number corresponding to object
 *
 *
 * ignore complexity in Z dimension for now
 *
 * **************************************************************************/
PROJAPI_CALL SHPObject * SHPUnCompound  ( SHPObject *psCShape, int * ringNumber ) 
{
	int ringDir, ring, lRing;

	if ( (*ringNumber >= psCShape->nParts) || *ringNumber == -1 ) 
	{
		*ringNumber = -1;	
		return (NULL);
	}
   
	if ( *ringNumber == (psCShape->nParts - 1) )  
	{
		*ringNumber =  -1;
		return ( SHPClone(psCShape, (psCShape->nParts - 1), -1) );
	}

	lRing = *ringNumber;
	ringDir = -1;
	for ( ring = (lRing + 1); (ring < psCShape->nParts) && ( ringDir < 0 ); ring ++)
		ringDir = SHPRingDir_2d ( psCShape, ring);
   
	if ( ring ==  psCShape->nParts )
		*ringNumber = -1; 
	else
		*ringNumber = ring;
/*    I am strictly assuming that all R- parts of a complex object 
 *	   directly follow their R+, so when we hit a new R+ its a
 *	   new part of a compound object
 *         a SHPClean may be needed to enforce this as it is not part
 *	   of ESRI's definition of a SHPfile
 */

	return ( SHPClone(psCShape, lRing, ring ) );
}


/* **************************************************************************
 * SHPIntersect_2d
 *
 * 
 * prototype only for now
 *
 * return object with lowest common dimensionality of objects
 * 
 * **************************************************************************/ 
PROJAPI_CALL SHPObject * SHPIntersect_2d ( SHPObject* a, SHPObject* b ) 
{
  SHPObject	*C;
  
  if ( (SHPDimension(a->nSHPType) && SHPD_POINT) || ( SHPDimension(b->nSHPType) && SHPD_POINT ) )
    return ( NULL );
  /* there is no intersect function like this for points  */
    
  C = SHPClone ( a, 0 , -1 );

  return ( C);
}



/* **************************************************************************
 * SHPClean
 *
 * Test and fix normalization problems in shapes
 * Different tests need to be implemented for different SHPTypes
 * 	SHPT_POLYGON	check ring directions CW / CCW   ( R+ / R- )
 *				put all R- after the R+ they are members of
 *				i.e. each complex object is completed before the
 *     				next object is started
 *				check for closed rings
 *				ring must not intersect itself, even on edge
 *
 *  no other types implemented yet
 * 
 * not sure why but return object in place
 * use for object casting and object verification						 
 * **************************************************************************/ 
PROJAPI_CALL int SHPClean ( SHPObject *psCShape ) 
{


    return (0);
}


/* **************************************************************************
 * SHPClone
 *
 * Clone a SHPObject, replicating all data
 * 
 * **************************************************************************/ 
PROJAPI_CALL SHPObject * SHPClone ( SHPObject *psCShape, int lowPart, int highPart ) 
{
    SHPObject	*psObject;
    int		newParts, newVertices;

    if ( highPart >= psCShape->nParts || highPart == -1 )
	highPart = psCShape->nParts ;

    newParts = highPart - lowPart;
    if ( newParts == 0 ) { return ( NULL ); }
    
    psObject = (SHPObject *) calloc(1,sizeof(SHPObject));
    psObject->nSHPType = psCShape->nSHPType;
    psObject->nShapeId = psCShape->nShapeId;
    
    psObject->nParts = newParts;
    if ( psCShape->padfX ) {
        psObject->panPartStart = (int*) calloc (newParts, sizeof (int)); 
        memcpy ( psObject->panPartStart, psCShape->panPartStart, 
      		newParts * sizeof (int) );
     }
    if ( psCShape->padfX ) {
      psObject->panPartType = (int*) calloc (newParts, sizeof (int)); 
      memcpy ( psObject->panPartType, 
  		(int *) &(psCShape->panPartType[lowPart]), 
      		newParts * sizeof (int) );
     }

    if ( highPart != psCShape->nParts ) {
      newVertices = psCShape->panPartStart[highPart] -
	 psCShape->panPartStart[lowPart]; 
     }
    else 
     { newVertices = psCShape->nVertices - psCShape->panPartStart[lowPart]; }


    psObject->nVertices = newVertices;
    if ( psCShape->padfX ) {
      psObject->padfX = (double*) calloc (newVertices, sizeof (double)); 
      memcpy ( psObject->padfX,
	 (double *) &(psCShape->padfX[psCShape->panPartStart[lowPart]]),
      		newVertices * sizeof (double) );
     } 
    if ( psCShape->padfY ) {
      psObject->padfY = (double*) calloc (newVertices, sizeof (double)); 
      memcpy ( psObject->padfY, 
	 (double *) &(psCShape->padfY[psCShape->panPartStart[lowPart]]),
      		newVertices * sizeof (double) );
     }
    if ( psCShape->padfZ ) {
      psObject->padfZ = (double*) calloc (newVertices, sizeof (double)); 
      memcpy ( psObject->padfZ, 
	 (double *) &(psCShape->padfZ[psCShape->panPartStart[lowPart]]),
      		newVertices * sizeof (double) );
     }
    if ( psCShape->padfM ) {
      psObject->padfM = (double*) calloc (newVertices, sizeof (double)); 
      memcpy ( psObject->padfM,
	(double *) &(psCShape->padfM[psCShape->panPartStart[lowPart]]),
      		newVertices * sizeof (double) );
     }

    psObject->dfXMin = psCShape->dfXMin;
    psObject->dfYMin = psCShape->dfYMin;
    psObject->dfZMin = psCShape->dfZMin;
    psObject->dfMMin = psCShape->dfMMin;

    psObject->dfXMax = psCShape->dfXMax;
    psObject->dfYMax = psCShape->dfYMax;
    psObject->dfZMax = psCShape->dfZMax;
    psObject->dfMMax = psCShape->dfMMax;

    SHPComputeExtents ( psObject );
    return ( psObject );
}