shpgeo.cpp

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

#include "stdafx.h"

#define PROJ4
#include "shapefil.h"
#ifdef PROJ4
  #include "projects.h"
#else
  #define PJ void*
#endif

#include "shpgeo.h"
#include <math.h>

typedef projUV UV; 

/* **************************************************************************
 * asFileName
 *
 * utility function, toss part of filename after last dot
 *
 * **************************************************************************/ 
PROJAPI_CALL 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
 * 
 * **************************************************************************/ 
PROJAPI_CALL int SHPProject ( SHPObject *psCShape, PJ *inproj, PJ *outproj ) {
#ifdef	PROJ4
   int	j;
   UV   p;    /* struct { double u, double v } */

   /* for each vertex project it and stuff the projeted point back into 	*/
   /*	same SHPObject.  Proj assumes data is in radians so convert it.		*/
   /*   Proj will convert Geographic -> <proj> and <proj> -> Geographic		*/
   /*	so <proj1> -> <proj2> requires bouncing though geographic			*/
 /*  
   for ( j=0; j < psCShape->nVertices; j++ ) {
     p.u = psCShape->padfX[j];
     p.v = psCShape->padfY[j];

     if ( inproj )
       p = pj_inv ( p, inproj );
     else
      { p.u *= DEG_TO_RAD;
        p.v *= DEG_TO_RAD;
      }

     if ( outproj ) 
       p = pj_fwd ( p, outproj );
     else
     { p.u *= RAD_TO_DEG;
       p.v *= RAD_TO_DEG;
     } 

     psCShape->padfX[j] = p.u;
     psCShape->padfY[j] = p.v;
   }
  */
	int	error;

	for ( j=0; j < psCShape->nVertices; j++ ) 
	{
		p.u = psCShape->padfX[j];
		p.v = psCShape->padfY[j];

		if( inproj && outproj )
		{
			double	z = 0.0;

			if( pj_is_latlong(inproj) )
			{
			  p.u *= DEG_TO_RAD;
			  p.v *= DEG_TO_RAD;
			}

			error = pj_transform( inproj, outproj, 1, 0, 
								&p.u, &p.v, &z );

			if( error )
			  return 0;

			if( pj_is_latlong(outproj) )
			{
			  p.u *= RAD_TO_DEG;
			  p.v *= RAD_TO_DEG;
			}
		}
		else
		{
			if(inproj == NULL) /* input coordinates are lat/lon */
			{ 
			  p.u *= DEG_TO_RAD; /* convert to radians */
			  p.v *= DEG_TO_RAD;  
			  p = pj_fwd(p, outproj);
			} 
			else 
			{
			  if(outproj == NULL) /* output coordinates are lat/lon */
			  { 
				  p = pj_inv(p, inproj);
				  p.u *= RAD_TO_DEG; /* convert to decimal degrees */
				  p.v *= RAD_TO_DEG;
			  } 
			  else 
			  { /* need to go from one projection to another */
				  p = pj_inv(p, inproj);
				  p = pj_fwd(p, outproj);
			  }
			}
		}

		psCShape->padfX[j] = p.u;
		psCShape->padfY[j] = p.v;
	}

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

	return ( 1 );
}

PROJAPI_CALL int SHPProjectPoint ( double *x, double *y, PJ *inproj, PJ *outproj )
{
#ifdef	PROJ4
   UV   p;    /* struct { double u, double v } */
	int	error;

	p.u = *x;
	p.v = *y;

	if( inproj && outproj )
	{
		double	z = 0.0;

		if( pj_is_latlong(inproj) )
		{
		  p.u *= DEG_TO_RAD;
		  p.v *= DEG_TO_RAD;
		}

		error = pj_transform( inproj, outproj, 1, 0, 
							&p.u, &p.v, &z );
		if( error )
			return 0;

		if( pj_is_latlong(outproj) )
		{
		  p.u *= RAD_TO_DEG;
		  p.v *= RAD_TO_DEG;
		}
	}
	else
	{
		if(inproj == NULL) /* input coordinates are lat/lon */
		{ 
		  p.u *= DEG_TO_RAD; /* convert to radians */
		  p.v *= DEG_TO_RAD;  
		  p = pj_fwd(p, outproj);
		} 
		else 
		{
		  if(outproj == NULL) /* output coordinates are lat/lon */
		  { 
			  p = pj_inv(p, inproj);
			  p.u *= RAD_TO_DEG; /* convert to decimal degrees */
			  p.v *= RAD_TO_DEG;
		  } 
		  else 
		  { /* need to go from one projection to another */
			  p = pj_inv(p, inproj);
			  p = pj_fwd(p, outproj);
		  }
		}
	}

	*x = p.u;
	*y = p.v;
#endif  

	return ( 1 );
}

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

  if ( param_cnt > 0 && params[0] )
  { p = pj_init ( param_cnt, params ); }

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

//such as "+proj=utm +zone=11 +ellps=WGS84"
PROJAPI_CALL PJ *SHPSetProjectionEx (const char *definition )
{
	if(definition)
		return pj_init_plus(definition);
	else 
		return NULL;
}


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


/* **************************************************************************
 * SHPDimension
 *
 * Return the Dimensionality of the SHPObject 
 * a handy utility function
 * 
 * **************************************************************************/
PROJAPI_CALL 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
 * 
 * **************************************************************************/
PROJAPI_CALL PT	SHPPointinPoly_2d ( SHPObject *psCShape ) 
{
	PT	*sPT, rPT;

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

	sPT = SHPPointsinPoly_2d ( psCShape );
   
	if ( sPT ) 
	{
		rPT.x = sPT[0].x;
		rPT.y = sPT[0].y; 
	} 
	else 
	{
		rPT =  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
 * 
 * **************************************************************************/
PROJAPI_CALL 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 = (double *)calloc ( 4, sizeof(double));
     CLy = (double *)calloc ( 4, sizeof(double));
     CLst = (int *)calloc ( 2, sizeof(int));
     CLstt = (int *)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
 * 
 * **************************************************************************/
PROJAPI_CALL PT SHPCentrd_2d ( SHPObject *psCShape ) 
{
    int		ring, ringPrev, ring_nVertices, rStart;
    double	Area, ringArea;
    PT		ringCentrd, C;
    
  
   if ( !(SHPDimension (psCShape->nSHPType) & SHPD_AREA) )  
       return NAN;

   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;