shpgeo.c

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

/******************************************************************************
 * 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.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 "shapefil.h"
#ifdef PROJ4
  #define projUV UV
  #include <projects.h>
#else
  #define PJ void*
#endif

#include "my_nan.h" 
#include "shpgeo.h"


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


/* **************************************************************************
 * SHPFreeProjection
 *
 * release a projection handle for use with PROJ4.3
 *
 * act as a wrapper to protect against library changes in PROJ
 * 
 * **************************************************************************/
int SHPFreeProjection ( PJ *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;
	  } 	  	  	  

}


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



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