1 /******************************************************************************
2 * $Id: ogrpgeogeometry.cpp 27959 2014-11-14 18:29:21Z rouault $
3 *
4 * Project: OpenGIS Simple Features Reference Implementation
5 * Purpose: Implements decoder of shapebin geometry for PGeo
6 * Author: Frank Warmerdam, warmerdam@pobox.com
7 * Paul Ramsey, pramsey at cleverelephant.ca
8 *
9 ******************************************************************************
10 * Copyright (c) 2005, Frank Warmerdam <warmerdam@pobox.com>
11 * Copyright (c) 2011, Paul Ramsey <pramsey at cleverelephant.ca>
12 * Copyright (c) 2011-2014, Even Rouault <even dot rouault at mines-paris dot org>
13 *
14 * Permission is hereby granted, free of charge, to any person obtaining a
15 * copy of this software and associated documentation files (the "Software"),
16 * to deal in the Software without restriction, including without limitation
17 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 * and/or sell copies of the Software, and to permit persons to whom the
19 * Software is furnished to do so, subject to the following conditions:
20 *
21 * The above copyright notice and this permission notice shall be included
22 * in all copies or substantial portions of the Software.
23 *
24 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 * DEALINGS IN THE SOFTWARE.
31 ****************************************************************************/
32
33 #include "ogrpgeogeometry.h"
34 #include "ogr_p.h"
35 #include "cpl_string.h"
36
37 CPL_CVSID("$Id: ogrpgeogeometry.cpp 27959 2014-11-14 18:29:21Z rouault $");
38
39 #define SHPP_TRISTRIP 0
40 #define SHPP_TRIFAN 1
41 #define SHPP_OUTERRING 2
42 #define SHPP_INNERRING 3
43 #define SHPP_FIRSTRING 4
44 #define SHPP_RING 5
45 #define SHPP_TRIANGLES 6 /* Multipatch 9.0 specific */
46
47
48 /************************************************************************/
49 /* OGRCreateFromMultiPatchPart() */
50 /************************************************************************/
51
OGRCreateFromMultiPatchPart(OGRMultiPolygon * poMP,OGRPolygon * & poLastPoly,int nPartType,int nPartPoints,double * padfX,double * padfY,double * padfZ)52 void OGRCreateFromMultiPatchPart(OGRMultiPolygon *poMP,
53 OGRPolygon*& poLastPoly,
54 int nPartType,
55 int nPartPoints,
56 double* padfX,
57 double* padfY,
58 double* padfZ)
59 {
60 nPartType &= 0xf;
61
62 if( nPartType == SHPP_TRISTRIP )
63 {
64 int iBaseVert;
65
66 if( poLastPoly != NULL )
67 {
68 poMP->addGeometryDirectly( poLastPoly );
69 poLastPoly = NULL;
70 }
71
72 for( iBaseVert = 0; iBaseVert < nPartPoints-2; iBaseVert++ )
73 {
74 OGRPolygon *poPoly = new OGRPolygon();
75 OGRLinearRing *poRing = new OGRLinearRing();
76 int iSrcVert = iBaseVert;
77
78 poRing->setPoint( 0,
79 padfX[iSrcVert],
80 padfY[iSrcVert],
81 padfZ[iSrcVert] );
82 poRing->setPoint( 1,
83 padfX[iSrcVert+1],
84 padfY[iSrcVert+1],
85 padfZ[iSrcVert+1] );
86
87 poRing->setPoint( 2,
88 padfX[iSrcVert+2],
89 padfY[iSrcVert+2],
90 padfZ[iSrcVert+2] );
91 poRing->setPoint( 3,
92 padfX[iSrcVert],
93 padfY[iSrcVert],
94 padfZ[iSrcVert] );
95
96 poPoly->addRingDirectly( poRing );
97 poMP->addGeometryDirectly( poPoly );
98 }
99 }
100 else if( nPartType == SHPP_TRIFAN )
101 {
102 int iBaseVert;
103
104 if( poLastPoly != NULL )
105 {
106 poMP->addGeometryDirectly( poLastPoly );
107 poLastPoly = NULL;
108 }
109
110 for( iBaseVert = 0; iBaseVert < nPartPoints-2; iBaseVert++ )
111 {
112 OGRPolygon *poPoly = new OGRPolygon();
113 OGRLinearRing *poRing = new OGRLinearRing();
114 int iSrcVert = iBaseVert;
115
116 poRing->setPoint( 0,
117 padfX[0],
118 padfY[0],
119 padfZ[0] );
120 poRing->setPoint( 1,
121 padfX[iSrcVert+1],
122 padfY[iSrcVert+1],
123 padfZ[iSrcVert+1] );
124
125 poRing->setPoint( 2,
126 padfX[iSrcVert+2],
127 padfY[iSrcVert+2],
128 padfZ[iSrcVert+2] );
129 poRing->setPoint( 3,
130 padfX[0],
131 padfY[0],
132 padfZ[0] );
133
134 poPoly->addRingDirectly( poRing );
135 poMP->addGeometryDirectly( poPoly );
136 }
137 }
138 else if( nPartType == SHPP_OUTERRING
139 || nPartType == SHPP_INNERRING
140 || nPartType == SHPP_FIRSTRING
141 || nPartType == SHPP_RING )
142 {
143 if( poLastPoly != NULL
144 && (nPartType == SHPP_OUTERRING
145 || nPartType == SHPP_FIRSTRING) )
146 {
147 poMP->addGeometryDirectly( poLastPoly );
148 poLastPoly = NULL;
149 }
150
151 if( poLastPoly == NULL )
152 poLastPoly = new OGRPolygon();
153
154 OGRLinearRing *poRing = new OGRLinearRing;
155
156 poRing->setPoints( nPartPoints,
157 padfX,
158 padfY,
159 padfZ );
160
161 poRing->closeRings();
162
163 poLastPoly->addRingDirectly( poRing );
164 }
165 else if ( nPartType == SHPP_TRIANGLES )
166 {
167 int iBaseVert;
168
169 if( poLastPoly != NULL )
170 {
171 poMP->addGeometryDirectly( poLastPoly );
172 poLastPoly = NULL;
173 }
174
175 for( iBaseVert = 0; iBaseVert < nPartPoints-2; iBaseVert+=3 )
176 {
177 OGRPolygon *poPoly = new OGRPolygon();
178 OGRLinearRing *poRing = new OGRLinearRing();
179 int iSrcVert = iBaseVert;
180
181 poRing->setPoint( 0,
182 padfX[iSrcVert],
183 padfY[iSrcVert],
184 padfZ[iSrcVert] );
185 poRing->setPoint( 1,
186 padfX[iSrcVert+1],
187 padfY[iSrcVert+1],
188 padfZ[iSrcVert+1] );
189
190 poRing->setPoint( 2,
191 padfX[iSrcVert+2],
192 padfY[iSrcVert+2],
193 padfZ[iSrcVert+2] );
194 poRing->setPoint( 3,
195 padfX[iSrcVert],
196 padfY[iSrcVert],
197 padfZ[iSrcVert] );
198
199 poPoly->addRingDirectly( poRing );
200 poMP->addGeometryDirectly( poPoly );
201 }
202 }
203 else
204 CPLDebug( "OGR", "Unrecognised parttype %d, ignored.",
205 nPartType );
206 }
207
208 /************************************************************************/
209 /* OGRCreateFromMultiPatch() */
210 /* */
211 /* Translate a multipatch representation to an OGR geometry */
212 /* Mostly copied from shape2ogr.cpp */
213 /************************************************************************/
214
OGRCreateFromMultiPatch(int nParts,GInt32 * panPartStart,GInt32 * panPartType,int nPoints,double * padfX,double * padfY,double * padfZ)215 static OGRGeometry* OGRCreateFromMultiPatch(int nParts,
216 GInt32* panPartStart,
217 GInt32* panPartType,
218 int nPoints,
219 double* padfX,
220 double* padfY,
221 double* padfZ)
222 {
223 OGRMultiPolygon *poMP = new OGRMultiPolygon();
224 int iPart;
225 OGRPolygon *poLastPoly = NULL;
226
227 for( iPart = 0; iPart < nParts; iPart++ )
228 {
229 int nPartPoints, nPartStart;
230
231 // Figure out details about this part's vertex list.
232 if( panPartStart == NULL )
233 {
234 nPartPoints = nPoints;
235 nPartStart = 0;
236 }
237 else
238 {
239
240 if( iPart == nParts - 1 )
241 nPartPoints =
242 nPoints - panPartStart[iPart];
243 else
244 nPartPoints = panPartStart[iPart+1]
245 - panPartStart[iPart];
246 nPartStart = panPartStart[iPart];
247 }
248
249 OGRCreateFromMultiPatchPart(poMP,
250 poLastPoly,
251 panPartType[iPart],
252 nPartPoints,
253 padfX + nPartStart,
254 padfY + nPartStart,
255 padfZ + nPartStart);
256 }
257
258 if( poLastPoly != NULL )
259 {
260 poMP->addGeometryDirectly( poLastPoly );
261 poLastPoly = NULL;
262 }
263
264 return poMP;
265 }
266
267
268 /************************************************************************/
269 /* OGRWriteToShapeBin() */
270 /* */
271 /* Translate OGR geometry to a shapefile binary representation */
272 /************************************************************************/
273
OGRWriteToShapeBin(OGRGeometry * poGeom,GByte ** ppabyShape,int * pnBytes)274 OGRErr OGRWriteToShapeBin( OGRGeometry *poGeom,
275 GByte **ppabyShape,
276 int *pnBytes )
277 {
278 GUInt32 nGType = SHPT_NULL;
279 int nShpSize = 4; /* All types start with integer type number */
280 int nShpZSize = 0; /* Z gets tacked onto the end */
281 GUInt32 nPoints = 0;
282 GUInt32 nParts = 0;
283
284 /* -------------------------------------------------------------------- */
285 /* Null or Empty input maps to SHPT_NULL. */
286 /* -------------------------------------------------------------------- */
287 if ( ! poGeom || poGeom->IsEmpty() )
288 {
289 *ppabyShape = (GByte*)VSIMalloc(nShpSize);
290 GUInt32 zero = SHPT_NULL;
291 memcpy(*ppabyShape, &zero, nShpSize);
292 *pnBytes = nShpSize;
293 return OGRERR_NONE;
294 }
295
296 OGRwkbGeometryType nOGRType = wkbFlatten(poGeom->getGeometryType());
297 int b3d = wkbHasZ(poGeom->getGeometryType());
298 int nCoordDims = b3d ? 3 : 2;
299
300 /* -------------------------------------------------------------------- */
301 /* Calculate the shape buffer size */
302 /* -------------------------------------------------------------------- */
303 if ( nOGRType == wkbPoint )
304 {
305 nShpSize += 8 * nCoordDims;
306 }
307 else if ( nOGRType == wkbLineString )
308 {
309 OGRLineString *poLine = (OGRLineString*)poGeom;
310 nPoints = poLine->getNumPoints();
311 nParts = 1;
312 nShpSize += 16 * nCoordDims; /* xy(z) box */
313 nShpSize += 4; /* nparts */
314 nShpSize += 4; /* npoints */
315 nShpSize += 4; /* parts[1] */
316 nShpSize += 8 * nCoordDims * nPoints; /* points */
317 nShpZSize = 16 + 8 * nPoints;
318 }
319 else if ( nOGRType == wkbPolygon )
320 {
321 poGeom->closeRings();
322 OGRPolygon *poPoly = (OGRPolygon*)poGeom;
323 nParts = poPoly->getNumInteriorRings() + 1;
324 for ( GUInt32 i = 0; i < nParts; i++ )
325 {
326 OGRLinearRing *poRing;
327 if ( i == 0 )
328 poRing = poPoly->getExteriorRing();
329 else
330 poRing = poPoly->getInteriorRing(i-1);
331 nPoints += poRing->getNumPoints();
332 }
333 nShpSize += 16 * nCoordDims; /* xy(z) box */
334 nShpSize += 4; /* nparts */
335 nShpSize += 4; /* npoints */
336 nShpSize += 4 * nParts; /* parts[nparts] */
337 nShpSize += 8 * nCoordDims * nPoints; /* points */
338 nShpZSize = 16 + 8 * nPoints;
339 }
340 else if ( nOGRType == wkbMultiPoint )
341 {
342 OGRMultiPoint *poMPoint = (OGRMultiPoint*)poGeom;
343 for ( int i = 0; i < poMPoint->getNumGeometries(); i++ )
344 {
345 OGRPoint *poPoint = (OGRPoint*)(poMPoint->getGeometryRef(i));
346 if ( poPoint->IsEmpty() )
347 continue;
348 nPoints++;
349 }
350 nShpSize += 16 * nCoordDims; /* xy(z) box */
351 nShpSize += 4; /* npoints */
352 nShpSize += 8 * nCoordDims * nPoints; /* points */
353 nShpZSize = 16 + 8 * nPoints;
354 }
355 else if ( nOGRType == wkbMultiLineString )
356 {
357 OGRMultiLineString *poMLine = (OGRMultiLineString*)poGeom;
358 for ( int i = 0; i < poMLine->getNumGeometries(); i++ )
359 {
360 OGRLineString *poLine = (OGRLineString*)(poMLine->getGeometryRef(i));
361 /* Skip empties */
362 if ( poLine->IsEmpty() )
363 continue;
364 nParts++;
365 nPoints += poLine->getNumPoints();
366 }
367 nShpSize += 16 * nCoordDims; /* xy(z) box */
368 nShpSize += 4; /* nparts */
369 nShpSize += 4; /* npoints */
370 nShpSize += 4 * nParts; /* parts[nparts] */
371 nShpSize += 8 * nCoordDims * nPoints ; /* points */
372 nShpZSize = 16 + 8 * nPoints;
373 }
374 else if ( nOGRType == wkbMultiPolygon )
375 {
376 poGeom->closeRings();
377 OGRMultiPolygon *poMPoly = (OGRMultiPolygon*)poGeom;
378 for( int j = 0; j < poMPoly->getNumGeometries(); j++ )
379 {
380 OGRPolygon *poPoly = (OGRPolygon*)(poMPoly->getGeometryRef(j));
381 int nRings = poPoly->getNumInteriorRings() + 1;
382
383 /* Skip empties */
384 if ( poPoly->IsEmpty() )
385 continue;
386
387 nParts += nRings;
388 for ( int i = 0; i < nRings; i++ )
389 {
390 OGRLinearRing *poRing;
391 if ( i == 0 )
392 poRing = poPoly->getExteriorRing();
393 else
394 poRing = poPoly->getInteriorRing(i-1);
395 nPoints += poRing->getNumPoints();
396 }
397 }
398 nShpSize += 16 * nCoordDims; /* xy(z) box */
399 nShpSize += 4; /* nparts */
400 nShpSize += 4; /* npoints */
401 nShpSize += 4 * nParts; /* parts[nparts] */
402 nShpSize += 8 * nCoordDims * nPoints ; /* points */
403 nShpZSize = 16 + 8 * nPoints;
404 }
405 else
406 {
407 return OGRERR_UNSUPPORTED_OPERATION;
408 }
409
410 /* Allocate our shape buffer */
411 *ppabyShape = (GByte*)VSIMalloc(nShpSize);
412 if ( ! *ppabyShape )
413 return OGRERR_FAILURE;
414
415 /* Fill in the output size. */
416 *pnBytes = nShpSize;
417
418 /* Set up write pointers */
419 unsigned char *pabyPtr = *ppabyShape;
420 unsigned char *pabyPtrZ = NULL;
421 if ( b3d )
422 pabyPtrZ = pabyPtr + nShpSize - nShpZSize;
423
424 /* -------------------------------------------------------------------- */
425 /* Write in the Shape type number now */
426 /* -------------------------------------------------------------------- */
427 switch(nOGRType)
428 {
429 case wkbPoint:
430 {
431 nGType = b3d ? SHPT_POINTZ : SHPT_POINT;
432 break;
433 }
434 case wkbMultiPoint:
435 {
436 nGType = b3d ? SHPT_MULTIPOINTZ : SHPT_MULTIPOINT;
437 break;
438 }
439 case wkbLineString:
440 case wkbMultiLineString:
441 {
442 nGType = b3d ? SHPT_ARCZ : SHPT_ARC;
443 break;
444 }
445 case wkbPolygon:
446 case wkbMultiPolygon:
447 {
448 nGType = b3d ? SHPT_POLYGONZ : SHPT_POLYGON;
449 break;
450 }
451 default:
452 {
453 return OGRERR_UNSUPPORTED_OPERATION;
454 }
455 }
456 /* Write in the type number and advance the pointer */
457 nGType = CPL_LSBWORD32( nGType );
458 memcpy( pabyPtr, &nGType, 4 );
459 pabyPtr += 4;
460
461
462 /* -------------------------------------------------------------------- */
463 /* POINT and POINTZ */
464 /* -------------------------------------------------------------------- */
465 if ( nOGRType == wkbPoint )
466 {
467 OGRPoint *poPoint = (OGRPoint*)poGeom;
468 double x = poPoint->getX();
469 double y = poPoint->getY();
470
471 /* Copy in the raw data. */
472 memcpy( pabyPtr, &x, 8 );
473 memcpy( pabyPtr+8, &y, 8 );
474 if( b3d )
475 {
476 double z = poPoint->getZ();
477 memcpy( pabyPtr+8+8, &z, 8 );
478 }
479
480 /* Swap if needed. Shape doubles always LSB */
481 if( OGR_SWAP( wkbNDR ) )
482 {
483 CPL_SWAPDOUBLE( pabyPtr );
484 CPL_SWAPDOUBLE( pabyPtr+8 );
485 if( b3d )
486 CPL_SWAPDOUBLE( pabyPtr+8+8 );
487 }
488
489 return OGRERR_NONE;
490 }
491
492 /* -------------------------------------------------------------------- */
493 /* All the non-POINT types require an envelope next */
494 /* -------------------------------------------------------------------- */
495 OGREnvelope3D envelope;
496 poGeom->getEnvelope(&envelope);
497 memcpy( pabyPtr, &(envelope.MinX), 8 );
498 memcpy( pabyPtr+8, &(envelope.MinY), 8 );
499 memcpy( pabyPtr+8+8, &(envelope.MaxX), 8 );
500 memcpy( pabyPtr+8+8+8, &(envelope.MaxY), 8 );
501
502 /* Swap box if needed. Shape doubles are always LSB */
503 if( OGR_SWAP( wkbNDR ) )
504 {
505 for ( int i = 0; i < 4; i++ )
506 CPL_SWAPDOUBLE( pabyPtr + 8*i );
507 }
508 pabyPtr += 32;
509
510 /* Write in the Z bounds at the end of the XY buffer */
511 if ( b3d )
512 {
513 memcpy( pabyPtrZ, &(envelope.MinZ), 8 );
514 memcpy( pabyPtrZ+8, &(envelope.MaxZ), 8 );
515
516 /* Swap Z bounds if necessary */
517 if( OGR_SWAP( wkbNDR ) )
518 {
519 for ( int i = 0; i < 2; i++ )
520 CPL_SWAPDOUBLE( pabyPtrZ + 8*i );
521 }
522 pabyPtrZ += 16;
523 }
524
525 /* -------------------------------------------------------------------- */
526 /* LINESTRING and LINESTRINGZ */
527 /* -------------------------------------------------------------------- */
528 if ( nOGRType == wkbLineString )
529 {
530 const OGRLineString *poLine = (OGRLineString*)poGeom;
531
532 /* Write in the nparts (1) */
533 GUInt32 nPartsLsb = CPL_LSBWORD32( nParts );
534 memcpy( pabyPtr, &nPartsLsb, 4 );
535 pabyPtr += 4;
536
537 /* Write in the npoints */
538 GUInt32 nPointsLsb = CPL_LSBWORD32( nPoints );
539 memcpy( pabyPtr, &nPointsLsb, 4 );
540 pabyPtr += 4;
541
542 /* Write in the part index (0) */
543 GUInt32 nPartIndex = 0;
544 memcpy( pabyPtr, &nPartIndex, 4 );
545 pabyPtr += 4;
546
547 /* Write in the point data */
548 poLine->getPoints((OGRRawPoint*)pabyPtr, (double*)pabyPtrZ);
549
550 /* Swap if necessary */
551 if( OGR_SWAP( wkbNDR ) )
552 {
553 for( GUInt32 k = 0; k < nPoints; k++ )
554 {
555 CPL_SWAPDOUBLE( pabyPtr + 16*k );
556 CPL_SWAPDOUBLE( pabyPtr + 16*k + 8 );
557 if( b3d )
558 CPL_SWAPDOUBLE( pabyPtrZ + 8*k );
559 }
560 }
561 return OGRERR_NONE;
562
563 }
564 /* -------------------------------------------------------------------- */
565 /* POLYGON and POLYGONZ */
566 /* -------------------------------------------------------------------- */
567 else if ( nOGRType == wkbPolygon )
568 {
569 OGRPolygon *poPoly = (OGRPolygon*)poGeom;
570
571 /* Write in the part count */
572 GUInt32 nPartsLsb = CPL_LSBWORD32( nParts );
573 memcpy( pabyPtr, &nPartsLsb, 4 );
574 pabyPtr += 4;
575
576 /* Write in the total point count */
577 GUInt32 nPointsLsb = CPL_LSBWORD32( nPoints );
578 memcpy( pabyPtr, &nPointsLsb, 4 );
579 pabyPtr += 4;
580
581 /* -------------------------------------------------------------------- */
582 /* Now we have to visit each ring and write an index number into */
583 /* the parts list, and the coordinates into the points list. */
584 /* to do it in one pass, we will use three write pointers. */
585 /* pabyPtr writes the part indexes */
586 /* pabyPoints writes the xy coordinates */
587 /* pabyPtrZ writes the z coordinates */
588 /* -------------------------------------------------------------------- */
589
590 /* Just past the partindex[nparts] array */
591 unsigned char* pabyPoints = pabyPtr + 4*nParts;
592
593 int nPointIndexCount = 0;
594
595 for( GUInt32 i = 0; i < nParts; i++ )
596 {
597 /* Check our Ring and condition it */
598 OGRLinearRing *poRing;
599 if ( i == 0 )
600 {
601 poRing = poPoly->getExteriorRing();
602 /* Outer ring must be clockwise */
603 if ( ! poRing->isClockwise() )
604 poRing->reverseWindingOrder();
605 }
606 else
607 {
608 poRing = poPoly->getInteriorRing(i-1);
609 /* Inner rings should be anti-clockwise */
610 if ( poRing->isClockwise() )
611 poRing->reverseWindingOrder();
612 }
613
614 int nRingNumPoints = poRing->getNumPoints();
615
616 /* Cannot write un-closed rings to shape */
617 if( nRingNumPoints <= 2 || ! poRing->get_IsClosed() )
618 return OGRERR_FAILURE;
619
620 /* Write in the part index */
621 GUInt32 nPartIndex = CPL_LSBWORD32( nPointIndexCount );
622 memcpy( pabyPtr, &nPartIndex, 4 );
623
624 /* Write in the point data */
625 poRing->getPoints((OGRRawPoint*)pabyPoints, (double*)pabyPtrZ);
626
627 /* Swap if necessary */
628 if( OGR_SWAP( wkbNDR ) )
629 {
630 for( int k = 0; k < nRingNumPoints; k++ )
631 {
632 CPL_SWAPDOUBLE( pabyPoints + 16*k );
633 CPL_SWAPDOUBLE( pabyPoints + 16*k + 8 );
634 if( b3d )
635 CPL_SWAPDOUBLE( pabyPtrZ + 8*k );
636 }
637 }
638
639 nPointIndexCount += nRingNumPoints;
640 /* Advance the write pointers */
641 pabyPtr += 4;
642 pabyPoints += 16 * nRingNumPoints;
643 if ( b3d )
644 pabyPtrZ += 8 * nRingNumPoints;
645 }
646
647 return OGRERR_NONE;
648
649 }
650 /* -------------------------------------------------------------------- */
651 /* MULTIPOINT and MULTIPOINTZ */
652 /* -------------------------------------------------------------------- */
653 else if ( nOGRType == wkbMultiPoint )
654 {
655 OGRMultiPoint *poMPoint = (OGRMultiPoint*)poGeom;
656
657 /* Write in the total point count */
658 GUInt32 nPointsLsb = CPL_LSBWORD32( nPoints );
659 memcpy( pabyPtr, &nPointsLsb, 4 );
660 pabyPtr += 4;
661
662 /* -------------------------------------------------------------------- */
663 /* Now we have to visit each point write it into the points list */
664 /* We will use two write pointers. */
665 /* pabyPtr writes the xy coordinates */
666 /* pabyPtrZ writes the z coordinates */
667 /* -------------------------------------------------------------------- */
668
669 for( GUInt32 i = 0; i < nPoints; i++ )
670 {
671 const OGRPoint *poPt = (OGRPoint*)(poMPoint->getGeometryRef(i));
672
673 /* Skip empties */
674 if ( poPt->IsEmpty() )
675 continue;
676
677 /* Write the coordinates */
678 double x = poPt->getX();
679 double y = poPt->getY();
680 memcpy(pabyPtr, &x, 8);
681 memcpy(pabyPtr+8, &y, 8);
682 if ( b3d )
683 {
684 double z = poPt->getZ();
685 memcpy(pabyPtrZ, &z, 8);
686 }
687
688 /* Swap if necessary */
689 if( OGR_SWAP( wkbNDR ) )
690 {
691 CPL_SWAPDOUBLE( pabyPtr );
692 CPL_SWAPDOUBLE( pabyPtr + 8 );
693 if( b3d )
694 CPL_SWAPDOUBLE( pabyPtrZ );
695 }
696
697 /* Advance the write pointers */
698 pabyPtr += 16;
699 if ( b3d )
700 pabyPtrZ += 8;
701 }
702
703 return OGRERR_NONE;
704 }
705
706 /* -------------------------------------------------------------------- */
707 /* MULTILINESTRING and MULTILINESTRINGZ */
708 /* -------------------------------------------------------------------- */
709 else if ( nOGRType == wkbMultiLineString )
710 {
711 OGRMultiLineString *poMLine = (OGRMultiLineString*)poGeom;
712
713 /* Write in the part count */
714 GUInt32 nPartsLsb = CPL_LSBWORD32( nParts );
715 memcpy( pabyPtr, &nPartsLsb, 4 );
716 pabyPtr += 4;
717
718 /* Write in the total point count */
719 GUInt32 nPointsLsb = CPL_LSBWORD32( nPoints );
720 memcpy( pabyPtr, &nPointsLsb, 4 );
721 pabyPtr += 4;
722
723 /* Just past the partindex[nparts] array */
724 unsigned char* pabyPoints = pabyPtr + 4*nParts;
725
726 int nPointIndexCount = 0;
727
728 for( GUInt32 i = 0; i < nParts; i++ )
729 {
730 const OGRLineString *poLine = (OGRLineString*)(poMLine->getGeometryRef(i));
731
732 /* Skip empties */
733 if ( poLine->IsEmpty() )
734 continue;
735
736 int nLineNumPoints = poLine->getNumPoints();
737
738 /* Write in the part index */
739 GUInt32 nPartIndex = CPL_LSBWORD32( nPointIndexCount );
740 memcpy( pabyPtr, &nPartIndex, 4 );
741
742 /* Write in the point data */
743 poLine->getPoints((OGRRawPoint*)pabyPoints, (double*)pabyPtrZ);
744
745 /* Swap if necessary */
746 if( OGR_SWAP( wkbNDR ) )
747 {
748 for( int k = 0; k < nLineNumPoints; k++ )
749 {
750 CPL_SWAPDOUBLE( pabyPoints + 16*k );
751 CPL_SWAPDOUBLE( pabyPoints + 16*k + 8 );
752 if( b3d )
753 CPL_SWAPDOUBLE( pabyPtrZ + 8*k );
754 }
755 }
756
757 nPointIndexCount += nLineNumPoints;
758
759 /* Advance the write pointers */
760 pabyPtr += 4;
761 pabyPoints += 16 * nLineNumPoints;
762 if ( b3d )
763 pabyPtrZ += 8 * nLineNumPoints;
764 }
765
766 return OGRERR_NONE;
767
768 }
769 /* -------------------------------------------------------------------- */
770 /* MULTIPOLYGON and MULTIPOLYGONZ */
771 /* -------------------------------------------------------------------- */
772 else if ( nOGRType == wkbMultiPolygon )
773 {
774 OGRMultiPolygon *poMPoly = (OGRMultiPolygon*)poGeom;
775
776 /* Write in the part count */
777 GUInt32 nPartsLsb = CPL_LSBWORD32( nParts );
778 memcpy( pabyPtr, &nPartsLsb, 4 );
779 pabyPtr += 4;
780
781 /* Write in the total point count */
782 GUInt32 nPointsLsb = CPL_LSBWORD32( nPoints );
783 memcpy( pabyPtr, &nPointsLsb, 4 );
784 pabyPtr += 4;
785
786 /* -------------------------------------------------------------------- */
787 /* Now we have to visit each ring and write an index number into */
788 /* the parts list, and the coordinates into the points list. */
789 /* to do it in one pass, we will use three write pointers. */
790 /* pabyPtr writes the part indexes */
791 /* pabyPoints writes the xy coordinates */
792 /* pabyPtrZ writes the z coordinates */
793 /* -------------------------------------------------------------------- */
794
795 /* Just past the partindex[nparts] array */
796 unsigned char* pabyPoints = pabyPtr + 4*nParts;
797
798 int nPointIndexCount = 0;
799
800 for( int i = 0; i < poMPoly->getNumGeometries(); i++ )
801 {
802 OGRPolygon *poPoly = (OGRPolygon*)(poMPoly->getGeometryRef(i));
803
804 /* Skip empties */
805 if ( poPoly->IsEmpty() )
806 continue;
807
808 int nRings = 1 + poPoly->getNumInteriorRings();
809
810 for( int j = 0; j < nRings; j++ )
811 {
812 /* Check our Ring and condition it */
813 OGRLinearRing *poRing;
814 if ( j == 0 )
815 {
816 poRing = poPoly->getExteriorRing();
817 /* Outer ring must be clockwise */
818 if ( ! poRing->isClockwise() )
819 poRing->reverseWindingOrder();
820 }
821 else
822 {
823 poRing = poPoly->getInteriorRing(j-1);
824 /* Inner rings should be anti-clockwise */
825 if ( poRing->isClockwise() )
826 poRing->reverseWindingOrder();
827 }
828
829 int nRingNumPoints = poRing->getNumPoints();
830
831 /* Cannot write closed rings to shape */
832 if( nRingNumPoints <= 2 || ! poRing->get_IsClosed() )
833 return OGRERR_FAILURE;
834
835 /* Write in the part index */
836 GUInt32 nPartIndex = CPL_LSBWORD32( nPointIndexCount );
837 memcpy( pabyPtr, &nPartIndex, 4 );
838
839 /* Write in the point data */
840 poRing->getPoints((OGRRawPoint*)pabyPoints, (double*)pabyPtrZ);
841
842 /* Swap if necessary */
843 if( OGR_SWAP( wkbNDR ) )
844 {
845 for( int k = 0; k < nRingNumPoints; k++ )
846 {
847 CPL_SWAPDOUBLE( pabyPoints + 16*k );
848 CPL_SWAPDOUBLE( pabyPoints + 16*k + 8 );
849 if( b3d )
850 CPL_SWAPDOUBLE( pabyPtrZ + 8*k );
851 }
852 }
853
854 nPointIndexCount += nRingNumPoints;
855 /* Advance the write pointers */
856 pabyPtr += 4;
857 pabyPoints += 16 * nRingNumPoints;
858 if ( b3d )
859 pabyPtrZ += 8 * nRingNumPoints;
860 }
861 }
862
863 return OGRERR_NONE;
864
865 }
866 else
867 {
868 return OGRERR_UNSUPPORTED_OPERATION;
869 }
870
871 }
872
873
874 /************************************************************************/
875 /* OGRWriteMultiPatchToShapeBin() */
876 /************************************************************************/
877
OGRWriteMultiPatchToShapeBin(OGRGeometry * poGeom,GByte ** ppabyShape,int * pnBytes)878 OGRErr OGRWriteMultiPatchToShapeBin( OGRGeometry *poGeom,
879 GByte **ppabyShape,
880 int *pnBytes )
881 {
882 if( wkbFlatten(poGeom->getGeometryType()) != wkbMultiPolygon )
883 return OGRERR_UNSUPPORTED_OPERATION;
884
885 poGeom->closeRings();
886 OGRMultiPolygon *poMPoly = (OGRMultiPolygon*)poGeom;
887 int nParts = 0;
888 int* panPartStart = NULL;
889 int* panPartType = NULL;
890 int nPoints = 0;
891 OGRRawPoint* poPoints = NULL;
892 double* padfZ = NULL;
893 int nBeginLastPart = 0;
894
895 for( int j = 0; j < poMPoly->getNumGeometries(); j++ )
896 {
897 OGRPolygon *poPoly = (OGRPolygon*)(poMPoly->getGeometryRef(j));
898 int nRings = poPoly->getNumInteriorRings() + 1;
899
900 /* Skip empties */
901 if ( poPoly->IsEmpty() )
902 continue;
903
904 OGRLinearRing *poRing = poPoly->getExteriorRing();
905 if( nRings == 1 && poRing->getNumPoints() == 4 )
906 {
907 if( nParts > 0 &&
908 ((panPartType[nParts-1] == SHPP_TRIANGLES && nPoints - panPartStart[nParts-1] == 3) ||
909 panPartType[nParts-1] == SHPP_TRIFAN) &&
910 poRing->getX(0) == poPoints[nBeginLastPart].x &&
911 poRing->getY(0) == poPoints[nBeginLastPart].y &&
912 poRing->getZ(0) == padfZ[nBeginLastPart] &&
913 poRing->getX(1) == poPoints[nPoints-1].x &&
914 poRing->getY(1) == poPoints[nPoints-1].y &&
915 poRing->getZ(1) == padfZ[nPoints-1] )
916 {
917 panPartType[nParts-1] = SHPP_TRIFAN;
918
919 poPoints = (OGRRawPoint*)CPLRealloc(poPoints,
920 (nPoints + 1) * sizeof(OGRRawPoint));
921 padfZ = (double*)CPLRealloc(padfZ, (nPoints + 1) * sizeof(double));
922 poPoints[nPoints].x = poRing->getX(2);
923 poPoints[nPoints].y = poRing->getY(2);
924 padfZ[nPoints] = poRing->getZ(2);
925 nPoints ++;
926 }
927 else if( nParts > 0 &&
928 ((panPartType[nParts-1] == SHPP_TRIANGLES && nPoints - panPartStart[nParts-1] == 3) ||
929 panPartType[nParts-1] == SHPP_TRISTRIP) &&
930 poRing->getX(0) == poPoints[nPoints-2].x &&
931 poRing->getY(0) == poPoints[nPoints-2].y &&
932 poRing->getZ(0) == padfZ[nPoints-2] &&
933 poRing->getX(1) == poPoints[nPoints-1].x &&
934 poRing->getY(1) == poPoints[nPoints-1].y &&
935 poRing->getZ(1) == padfZ[nPoints-1] )
936 {
937 panPartType[nParts-1] = SHPP_TRISTRIP;
938
939 poPoints = (OGRRawPoint*)CPLRealloc(poPoints,
940 (nPoints + 1) * sizeof(OGRRawPoint));
941 padfZ = (double*)CPLRealloc(padfZ, (nPoints + 1) * sizeof(double));
942 poPoints[nPoints].x = poRing->getX(2);
943 poPoints[nPoints].y = poRing->getY(2);
944 padfZ[nPoints] = poRing->getZ(2);
945 nPoints ++;
946 }
947 else
948 {
949 if( nParts == 0 || panPartType[nParts-1] != SHPP_TRIANGLES )
950 {
951 nBeginLastPart = nPoints;
952
953 panPartStart = (int*)CPLRealloc(panPartStart, (nParts + 1) * sizeof(int));
954 panPartType = (int*)CPLRealloc(panPartType, (nParts + 1) * sizeof(int));
955 panPartStart[nParts] = nPoints;
956 panPartType[nParts] = SHPP_TRIANGLES;
957 nParts ++;
958 }
959
960 poPoints = (OGRRawPoint*)CPLRealloc(poPoints,
961 (nPoints + 3) * sizeof(OGRRawPoint));
962 padfZ = (double*)CPLRealloc(padfZ, (nPoints + 3) * sizeof(double));
963 for(int i=0;i<3;i++)
964 {
965 poPoints[nPoints+i].x = poRing->getX(i);
966 poPoints[nPoints+i].y = poRing->getY(i);
967 padfZ[nPoints+i] = poRing->getZ(i);
968 }
969 nPoints += 3;
970 }
971 }
972 else
973 {
974 panPartStart = (int*)CPLRealloc(panPartStart, (nParts + nRings) * sizeof(int));
975 panPartType = (int*)CPLRealloc(panPartType, (nParts + nRings) * sizeof(int));
976
977 for ( int i = 0; i < nRings; i++ )
978 {
979 panPartStart[nParts + i] = nPoints;
980 if ( i == 0 )
981 {
982 poRing = poPoly->getExteriorRing();
983 panPartType[nParts + i] = SHPP_OUTERRING;
984 }
985 else
986 {
987 poRing = poPoly->getInteriorRing(i-1);
988 panPartType[nParts + i] = SHPP_INNERRING;
989 }
990 poPoints = (OGRRawPoint*)CPLRealloc(poPoints,
991 (nPoints + poRing->getNumPoints()) * sizeof(OGRRawPoint));
992 padfZ = (double*)CPLRealloc(padfZ,
993 (nPoints + poRing->getNumPoints()) * sizeof(double));
994 for( int k = 0; k < poRing->getNumPoints(); k++ )
995 {
996 poPoints[nPoints+k].x = poRing->getX(k);
997 poPoints[nPoints+k].y = poRing->getY(k);
998 padfZ[nPoints+k] = poRing->getZ(k);
999 }
1000 nPoints += poRing->getNumPoints();
1001 }
1002
1003 nParts += nRings;
1004 }
1005 }
1006
1007 int nShpSize = 4; /* All types start with integer type number */
1008 nShpSize += 16 * 2; /* xy bbox */
1009 nShpSize += 4; /* nparts */
1010 nShpSize += 4; /* npoints */
1011 nShpSize += 4 * nParts; /* panPartStart[nparts] */
1012 nShpSize += 4 * nParts; /* panPartType[nparts] */
1013 nShpSize += 8 * 2 * nPoints; /* xy points */
1014 nShpSize += 16; /* z bbox */
1015 nShpSize += 8 * nPoints; /* z points */
1016
1017 *pnBytes = nShpSize;
1018 *ppabyShape = (GByte*) CPLMalloc(nShpSize);
1019
1020 GByte* pabyPtr = *ppabyShape;
1021
1022 /* Write in the type number and advance the pointer */
1023 GUInt32 nGType = CPL_LSBWORD32( SHPT_MULTIPATCH );
1024 memcpy( pabyPtr, &nGType, 4 );
1025 pabyPtr += 4;
1026
1027 OGREnvelope3D envelope;
1028 poGeom->getEnvelope(&envelope);
1029 memcpy( pabyPtr, &(envelope.MinX), 8 );
1030 memcpy( pabyPtr+8, &(envelope.MinY), 8 );
1031 memcpy( pabyPtr+8+8, &(envelope.MaxX), 8 );
1032 memcpy( pabyPtr+8+8+8, &(envelope.MaxY), 8 );
1033
1034 int i;
1035 /* Swap box if needed. Shape doubles are always LSB */
1036 if( OGR_SWAP( wkbNDR ) )
1037 {
1038 for ( i = 0; i < 4; i++ )
1039 CPL_SWAPDOUBLE( pabyPtr + 8*i );
1040 }
1041 pabyPtr += 32;
1042
1043 /* Write in the part count */
1044 GUInt32 nPartsLsb = CPL_LSBWORD32( nParts );
1045 memcpy( pabyPtr, &nPartsLsb, 4 );
1046 pabyPtr += 4;
1047
1048 /* Write in the total point count */
1049 GUInt32 nPointsLsb = CPL_LSBWORD32( nPoints );
1050 memcpy( pabyPtr, &nPointsLsb, 4 );
1051 pabyPtr += 4;
1052
1053 for( i = 0; i < nParts; i ++ )
1054 {
1055 int nPartStart = CPL_LSBWORD32(panPartStart[i]);
1056 memcpy( pabyPtr, &nPartStart, 4 );
1057 pabyPtr += 4;
1058 }
1059 for( i = 0; i < nParts; i ++ )
1060 {
1061 int nPartType = CPL_LSBWORD32(panPartType[i]);
1062 memcpy( pabyPtr, &nPartType, 4 );
1063 pabyPtr += 4;
1064 }
1065
1066 memcpy(pabyPtr, poPoints, 2 * 8 * nPoints);
1067 /* Swap box if needed. Shape doubles are always LSB */
1068 if( OGR_SWAP( wkbNDR ) )
1069 {
1070 for ( i = 0; i < 2 * nPoints; i++ )
1071 CPL_SWAPDOUBLE( pabyPtr + 8*i );
1072 }
1073 pabyPtr += 2 * 8 * nPoints;
1074
1075 memcpy( pabyPtr, &(envelope.MinZ), 8 );
1076 memcpy( pabyPtr+8, &(envelope.MaxZ), 8 );
1077 if( OGR_SWAP( wkbNDR ) )
1078 {
1079 for ( i = 0; i < 2; i++ )
1080 CPL_SWAPDOUBLE( pabyPtr + 8*i );
1081 }
1082 pabyPtr += 16;
1083
1084 memcpy(pabyPtr, padfZ, 8 * nPoints);
1085 /* Swap box if needed. Shape doubles are always LSB */
1086 if( OGR_SWAP( wkbNDR ) )
1087 {
1088 for ( i = 0; i < nPoints; i++ )
1089 CPL_SWAPDOUBLE( pabyPtr + 8*i );
1090 }
1091 pabyPtr += 8 * nPoints;
1092
1093 CPLFree(panPartStart);
1094 CPLFree(panPartType);
1095 CPLFree(poPoints);
1096 CPLFree(padfZ);
1097
1098 return OGRERR_NONE;
1099 }
1100
1101 /************************************************************************/
1102 /* OGRCreateFromShapeBin() */
1103 /* */
1104 /* Translate shapefile binary representation to an OGR */
1105 /* geometry. */
1106 /************************************************************************/
1107
OGRCreateFromShapeBin(GByte * pabyShape,OGRGeometry ** ppoGeom,int nBytes)1108 OGRErr OGRCreateFromShapeBin( GByte *pabyShape,
1109 OGRGeometry **ppoGeom,
1110 int nBytes )
1111
1112 {
1113 *ppoGeom = NULL;
1114
1115 if( nBytes < 4 )
1116 {
1117 CPLError(CE_Failure, CPLE_AppDefined,
1118 "Shape buffer size (%d) too small",
1119 nBytes);
1120 return OGRERR_FAILURE;
1121 }
1122
1123 /* -------------------------------------------------------------------- */
1124 /* Detect zlib compressed shapes and uncompress buffer if necessary */
1125 /* NOTE: this seems to be an undocumented feature, even in the */
1126 /* extended_shapefile_format.pdf found in the FileGDB API documentation*/
1127 /* -------------------------------------------------------------------- */
1128 if( nBytes >= 14 &&
1129 pabyShape[12] == 0x78 && pabyShape[13] == 0xDA /* zlib marker */)
1130 {
1131 GInt32 nUncompressedSize, nCompressedSize;
1132 memcpy( &nUncompressedSize, pabyShape + 4, 4 );
1133 memcpy( &nCompressedSize, pabyShape + 8, 4 );
1134 CPL_LSBPTR32( &nUncompressedSize );
1135 CPL_LSBPTR32( &nCompressedSize );
1136 if (nCompressedSize + 12 == nBytes &&
1137 nUncompressedSize > 0)
1138 {
1139 GByte* pabyUncompressedBuffer = (GByte*)VSIMalloc(nUncompressedSize);
1140 if (pabyUncompressedBuffer == NULL)
1141 {
1142 CPLError(CE_Failure, CPLE_OutOfMemory,
1143 "Cannot allocate %d bytes to uncompress zlib buffer",
1144 nUncompressedSize);
1145 return OGRERR_FAILURE;
1146 }
1147
1148 size_t nRealUncompressedSize = 0;
1149 if( CPLZLibInflate( pabyShape + 12, nCompressedSize,
1150 pabyUncompressedBuffer, nUncompressedSize,
1151 &nRealUncompressedSize ) == NULL )
1152 {
1153 CPLError(CE_Failure, CPLE_AppDefined,
1154 "CPLZLibInflate() failed");
1155 VSIFree(pabyUncompressedBuffer);
1156 return OGRERR_FAILURE;
1157 }
1158
1159 OGRErr eErr = OGRCreateFromShapeBin(pabyUncompressedBuffer,
1160 ppoGeom,
1161 nRealUncompressedSize);
1162
1163 VSIFree(pabyUncompressedBuffer);
1164
1165 return eErr;
1166 }
1167 }
1168
1169 int nSHPType = pabyShape[0];
1170
1171 /* -------------------------------------------------------------------- */
1172 /* Return a NULL geometry when SHPT_NULL is encountered. */
1173 /* Watch out, null return does not mean "bad data" it means */
1174 /* "no geometry here". Watch the OGRErr for the error status */
1175 /* -------------------------------------------------------------------- */
1176 if ( nSHPType == SHPT_NULL )
1177 {
1178 *ppoGeom = NULL;
1179 return OGRERR_NONE;
1180 }
1181
1182 // CPLDebug( "PGeo",
1183 // "Shape type read from PGeo data is nSHPType = %d",
1184 // nSHPType );
1185
1186 /* -------------------------------------------------------------------- */
1187 /* TODO: These types include additional attributes including */
1188 /* non-linear segments and such. They should be handled. */
1189 /* This is documented in the extended_shapefile_format.pdf */
1190 /* from the FileGDB API */
1191 /* -------------------------------------------------------------------- */
1192 switch( nSHPType )
1193 {
1194 case SHPT_GENERALPOLYLINE:
1195 nSHPType = SHPT_ARC;
1196 break;
1197 case SHPT_GENERALPOLYGON:
1198 nSHPType = SHPT_POLYGON;
1199 break;
1200 case SHPT_GENERALPOINT:
1201 nSHPType = SHPT_POINT;
1202 break;
1203 case SHPT_GENERALMULTIPOINT:
1204 nSHPType = SHPT_MULTIPOINT;
1205 break;
1206 case SHPT_GENERALMULTIPATCH:
1207 nSHPType = SHPT_MULTIPATCH;
1208 }
1209
1210 /* ==================================================================== */
1211 /* Extract vertices for a Polygon or Arc. */
1212 /* ==================================================================== */
1213 if( nSHPType == SHPT_ARC
1214 || nSHPType == SHPT_ARCZ
1215 || nSHPType == SHPT_ARCM
1216 || nSHPType == SHPT_ARCZM
1217 || nSHPType == SHPT_POLYGON
1218 || nSHPType == SHPT_POLYGONZ
1219 || nSHPType == SHPT_POLYGONM
1220 || nSHPType == SHPT_POLYGONZM
1221 || nSHPType == SHPT_MULTIPATCH
1222 || nSHPType == SHPT_MULTIPATCHM)
1223 {
1224 GInt32 nPoints, nParts;
1225 int i, nOffset;
1226 GInt32 *panPartStart;
1227 GInt32 *panPartType = NULL;
1228
1229 if (nBytes < 44)
1230 {
1231 CPLError(CE_Failure, CPLE_AppDefined,
1232 "Corrupted Shape : nBytes=%d, nSHPType=%d", nBytes, nSHPType);
1233 return OGRERR_FAILURE;
1234 }
1235
1236 /* -------------------------------------------------------------------- */
1237 /* Extract part/point count, and build vertex and part arrays */
1238 /* to proper size. */
1239 /* -------------------------------------------------------------------- */
1240 memcpy( &nPoints, pabyShape + 40, 4 );
1241 memcpy( &nParts, pabyShape + 36, 4 );
1242
1243 CPL_LSBPTR32( &nPoints );
1244 CPL_LSBPTR32( &nParts );
1245
1246 if (nPoints < 0 || nParts < 0 ||
1247 nPoints > 50 * 1000 * 1000 || nParts > 10 * 1000 * 1000)
1248 {
1249 CPLError(CE_Failure, CPLE_AppDefined, "Corrupted Shape : nPoints=%d, nParts=%d.",
1250 nPoints, nParts);
1251 return OGRERR_FAILURE;
1252 }
1253
1254 int bHasZ = ( nSHPType == SHPT_POLYGONZ
1255 || nSHPType == SHPT_POLYGONZM
1256 || nSHPType == SHPT_ARCZ
1257 || nSHPType == SHPT_ARCZM
1258 || nSHPType == SHPT_MULTIPATCH
1259 || nSHPType == SHPT_MULTIPATCHM );
1260
1261 int bIsMultiPatch = ( nSHPType == SHPT_MULTIPATCH || nSHPType == SHPT_MULTIPATCHM );
1262
1263 /* With the previous checks on nPoints and nParts, */
1264 /* we should not overflow here and after */
1265 /* since 50 M * (16 + 8 + 8) = 1 600 MB */
1266 int nRequiredSize = 44 + 4 * nParts + 16 * nPoints;
1267 if ( bHasZ )
1268 {
1269 nRequiredSize += 16 + 8 * nPoints;
1270 }
1271 if( bIsMultiPatch )
1272 {
1273 nRequiredSize += 4 * nParts;
1274 }
1275 if (nRequiredSize > nBytes)
1276 {
1277 CPLError(CE_Failure, CPLE_AppDefined,
1278 "Corrupted Shape : nPoints=%d, nParts=%d, nBytes=%d, nSHPType=%d, nRequiredSize=%d",
1279 nPoints, nParts, nBytes, nSHPType, nRequiredSize);
1280 return OGRERR_FAILURE;
1281 }
1282
1283 panPartStart = (GInt32 *) VSICalloc(nParts,sizeof(GInt32));
1284 if (panPartStart == NULL)
1285 {
1286 CPLError(CE_Failure, CPLE_OutOfMemory,
1287 "Not enough memory for shape (nPoints=%d, nParts=%d)", nPoints, nParts);
1288 return OGRERR_FAILURE;
1289 }
1290
1291 /* -------------------------------------------------------------------- */
1292 /* Copy out the part array from the record. */
1293 /* -------------------------------------------------------------------- */
1294 memcpy( panPartStart, pabyShape + 44, 4 * nParts );
1295 for( i = 0; i < nParts; i++ )
1296 {
1297 CPL_LSBPTR32( panPartStart + i );
1298
1299 /* We check that the offset is inside the vertex array */
1300 if (panPartStart[i] < 0 ||
1301 panPartStart[i] >= nPoints)
1302 {
1303 CPLError(CE_Failure, CPLE_AppDefined,
1304 "Corrupted Shape : panPartStart[%d] = %d, nPoints = %d",
1305 i, panPartStart[i], nPoints);
1306 CPLFree(panPartStart);
1307 return OGRERR_FAILURE;
1308 }
1309 if (i > 0 && panPartStart[i] <= panPartStart[i-1])
1310 {
1311 CPLError(CE_Failure, CPLE_AppDefined,
1312 "Corrupted Shape : panPartStart[%d] = %d, panPartStart[%d] = %d",
1313 i, panPartStart[i], i - 1, panPartStart[i - 1]);
1314 CPLFree(panPartStart);
1315 return OGRERR_FAILURE;
1316 }
1317 }
1318
1319 nOffset = 44 + 4*nParts;
1320
1321 /* -------------------------------------------------------------------- */
1322 /* If this is a multipatch, we will also have parts types. */
1323 /* -------------------------------------------------------------------- */
1324 if( bIsMultiPatch )
1325 {
1326 panPartType = (GInt32 *) VSICalloc(nParts,sizeof(GInt32));
1327 if (panPartType == NULL)
1328 {
1329 CPLError(CE_Failure, CPLE_OutOfMemory,
1330 "Not enough memory for panPartType for shape (nPoints=%d, nParts=%d)", nPoints, nParts);
1331 CPLFree(panPartStart);
1332 return OGRERR_FAILURE;
1333 }
1334
1335 memcpy( panPartType, pabyShape + nOffset, 4*nParts );
1336 for( i = 0; i < nParts; i++ )
1337 {
1338 CPL_LSBPTR32( panPartType + i );
1339 }
1340 nOffset += 4*nParts;
1341 }
1342
1343 /* -------------------------------------------------------------------- */
1344 /* Copy out the vertices from the record. */
1345 /* -------------------------------------------------------------------- */
1346 double *padfX = (double *) VSIMalloc(sizeof(double)*nPoints);
1347 double *padfY = (double *) VSIMalloc(sizeof(double)*nPoints);
1348 double *padfZ = (double *) VSICalloc(sizeof(double),nPoints);
1349 if (padfX == NULL || padfY == NULL || padfZ == NULL)
1350 {
1351 CPLFree( panPartStart );
1352 CPLFree( panPartType );
1353 CPLFree( padfX );
1354 CPLFree( padfY );
1355 CPLFree( padfZ );
1356 CPLError(CE_Failure, CPLE_OutOfMemory,
1357 "Not enough memory for shape (nPoints=%d, nParts=%d)", nPoints, nParts);
1358 return OGRERR_FAILURE;
1359 }
1360
1361 for( i = 0; i < nPoints; i++ )
1362 {
1363 memcpy(padfX + i, pabyShape + nOffset + i * 16, 8 );
1364 memcpy(padfY + i, pabyShape + nOffset + i * 16 + 8, 8 );
1365 CPL_LSBPTR64( padfX + i );
1366 CPL_LSBPTR64( padfY + i );
1367 }
1368
1369 nOffset += 16*nPoints;
1370
1371 /* -------------------------------------------------------------------- */
1372 /* If we have a Z coordinate, collect that now. */
1373 /* -------------------------------------------------------------------- */
1374 if( bHasZ )
1375 {
1376 for( i = 0; i < nPoints; i++ )
1377 {
1378 memcpy( padfZ + i, pabyShape + nOffset + 16 + i*8, 8 );
1379 CPL_LSBPTR64( padfZ + i );
1380 }
1381
1382 nOffset += 16 + 8*nPoints;
1383 }
1384
1385 /* -------------------------------------------------------------------- */
1386 /* Build corresponding OGR objects. */
1387 /* -------------------------------------------------------------------- */
1388 if( nSHPType == SHPT_ARC
1389 || nSHPType == SHPT_ARCZ
1390 || nSHPType == SHPT_ARCM
1391 || nSHPType == SHPT_ARCZM )
1392 {
1393 /* -------------------------------------------------------------------- */
1394 /* Arc - As LineString */
1395 /* -------------------------------------------------------------------- */
1396 if( nParts == 1 )
1397 {
1398 OGRLineString *poLine = new OGRLineString();
1399 *ppoGeom = poLine;
1400
1401 poLine->setPoints( nPoints, padfX, padfY, padfZ );
1402 }
1403
1404 /* -------------------------------------------------------------------- */
1405 /* Arc - As MultiLineString */
1406 /* -------------------------------------------------------------------- */
1407 else
1408 {
1409 OGRMultiLineString *poMulti = new OGRMultiLineString;
1410 *ppoGeom = poMulti;
1411
1412 for( i = 0; i < nParts; i++ )
1413 {
1414 OGRLineString *poLine = new OGRLineString;
1415 int nVerticesInThisPart;
1416
1417 if( i == nParts-1 )
1418 nVerticesInThisPart = nPoints - panPartStart[i];
1419 else
1420 nVerticesInThisPart =
1421 panPartStart[i+1] - panPartStart[i];
1422
1423 poLine->setPoints( nVerticesInThisPart,
1424 padfX + panPartStart[i],
1425 padfY + panPartStart[i],
1426 padfZ + panPartStart[i] );
1427
1428 poMulti->addGeometryDirectly( poLine );
1429 }
1430 }
1431 } /* ARC */
1432
1433 /* -------------------------------------------------------------------- */
1434 /* Polygon */
1435 /* -------------------------------------------------------------------- */
1436 else if( nSHPType == SHPT_POLYGON
1437 || nSHPType == SHPT_POLYGONZ
1438 || nSHPType == SHPT_POLYGONM
1439 || nSHPType == SHPT_POLYGONZM )
1440 {
1441 if (nParts != 0)
1442 {
1443 if (nParts == 1)
1444 {
1445 OGRPolygon *poOGRPoly = new OGRPolygon;
1446 *ppoGeom = poOGRPoly;
1447 OGRLinearRing *poRing = new OGRLinearRing;
1448 int nVerticesInThisPart = nPoints - panPartStart[0];
1449
1450 poRing->setPoints( nVerticesInThisPart,
1451 padfX + panPartStart[0],
1452 padfY + panPartStart[0],
1453 padfZ + panPartStart[0] );
1454
1455 poOGRPoly->addRingDirectly( poRing );
1456 }
1457 else
1458 {
1459 OGRGeometry *poOGR = NULL;
1460 OGRPolygon** tabPolygons = new OGRPolygon*[nParts];
1461
1462 for( i = 0; i < nParts; i++ )
1463 {
1464 tabPolygons[i] = new OGRPolygon();
1465 OGRLinearRing *poRing = new OGRLinearRing;
1466 int nVerticesInThisPart;
1467
1468 if( i == nParts-1 )
1469 nVerticesInThisPart = nPoints - panPartStart[i];
1470 else
1471 nVerticesInThisPart =
1472 panPartStart[i+1] - panPartStart[i];
1473
1474 poRing->setPoints( nVerticesInThisPart,
1475 padfX + panPartStart[i],
1476 padfY + panPartStart[i],
1477 padfZ + panPartStart[i] );
1478 tabPolygons[i]->addRingDirectly(poRing);
1479 }
1480
1481 int isValidGeometry;
1482 const char* papszOptions[] = { "METHOD=ONLY_CCW", NULL };
1483 poOGR = OGRGeometryFactory::organizePolygons(
1484 (OGRGeometry**)tabPolygons, nParts, &isValidGeometry, papszOptions );
1485
1486 if (!isValidGeometry)
1487 {
1488 CPLError(CE_Warning, CPLE_AppDefined,
1489 "Geometry of polygon cannot be translated to Simple Geometry. "
1490 "All polygons will be contained in a multipolygon.\n");
1491 }
1492
1493 *ppoGeom = poOGR;
1494 delete[] tabPolygons;
1495 }
1496 }
1497 } /* polygon */
1498
1499 /* -------------------------------------------------------------------- */
1500 /* Multipatch */
1501 /* -------------------------------------------------------------------- */
1502 else if( bIsMultiPatch )
1503 {
1504 *ppoGeom = OGRCreateFromMultiPatch( nParts,
1505 panPartStart,
1506 panPartType,
1507 nPoints,
1508 padfX,
1509 padfY,
1510 padfZ );
1511 }
1512
1513 CPLFree( panPartStart );
1514 CPLFree( panPartType );
1515 CPLFree( padfX );
1516 CPLFree( padfY );
1517 CPLFree( padfZ );
1518
1519 if (*ppoGeom != NULL)
1520 (*ppoGeom)->setCoordinateDimension( bHasZ ? 3 : 2 );
1521
1522 return OGRERR_NONE;
1523 }
1524
1525 /* ==================================================================== */
1526 /* Extract vertices for a MultiPoint. */
1527 /* ==================================================================== */
1528 else if( nSHPType == SHPT_MULTIPOINT
1529 || nSHPType == SHPT_MULTIPOINTM
1530 || nSHPType == SHPT_MULTIPOINTZ
1531 || nSHPType == SHPT_MULTIPOINTZM )
1532 {
1533 GInt32 nPoints;
1534 GInt32 nOffsetZ;
1535 int i;
1536
1537 int bHasZ = ( nSHPType == SHPT_MULTIPOINTZ
1538 || nSHPType == SHPT_MULTIPOINTZM );
1539
1540
1541 memcpy( &nPoints, pabyShape + 36, 4 );
1542 CPL_LSBPTR32( &nPoints );
1543
1544 if (nPoints < 0 || nPoints > 50 * 1000 * 1000 )
1545 {
1546 CPLError(CE_Failure, CPLE_AppDefined, "Corrupted Shape : nPoints=%d.",
1547 nPoints);
1548 return OGRERR_FAILURE;
1549 }
1550
1551 nOffsetZ = 40 + 2*8*nPoints + 2*8;
1552
1553 OGRMultiPoint *poMultiPt = new OGRMultiPoint;
1554 *ppoGeom = poMultiPt;
1555
1556 for( i = 0; i < nPoints; i++ )
1557 {
1558 double x, y, z;
1559 OGRPoint *poPt = new OGRPoint;
1560
1561 /* Copy X */
1562 memcpy(&x, pabyShape + 40 + i*16, 8);
1563 CPL_LSBPTR64(&x);
1564 poPt->setX(x);
1565
1566 /* Copy Y */
1567 memcpy(&y, pabyShape + 40 + i*16 + 8, 8);
1568 CPL_LSBPTR64(&y);
1569 poPt->setY(y);
1570
1571 /* Copy Z */
1572 if ( bHasZ )
1573 {
1574 memcpy(&z, pabyShape + nOffsetZ + i*8, 8);
1575 CPL_LSBPTR64(&z);
1576 poPt->setZ(z);
1577 }
1578
1579 poMultiPt->addGeometryDirectly( poPt );
1580 }
1581
1582 poMultiPt->setCoordinateDimension( bHasZ ? 3 : 2 );
1583
1584 return OGRERR_NONE;
1585 }
1586
1587 /* ==================================================================== */
1588 /* Extract vertices for a point. */
1589 /* ==================================================================== */
1590 else if( nSHPType == SHPT_POINT
1591 || nSHPType == SHPT_POINTM
1592 || nSHPType == SHPT_POINTZ
1593 || nSHPType == SHPT_POINTZM )
1594 {
1595 /* int nOffset; */
1596 double dfX, dfY, dfZ = 0;
1597
1598 int bHasZ = (nSHPType == SHPT_POINTZ || nSHPType == SHPT_POINTZM);
1599
1600 if (nBytes < 4 + 8 + 8 + ((bHasZ) ? 8 : 0))
1601 {
1602 CPLError(CE_Failure, CPLE_AppDefined,
1603 "Corrupted Shape : nBytes=%d, nSHPType=%d", nBytes, nSHPType);
1604 return OGRERR_FAILURE;
1605 }
1606
1607 memcpy( &dfX, pabyShape + 4, 8 );
1608 memcpy( &dfY, pabyShape + 4 + 8, 8 );
1609
1610 CPL_LSBPTR64( &dfX );
1611 CPL_LSBPTR64( &dfY );
1612 /* nOffset = 20 + 8; */
1613
1614 if( bHasZ )
1615 {
1616 memcpy( &dfZ, pabyShape + 4 + 16, 8 );
1617 CPL_LSBPTR64( &dfZ );
1618 }
1619
1620 *ppoGeom = new OGRPoint( dfX, dfY, dfZ );
1621 (*ppoGeom)->setCoordinateDimension( bHasZ ? 3 : 2 );
1622
1623 return OGRERR_NONE;
1624 }
1625
1626 CPLError(CE_Failure, CPLE_AppDefined,
1627 "Unsupported geometry type: %d",
1628 nSHPType );
1629
1630 return OGRERR_FAILURE;
1631 }
1632