/******************************************************************************
* 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.16 2017-07-10 18:01:35 erouault
* * contrib/shpgeo.c: fix compilation on _MSC_VER < 1800 regarding lack
* of NAN macro.
*
* Revision 1.15 2016-12-06 21:13:33 erouault
* * configure.ac: change soname to 2:1:0 to be in sync with Debian soname.
* http://bugzilla.maptools.org/show_bug.cgi?id=2628
* Patch by Bas Couwenberg
*
* * contrib/doc/Shape_PointInPoly_README.txt, contrib/shpgeo.c: typo fixes.
* http://bugzilla.maptools.org/show_bug.cgi?id=2629
* Patch by Bas Couwenberg
*
* * web/*: use a local .css file to avoid a privacy breach issue reported
* by the lintian QA tool.
* http://bugzilla.maptools.org/show_bug.cgi?id=2630
* Patch by Bas Couwenberg
*
*
* Contributed by Sandro Mani: https://github.com/manisandro/shapelib/tree/autotools
*
* Revision 1.14 2016-12-05 12:44:07 erouault
* * Major overhaul of Makefile build system to use autoconf/automake.
*
* * Warning fixes in contrib/
*
* Revision 1.13 2011-07-24 03:17:46 fwarmerdam
* include string.h and stdlib.h where needed in contrib (#2146)
*
* Revision 1.12 2007-09-03 23:17:46 fwarmerdam
* fix SHPDimension() function
*
* Revision 1.11 2006/11/06 20:45:58 fwarmerdam
* Fixed SHPProject.
*
* Revision 1.10 2006/11/06 20:44:58 fwarmerdam
* SHPProject() uses pj_transform now
*
* Revision 1.9 2006/01/25 15:33:50 fwarmerdam
* fixed ppsC assignment maptools bug 1263
*
* Revision 1.8 2002/01/15 14:36:56 warmerda
* upgrade to use proj_api.h
*
* Revision 1.7 2002/01/11 15:22:04 warmerda
* fix many warnings. Lots of this code is cruft.
*
* Revision 1.6 2001/08/30 13:42:31 warmerda
* avoid use of auto initialization of PT for VC++
*
* 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 <stdlib.h>
#include <string.h>
#include <math.h>
#include "shapefil.h"
#include "shpgeo.h"
#if defined(_MSC_VER) && _MSC_VER < 1800
#include <float.h>
#define INFINITY (DBL_MAX + DBL_MAX)
#define NAN (INFINITY - INFINITY)
#endif
/* 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, projPJ inproj, projPJ outproj ) {
#ifdef PROJ4
int j;
if ( pj_is_latlong(inproj) ) {
for(j=0; j < psCShape->nVertices; j++) {
psCShape->padfX[j] *= DEG_TO_RAD;
psCShape->padfY[j] *= DEG_TO_RAD;
}
}
pj_transform(inproj, outproj, psCShape->nVertices, 0, psCShape->padfX,
psCShape->padfY, NULL);
if ( pj_is_latlong(outproj) ) {
for(j=0; j < psCShape->nVertices; j++) {
psCShape->padfX[j] *= RAD_TO_DEG;
psCShape->padfY[j] *= RAD_TO_DEG;
}
}
/* 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
*
* **************************************************************************/
projPJ SHPSetProjection ( int param_cnt, char **params ) {
#ifdef PROJ4
projPJ *p = NULL;
if ( param_cnt > 0 && params[0] )
{
p = pj_init ( param_cnt, params );
}
else
{
char* params_local[] = { "+proj=longlat", NULL };
p = pj_init ( 1, params_local );
}
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 ( projPJ 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;
}
return 0;
}
/* **************************************************************************
* 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 = 0;*/
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;
return 0;
}
/* **************************************************************************
* 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;
return 0;
}
/* **************************************************************************
* SHPReadOGisWKB
*
* Encapsulate entire SHPObject for use with Postgresql
*
* **************************************************************************/
SHPObject* SHPReadOGisWKB ( WKBStreamObj *stream_obj) {
SHPObject *psCShape;
char WKB_order;
int need_swap = 0, my_order, GeoType = 0;
int use_Z = 0, use_M = 0;
int nSHPType, thisDim;
WKBStreamRead ( stream_obj, &WKB_order, 1, sizeof(char));
my_order = 1;
my_order = ((char*) (&my_order))[0];
stream_obj->NeedSwap = !(WKB_order & my_order);
/* convert OGis Types to SHP types */
nSHPType = SHPOGisType ( GeoType, 0 );
WKBStreamRead ( stream_obj, &GeoType, 1, sizeof(int));
thisDim = SHPDimension ( nSHPType );
if ( thisDim && SHPD_AREA )
{ psCShape = SHPReadOGisPolygon ( stream_obj ); }
else {
if ( thisDim && SHPD_LINE )
{ psCShape = SHPReadOGisLine ( stream_obj ); }
else {
if ( thisDim && SHPD_POINT )
{ psCShape = SHPReadOGisPoint ( stream_obj ); }
}
}
return (0);
}
/* **************************************************************************
* SHPWriteOGisWKB
*
* Encapsulate entire SHPObject for use with Postgresql
*
* **************************************************************************/
int SHPWriteOGisWKB ( WKBStreamObj* stream_obj, SHPObject *psCShape ) {
int need_swap = 0, my_order, GeoType, thisDim;
int use_Z = 0, use_M = 0;
char LSB = 1;
/* indicate that this WKB is in LSB Order */
/* OGis WKB can handle either byte order, but if I get to choose I'd
/* rather have it predicatable system-to-system */
if ( stream_obj ) {
if ( stream_obj->wStream )
free ( stream_obj->wStream );
} else
{ stream_obj = calloc ( 3, sizeof (int ) ); }
/* object size needs to be 9 bytes for the wrapper, and for each polygon */
/* another 9 bytes all plus twice the total number of vertices */
/* times the sizeof (double) and just pad with 10 more chars for fun */
stream_obj->wStream = calloc (1, (9 * (psCShape->nParts + 1)) +
( sizeof(double) * 2 * psCShape->nVertices ) + 10 );
#ifdef DEBUG2
printf (" I just allocated %d bytes to wkbObj \n",
(int)(sizeof (int) + sizeof (int) + sizeof(int) +
( sizeof(int) * psCShape->nParts + 1 ) +
( sizeof(double) * 2 * psCShape->nVertices ) + 10) );
#endif
my_order = 1;
my_order = ((char*) (&my_order))[0];
/* Need to swap if this system is not LSB (Intel Order) */
stream_obj->NeedSwap = ( my_order != LSB );
stream_obj->StreamPos = 0;
#ifdef DEBUG2
printf ("this system is (%d) LSB recorded as needSwap %d\n",my_order, stream_obj->NeedSwap);
#endif
WKBStreamWrite ( stream_obj, & LSB, 1, sizeof(char) );
#ifdef DEBUG2
printf ("this system in LSB \n");
#endif
/* convert SHP Types to OGis types */
GeoType = SHPOGisType ( psCShape->nSHPType, 1 );
WKBStreamWrite ( stream_obj, &GeoType, 1, sizeof(int) );
thisDim = SHPDimension ( psCShape->nSHPType );
if ( thisDim && SHPD_AREA )
{ SHPWriteOGisPolygon ( stream_obj, psCShape ); }
else {
if ( thisDim && SHPD_LINE )
{ SHPWriteOGisLine ( stream_obj, psCShape ); }
else {
if ( thisDim && SHPD_POINT )
{ SHPWriteOGisPoint ( stream_obj, psCShape ); }
}
}
#ifdef DEBUG2
printf("(SHPWriteOGisWKB) outta here when stream pos is %d \n", stream_obj->StreamPos);
#endif
return (0);
}
/* **************************************************************************
* SHPWriteOGisPolygon
*
* for this pass code to more generic OGis MultiPolygon Type
* later add support for OGis Polygon Type
*
* Encapsulate entire SHPObject for use with Postgresql
*
* **************************************************************************/
int SHPWriteOGisPolygon ( WKBStreamObj *stream_obj, SHPObject *psCShape ) {
SHPObject **ppsC;
SHPObject *psC;
int rPart, ring, rVertices, cpart, cParts, nextring, i, j;
char Flag = 1;
int GeoType = OGIST_POLYGON;
/* cant have more than nParts complex objects in this object */
ppsC = calloc ( psCShape->nParts, sizeof(int) );
nextring = 0;
cParts=0;
while ( nextring >= 0 ) {
ppsC[cParts] = SHPUnCompound ( psCShape, &nextring );
cParts++;
}
#ifdef DEBUG2
printf ("(SHPWriteOGisPolygon) Uncompounded into %d parts \n", cParts);
#endif
WKBStreamWrite ( stream_obj, &cParts, 1, sizeof(int) );
for ( cpart = 0; cpart < cParts; cpart++) {
WKBStreamWrite ( stream_obj, & Flag, 1, sizeof(char) );
WKBStreamWrite ( stream_obj, & GeoType, 1, sizeof(int) );
psC = (SHPObject*) ppsC[cpart];
WKBStreamWrite ( stream_obj, &(psC->nParts), 1, sizeof(int) );
for ( ring = 0; (ring < (psC->nParts)) && (psC->nParts > 0); ring ++) {
if ( ring < (psC->nParts-2) )
{ rVertices = psC->panPartStart[ring+1] - psC->panPartStart[ring]; }
else
{ rVertices = psC->nVertices - psC->panPartStart[ring]; }
#ifdef DEBUG2
printf ("(SHPWriteOGisPolygon) scanning part %d, ring %d %d vtxs \n",
cpart, ring, rVertices);
#endif
rPart = psC->panPartStart[ring];
WKBStreamWrite ( stream_obj, &rVertices, 1, sizeof(int) );
for ( j=rPart; j < (rPart + rVertices); j++ ) {
WKBStreamWrite ( stream_obj, &(psC->padfX[j]), 1, sizeof(double) );
WKBStreamWrite ( stream_obj, &(psC->padfY[j]), 1, sizeof(double) );
} /* for each vertex */
} /* for each ring */
} /* for each complex part */
#ifdef DEBUG2
printf ("(SHPWriteOGisPolygon) outta here \n");
#endif
return (1);
}
/* **************************************************************************
* SHPWriteOGisLine
*
* for this pass code to more generic OGis MultiXXXXXXX Type
* later add support for OGis LineString Type
*
* Encapsulate entire SHPObject for use with Postgresql
*
* **************************************************************************/
int SHPWriteOGisLine ( WKBStreamObj *stream_obj, SHPObject *psCShape ) {
return ( SHPWriteOGisPolygon( stream_obj, psCShape ));
}
/* **************************************************************************
* SHPWriteOGisPoint
*
* for this pass code to more generic OGis MultiPoint Type
* later add support for OGis Point Type
*
* Encapsulate entire SHPObject for use with Postgresql
*
* **************************************************************************/
int SHPWriteOGisPoint ( WKBStreamObj *stream_obj, SHPObject *psCShape ) {
int j;
WKBStreamWrite ( stream_obj, &(psCShape->nVertices), 1, sizeof(int) );
for ( j=0; j < psCShape->nVertices; j++ ) {
WKBStreamWrite ( stream_obj, &(psCShape->padfX[j]), 1, sizeof(double) );
WKBStreamWrite ( stream_obj, &(psCShape->padfY[j]), 1, sizeof(double) );
} /* for each vertex */
return (1);
}
/* **************************************************************************
* SHPReadOGisPolygon
*
* for this pass code to more generic OGis MultiPolygon Type
* later add support for OGis Polygon Type
*
* Encapsulate entire SHPObject for use with Postgresql
*
* **************************************************************************/
SHPObject* SHPReadOGisPolygon ( WKBStreamObj *stream_obj ) {
SHPObject **ppsC;
SHPObject *psC;
int rPart, ring, rVertices, cpart, cParts, nextring, i, j;
int totParts, totVertices, pRings, nParts;
psC = SHPCreateObject ( SHPT_POLYGON, -1, 0, NULL, NULL, 0,
NULL, NULL, NULL, NULL );
/* initialize a blank SHPObject */
WKBStreamRead ( stream_obj, &cParts, 1, sizeof(char) );
totParts = cParts;
totVertices = 0;
SfRealloc ( psC->panPartStart, cParts * sizeof(int));
SfRealloc ( psC->panPartType, cParts * sizeof(int));
for ( cpart = 0; cpart < cParts; cpart++) {
WKBStreamRead ( stream_obj, &nParts, 1, sizeof(int) );
pRings = nParts;
/* pRings is the number of rings prior to the Ring loop below */
if ( nParts > 1 ) {
totParts += nParts - 1;
SfRealloc ( psC->panPartStart, totParts * sizeof(int));
SfRealloc ( psC->panPartType, totParts * sizeof(int));
}
rPart = 0;
for ( ring = 0; ring < (nParts - 1); ring ++) {
WKBStreamRead ( stream_obj, &rVertices, 1, sizeof(int) );
totVertices += rVertices;
psC->panPartStart[ring+pRings] = rPart;
if ( ring == 0 )
{ psC->panPartType[ring + pRings] = SHPP_OUTERRING; }
else
{ psC->panPartType[ring + pRings] = SHPP_INNERRING; }
SfRealloc ( psC->padfX, totVertices * sizeof (double));
SfRealloc ( psC->padfY, totVertices * sizeof (double));
for ( j=rPart; j < (rPart + rVertices); j++ ) {
WKBStreamRead ( stream_obj, &(psC->padfX[j]), 1, sizeof(double) );
WKBStreamRead ( stream_obj, &(psC->padfY[j]), 1, sizeof(double) );
} /* for each vertex */
rPart += rVertices;
} /* for each ring */
} /* for each complex part */
return ( psC );
}
/* **************************************************************************
* SHPReadOGisLine
*
* for this pass code to more generic OGis MultiLineString Type
* later add support for OGis LineString Type
*
* Encapsulate entire SHPObject for use with Postgresql
*
* **************************************************************************/
SHPObject* SHPReadOGisLine ( WKBStreamObj *stream_obj ) {
SHPObject **ppsC;
SHPObject *psC;
int rPart, ring, rVertices, cpart, cParts, nextring, i, j;
int totParts, totVertices, pRings, nParts;
psC = SHPCreateObject ( SHPT_ARC, -1, 0, NULL, NULL, 0,
NULL, NULL, NULL, NULL );
/* initialize a blank SHPObject */
WKBStreamRead ( stream_obj, &cParts, 1, sizeof(int) );
totParts = cParts;
totVertices = 0;
SfRealloc ( psC->panPartStart, cParts * sizeof(int));
SfRealloc ( psC->panPartType, cParts * sizeof(int));
for ( cpart = 0; cpart < cParts; cpart++) {
WKBStreamRead ( stream_obj, &nParts, 1, sizeof(int) );
pRings = totParts;
/* pRings is the number of rings prior to the Ring loop below */
if ( nParts > 1 ) {
totParts += nParts - 1;
SfRealloc ( psC->panPartStart, totParts * sizeof(int));
SfRealloc ( psC->panPartType, totParts * sizeof(int));
}
rPart = 0;
for ( ring = 0; ring < (nParts - 1); ring ++) {
WKBStreamRead ( stream_obj, &rVertices, 1, sizeof(int) );
totVertices += rVertices;
psC->panPartStart[ring+pRings] = rPart;
if ( ring == 0 )
{ psC->panPartType[ring + pRings] = SHPP_OUTERRING; }
else
{ psC->panPartType[ring + pRings] = SHPP_INNERRING; }
SfRealloc ( psC->padfX, totVertices * sizeof (double));
SfRealloc ( psC->padfY, totVertices * sizeof (double));
for ( j=rPart; j < (rPart + rVertices); j++ ) {
WKBStreamRead ( stream_obj, &(psC->padfX[j]), 1, sizeof(double) );
WKBStreamRead ( stream_obj, &(psC->padfY[j]), 1, sizeof(double) );
} /* for each vertex */
rPart += rVertices;
} /* for each ring */
} /* for each complex part */
return ( psC );
}
/* **************************************************************************
* SHPReadOGisPoint
*
* Encapsulate entire SHPObject for use with Postgresql
*
* **************************************************************************/
SHPObject* SHPReadOGisPoint ( WKBStreamObj *stream_obj ) {
SHPObject *psC;
int nVertices, j;
psC = SHPCreateObject ( SHPT_MULTIPOINT, -1, 0, NULL, NULL, 0,
NULL, NULL, NULL, NULL );
/* initialize a blank SHPObject */
WKBStreamRead ( stream_obj, &nVertices, 1, sizeof(int) );
SfRealloc ( psC->padfX, nVertices * sizeof (double));
SfRealloc ( psC->padfY, nVertices * sizeof (double));
for ( j=0; j < nVertices; j++ ) {
WKBStreamRead ( stream_obj, &(psC->padfX[j]), 1, sizeof(double) );
WKBStreamRead ( stream_obj, &(psC->padfY[j]), 1, sizeof(double) );
} /* for each vertex */
return ( psC );
}
/* **************************************************************************
* RingReadOGisWKB
*
* this accepts OGisLineStrings which are basic building blocks
*
* Encapsulate entire SHPObject for use with Postgresql
*
* **************************************************************************/
int RingReadOgisWKB ( SHPObject *psCShape, char *stream_obj) {
return 0;
}
/* **************************************************************************
* RingWriteOGisWKB
*
* this emits OGisLineStrings which are basic building blocks
*
* Encapsulate entire SHPObject for use with Postgresql
*
* **************************************************************************/
int RingWriteOgisWKB ( SHPObject *psCShape, char *stream_obj) {
return 0;
}
/* **************************************************************************
* 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; break;
case SHPT_ARC : dimension = SHPD_LINE; break;
case SHPT_POLYGON : dimension = SHPD_AREA; break;
case SHPT_MULTIPOINT : dimension = SHPD_POINT; break;
case SHPT_POINTZ : dimension = SHPD_POINT | SHPD_Z; break;
case SHPT_ARCZ : dimension = SHPD_LINE | SHPD_Z; break;
case SHPT_POLYGONZ : dimension = SHPD_AREA | SHPD_Z; break;
case SHPT_MULTIPOINTZ : dimension = SHPD_POINT | SHPD_Z; break;
case SHPT_POINTM : dimension = SHPD_POINT | SHPD_MEASURE; break;
case SHPT_ARCM : dimension = SHPD_LINE | SHPD_MEASURE; break;
case SHPT_POLYGONM : dimension = SHPD_AREA | SHPD_MEASURE; break;
case SHPT_MULTIPOINTM : dimension = SHPD_POINT | SHPD_MEASURE; break;
case SHPT_MULTIPATCH : dimension = SHPD_AREA; break;
}
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) )
{
rPT.x = NAN;
rPT.y = NAN;
return rPT;
}
sPT = SHPPointsinPoly_2d ( psCShape );
if ( sPT ) {
rPT.x = sPT[0].x;
rPT.y = sPT[0].y;
} else {
rPT.x = NAN;
rPT.y = 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 = 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 = 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) )
{
C.x = NAN;
C.y = NAN;
return C;
}
#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 */
/* 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 to 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];
}