shpgeo.c

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

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

#ifdef DEBUG2
     printf("(shpgeo:SHPArea_2d) part %d, vtx %d \n", ring,  ring_nVertices);
#endif    
     cArea += RingArea_2d ( ring_nVertices, 
     	(double*) &(psCShape->padfX [psCShape->panPartStart[ring]]), 
     	(double*) &(psCShape->padfY [psCShape->panPartStart[ring]]) );

     ring_vtx = psCShape->panPartStart[ring];
    } 

#ifdef DEBUG2
    printf ("(shpgeo:SHPArea_2d) Area = %f \n", cArea);
#endif

    /* 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
 *
 * **************************************************************************/
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
 *
 * **************************************************************************/
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 
 *
 * **************************************************************************/
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;
#ifdef DEBUG2  
  printf("(shpgeo:RingArea) %d vertices \n", nVertices);
#endif  
  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;
#ifdef DEBUG2
    printf ("(shpgeo:RingArea)  dxArea %f  sArea %f for pt(%f, %f)\n", 
    		dx_Area, Area, x, y);
#endif
    		
    ppx = x;
    ppy = y;
  }

#ifdef DEBUG2
  printf ("(shpgeo:RingArea)  total RingArea %f \n", Area);
#endif  
  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
 *
 * **************************************************************************/
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
 */

#ifdef DEBUG2    
    printf ("(SHPUnCompound) asked for ring %d, lastring is %d \n", lRing, ring);
#endif
    return ( SHPClone(psCShape, lRing, ring ) );

}


/* **************************************************************************
 * SHPIntersect_2d
 *
 * 
 * prototype only for now
 *
 * return object with lowest common dimensionality of objects
 * 
 * **************************************************************************/ 
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						 
 * **************************************************************************/ 
int SHPClean ( SHPObject *psCShape ) {


    return (0);
}


/* **************************************************************************
 * SHPClone
 *
 * Clone a SHPObject, replicating all data
 * 
 * **************************************************************************/ 
SHPObject* SHPClone ( SHPObject *psCShape, int lowPart, int highPart ) {
    SHPObject	*psObject;
    int		newParts, newVertices;
#ifdef DEBUG
    int     	i;
#endif

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

#ifdef DEBUG
    printf (" cloning SHP (%d parts) from ring %d upto ring %d \n",
	 psCShape->nParts, lowPart, highPart);
#endif

    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]; }


#ifdef DEBUG
    if ( highPart = psCShape->nParts ) 
      i = psCShape->nVertices;
     else
      i = psCShape->panPartStart[highPart];
    printf (" from part %d (%d) to %d (%d) is %d vertices \n",
	 lowPart, psCShape->panPartStart[lowPart], highPart,
	 i, newVertices);
#endif
    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 );
}



/************************************************************************/
/*  SwapG 							                              	*/
/*                                                                      */
/*      Swap a 2, 4 or 8 byte word.                                     */
/************************************************************************/
void SwapG( void *so, void *in, int this_cnt, int this_size ) {
    int		i, j;
    unsigned char	temp;

/* return to a new pointer otherwise it would invalidate existing data	*/
/* as prevent further use of it											*/

    for( j=0; j < this_cnt; j++ )
     {
      for( i=0; i < this_size/2; i++ )
      {
	((unsigned char *) so)[i] = ((unsigned char *) in)[this_size-i-1];
	((unsigned char *) so)[this_size-i-1] = ((unsigned char *) in)[i];
      }
    }
}


/* **************************************************************************
 * SwapW
 *
 * change byte order on an array of 16 bit words
 * need to change this over to shapelib, Frank Warmerdam's functions
 *
 * **************************************************************************/ 
void swapW (void *so, unsigned char *in, long bytes) {
  int i, j;
  unsigned char map[4] = {3,2,1,0};
  unsigned char *out;

  out = so;
  for (i=0; i <= (bytes/4); i++)
   for (j=0; j < 4; j++)
      out[(i*4)+map[j]] = in[(i*4)+j]; 
}


/* **************************************************************************
 * SwapD
 *
 * change byte order on an array of (double) 32 bit words
 * need to change this over to shapelib, Frank Warmerdam's functons
 *
 * **************************************************************************/ 
void swapD (void *so, unsigned char *in, long bytes) {
  int i, j;
  unsigned char map[8] = {7,6,5,4,3,2,1,0};
  unsigned char *out;

  out = so;
  for (i=0; i <= (bytes/8); i++)
   for (j=0; j < 8; j++)
      out[(i*8)+map[j]] = in[(i*8)+j]; 
}