source: subversion/applications/utils/coastcheck/closeshp-ram.c @ 29269

Last change on this file since 29269 was 9804, checked in by jonb, 11 years ago

Experimental version which checks for overlap in returned results and caches shapefiles in ram. This branch is temporary until the ideas get merged back to the latest closeshp.c version

File size: 41.9 KB
Line 
1/* closeshp - Takes two shape files and bounding box and returns a new shape
2 * file with closed polygon as determined by the bounding box.
3 * By Martijn van Oosterhout <kleptog@svana.org> Copyright 2008
4 * Licence: GPL
5 *
6 * This is a similar idea to close-areas.pl for osmarender, except
7 * considerably easier due to the fact that all the segments are in order
8 * and connected already. There are two shapefiles: one where coast2shp has
9 * managed to close the way, these are type POLYGON. The second is where the
10 * closing didn't work, these are type ARC.
11 *
12 * The main differences between the two are:
13 * - If a closed polygon completely surrounds (but does not touch the edge
14 *   of) the bounding box it produce something, where as an ARC does not.
15 * - If a polygon intersects the box you always get something whereas if a
16 *   arc ends in the box, it gets dropped.
17 */
18
19#define _GNU_SOURCE
20#include <stdio.h>
21#include <unistd.h>
22#include <stdlib.h>
23#include <string.h>
24#include <assert.h>
25#include <getopt.h>
26#include <errno.h>
27#include <math.h>
28
29#include <shapefil.h>
30
31#define VERBOSE 0
32#define TEST 0
33#define UNUSED __attribute__((unused))
34
35#define MERC_MAX (20037508.34f)
36#define DIVISIONS 400
37#define MERC_BLOCK (2*MERC_MAX/DIVISIONS)
38/* Number of mercator metres the tiles overlap, so the antialising doesn't cause wierd effects */
39#define TILE_OVERLAP  10000
40#define MAX_NODES 500000
41#define MAX_SEGS  100
42int MAX_SUBAREAS;   /* Was a define, not anymore. Auto grown array, starting at... */
43#define INIT_MAX_SUBAREAS 1024
44#define MAX_NODES_PER_ARC 10000
45
46static double *v_x, *v_y, *v_z;
47static int *Parts, *Parttypes;
48
49static char *small_poly_list;
50static SHPHandle shp_out;
51static DBFHandle dbf_out;
52
53/* Column of error flag */
54#define DBF_OUT_ERROR   0
55#define DBF_OUT_TILE_X   1
56#define DBF_OUT_TILE_Y   2
57struct segment
58{
59  SHPHandle shp;   /* Shapefile of this segment */
60  int index;       /* Index of object in shapefile */
61  double sx, sy;   /* Where the segment enters the box */
62  double ex, ey;   /* Where the segment leaves the box */
63  int snode, enode; /* Numbers of nodes in polygon belonging to this segment */
64  double sa, ea;   /* Angle (in radians) of start/endpoint from centre of box */
65  int next;        /* In joining stage, have we used this one already */
66};
67
68/* This is used to track areas that fall completely within a box. After we
69 * have dealt will all the areas that intersect, we go through these to
70 * determine if they are contained. If a negative subarea is contained with
71 * a positive one, we append it to the shape. If a positive has no
72 * surrounding we just output it. And a positive within a positive is an
73 * error. */
74struct subarea
75{
76  SHPHandle shp;   /* Shapefile of this segment */
77  int index;       /* Index of this object in shapefile */
78  double x, y;     /* representative point on shape */
79  double areasize; /* Positive area is land, negative is lake */
80  int used;        /* Have we output this subarea yet */
81};
82
83struct state
84{
85  int x,y;
86 
87  double lb[2];
88  double rt[2];
89
90//  SHPHandle shp;
91//  DBFHandle dbf;
92//  SHPTree  *shx;
93   
94  struct segment seg_list[MAX_SEGS];
95  int seg_count;
96 
97  struct subarea *sub_areas;
98  int subarea_count;
99  int subarea_nodecount;
100 
101  int enclosed;
102};
103
104void OutputSegs( struct state *state );
105void Process( struct state *state, SHPHandle shp, DBFHandle dbf, SHPTree *shx, int polygon , SHPObject **objects);
106static double CalcArea( const SHPObject *obj );
107static void SplitCoastlines( SHPHandle shp_poly, DBFHandle dbf_poly, SHPHandle shp_arc, DBFHandle dbf_arc, char *out_filename );
108static int contains( double x, double y, double *v_x, double *v_y, int vertices );
109static int CopyShapeToArray( SHPHandle shp, int index, int snode, int enode, int node_count );
110static void ResizeSubareas( struct state *state, int count );
111void ProcessSubareas(struct state *state, int *pc, int *nc);
112
113int main( int argc, char *argv[] )
114{
115  char out_filename[256];
116 
117  if( argc != 4 )
118  {
119    fprintf( stderr, "closeshp poly_shape arc_shape output_shape\n" );
120    return 1;
121  }
122  char *poly_file = argv[1];
123  char *arc_file = argv[2];
124  char *out_file = argv[3];
125
126  int poly_count, arc_count;
127
128  /* open shapefiles and dbf files */
129  SHPHandle shp_poly = SHPOpen( poly_file, "rb" );
130  if( !shp_poly )
131  {
132    fprintf( stderr, "Couldn't open '%s': %s\n", poly_file, strerror(errno));
133    return 1;
134  }
135
136  SHPHandle shp_arc = SHPOpen( arc_file, "rb" );
137  if( !shp_arc )
138  {
139    fprintf( stderr, "Couldn't open '%s': %s\n", arc_file, strerror(errno));
140    return 1;
141  }
142 
143  DBFHandle dbf_poly = DBFOpen( poly_file, "rb" );
144  if( !dbf_poly )
145  {
146    fprintf( stderr, "Couldn't open DBF file for '%s'\n", poly_file );
147    return 1;
148  }
149  if( DBFGetFieldIndex( dbf_poly, "way_id" ) != 2 )
150  {
151    fprintf( stderr, "Unexpected format DBF file '%s'\n", poly_file );
152    return 1;
153  }
154
155  DBFHandle dbf_arc = DBFOpen( arc_file, "rb" );
156  if( !dbf_arc )
157  {
158    fprintf( stderr, "Couldn't open DBF file for '%s'\n", arc_file );
159    return 1;
160  }
161  if( DBFGetFieldIndex( dbf_arc, "way_id" ) != 2 )
162  {
163    fprintf( stderr, "Unexpected format DBF file '%s'\n", arc_file );
164    return 1;
165  }
166
167  /* Split coastlines into arc no longer than MAX_NODES_PER_ARC long */
168  sprintf( out_filename, "%s_i", out_file );
169  SplitCoastlines( shp_poly, dbf_poly, shp_arc, dbf_arc, out_filename );
170 
171  sprintf( out_filename, "%s_p", out_file );
172  shp_out = SHPCreate( out_filename, SHPT_POLYGON );
173  if( !shp_out )
174  {
175    fprintf( stderr, "Couldn't create shapefile '%s': %s\n", out_file, strerror(errno));
176    return 1;
177  }
178 
179  dbf_out = DBFCreate( out_filename );
180  if( !dbf_out )
181  {
182    fprintf( stderr, "Couldn't create DBF '%s': %s\n", out_file, strerror(errno));
183    return 1;
184  }
185//  DBFAddField( dbf_out, "way_id", FTInteger, 11, 0 );
186//  DBFAddField( dbf_out, "orientation", FTInteger, 1, 0 );
187  DBFAddField( dbf_out, "error", FTInteger, 1, 0 );
188  DBFAddField( dbf_out, "tile_x", FTInteger, 4, 0 );
189  DBFAddField( dbf_out, "tile_y", FTInteger, 4, 0 );
190
191  /* Check shapefiles are the right type */
192  {
193    int type;
194    SHPGetInfo( shp_poly, &poly_count, &type, NULL, NULL );
195    if( type != SHPT_POLYGON )
196    {
197      fprintf( stderr, "'%s' is not a POLYGON shapefile\n", poly_file );
198      return 1;
199    }
200    SHPGetInfo( shp_arc, &arc_count, &type, NULL, NULL );
201    if( type != SHPT_ARC )
202    {
203      fprintf( stderr, "'%s' is not a ARC shapefile\n", arc_file );
204      return 1;
205    }
206  }
207
208  // Cache poly & arc shapefiles
209  int i;
210  SHPObject **poly_objs = malloc(sizeof(SHPObject *) * poly_count);
211  if (!poly_objs) {
212      fprintf(stderr, "malloc for %d poly objects failed\n", poly_count);
213      exit(1);
214  }
215   for (i=0; i<poly_count; i++)
216      poly_objs[i] = SHPReadObject(shp_poly, i);
217
218  SHPObject **arc_objs = malloc(sizeof(SHPObject *) * arc_count);
219  if (!arc_objs) {
220      fprintf(stderr, "malloc for %d arc objects failed\n", arc_count);
221      exit(1);
222  }
223  for (i=0; i<arc_count; i++)
224      arc_objs[i] = SHPReadObject(shp_arc, i);
225
226  /* Build indexes on files, we need them... */
227  SHPTree *shx_poly = SHPCreateTree( shp_poly, 2, 0, NULL, NULL );
228  SHPTree *shx_arc  = SHPCreateTree( shp_arc, 2, 0, NULL, NULL );
229  if( !shx_poly || !shx_arc )
230  {
231    fprintf( stderr, "Couldn't open shape indexes\n" );
232    return 1;
233  }
234 
235  /* The poly tells us which polygons are big enough to need splitting */
236  small_poly_list = malloc( (poly_count>>3) + 1 );
237  if( !small_poly_list )
238  {
239    fprintf( stderr, "Couldn't allocate poly list\n" );
240    return 1;
241  }
242  memset( small_poly_list, 0, (poly_count>>3) + 1 );
243 
244  v_x = malloc( MAX_NODES * sizeof(double) );
245  v_y = malloc( MAX_NODES * sizeof(double) );
246  v_z = malloc( MAX_NODES * sizeof(double) );
247 
248  if( !v_x || !v_y || !v_z)
249  {
250    fprintf( stderr, "Couldn't allocate memory for nodes\n" );
251    return 1;
252  }
253  memset( v_z, 0, MAX_NODES * sizeof(double) );
254  struct state state;
255  memset( &state, 0, sizeof(state) );
256  ResizeSubareas(&state, INIT_MAX_SUBAREAS);
257
258#if !TEST 
259  for( int i=0; i<DIVISIONS; i++ )
260    for( int j=0; j<DIVISIONS; j++ )  //Divide the world into mercator blocks approx 100km x 100km
261#else
262  for( int i=307; i<=307; i++ )
263    for( int j=203; j<=203; j++ )  //Divide the world into mercator blocks approx 100km x 100km
264#endif
265    {
266      state.x = i;
267      state.y = j;
268     
269      double left   = -MERC_MAX + (i*MERC_BLOCK) - TILE_OVERLAP;
270      double right  = -MERC_MAX + ((i+1)*MERC_BLOCK) + TILE_OVERLAP;
271      double bottom = -MERC_MAX + (j*MERC_BLOCK) - TILE_OVERLAP;
272      double top    = -MERC_MAX + ((j+1)*MERC_BLOCK) + TILE_OVERLAP;
273     
274      if( left  < -MERC_MAX ) left  = -MERC_MAX;
275      if( right > +MERC_MAX ) right = +MERC_MAX;
276     
277      state.lb[0] = left;
278      state.lb[1] = bottom;
279      state.rt[0] = right;
280      state.rt[1] = top;
281     
282      if(isatty(STDERR_FILENO))
283//        fprintf( stderr, "\rProcessing (%d,%d)  (%.2f,%.2f)-(%.2f,%.2f)   ", i, j, left, bottom, right, top );
284        fprintf( stderr, "\rProcessing (%d,%d)  (%.2f,%.2f)-(%.2f,%.2f)   ", i, j, state.lb[0], state.lb[1], state.rt[0], state.rt[1] );
285      state.seg_count = 0;
286      state.subarea_count = 0;
287      state.subarea_nodecount = 0;
288      state.enclosed = 0;
289     
290      Process( &state, shp_poly, dbf_poly, shx_poly, 1 , poly_objs);
291      Process( &state, shp_arc,  dbf_arc,  shx_arc,  0 , arc_objs);
292     
293      OutputSegs( &state );
294    }
295   
296  SHPDestroyTree( shx_poly );
297  SHPDestroyTree( shx_arc );
298  DBFClose( dbf_poly );
299  DBFClose( dbf_arc );
300  DBFClose( dbf_out );
301  SHPClose( shp_poly );
302  SHPClose( shp_arc );
303  SHPClose( shp_out );
304 
305  printf("\n");
306  return 0;
307}
308
309static void ResizeSubareas( struct state *state, int count )
310{
311  if( count < MAX_SUBAREAS )
312  {
313    fprintf( stderr, "Tried to resize smaller??? (%d < %d)\n", count, MAX_SUBAREAS );
314    exit(1);
315  }
316  fprintf( stderr, "Resizing subarea array to %d\n", count );
317 
318  struct subarea *new_sa = malloc( count * sizeof(struct subarea) );
319  free(Parts);
320  Parts = malloc( 2 * count * sizeof(int) );  // parts and Parttypes are allocated in one chunk
321 
322  if( !new_sa || !Parts )
323  {
324    fprintf( stderr, "Out of memory resizing subarea array (count=%d)\n", count );
325    exit(1);
326  }
327 
328  MAX_SUBAREAS = count;
329  memcpy( new_sa, state->sub_areas, state->subarea_count * sizeof(struct subarea) );
330  free(state->sub_areas);
331  state->sub_areas = new_sa;
332 
333  Parttypes = Parts + MAX_SUBAREAS;
334}
335
336static void SplitCoastlines2( int show, SHPHandle shp, DBFHandle dbf, SHPHandle shp_arc_out, DBFHandle dbf_arc_out )
337{
338  int count;
339  SHPGetInfo( shp, &count, NULL, NULL, NULL );
340  for( int i=0; i<count; i++ )
341  {
342    SHPObject *obj = SHPReadObject( shp, i );
343    int way_id = DBFReadIntegerAttribute( dbf, i, 2 );
344   
345    if( obj->nVertices <= MAX_NODES_PER_ARC )
346    {
347      int new_id = SHPWriteObject( shp_arc_out, -1, obj );
348      if( new_id < 0 ) { fprintf( stderr, "Output failure: %m\n"); exit(1); }
349      DBFWriteIntegerAttribute( dbf_arc_out, new_id, 0, way_id );
350      DBFWriteIntegerAttribute( dbf_arc_out, new_id, 1, show );
351      DBFWriteIntegerAttribute( dbf_arc_out, new_id, 2, obj->nVertices < 4 );  /* Flag not real objects */
352      SHPDestroyObject(obj);
353      continue;
354    }
355    int arcs = (obj->nVertices / MAX_NODES_PER_ARC) + 1;
356    int len = (obj->nVertices / arcs) + 1;
357//    printf( "Splitting object with %d vertices, len=%d, arcs=%d\n", obj->nVertices, len, arcs );
358   
359    for( int j=0; j<arcs; j++ )
360    {
361      int this_len = (j==arcs-1)? obj->nVertices - (j*len): len+1;
362//      printf( "Subobject start=%d, length=%d\n", j*len, this_len );
363      SHPObject *new_obj = SHPCreateSimpleObject( SHPT_ARC, this_len, &obj->padfX[j*len], &obj->padfY[j*len], &obj->padfZ[j*len] );
364      int new_id = SHPWriteObject( shp_arc_out, -1, new_obj );
365      if( new_id < 0 ) { fprintf( stderr, "Output failure: %m\n"); exit(1); }
366      DBFWriteIntegerAttribute( dbf_arc_out, new_id, 0, way_id );
367      DBFWriteIntegerAttribute( dbf_arc_out, new_id, 1, show );
368      DBFWriteIntegerAttribute( dbf_arc_out, new_id, 2, 0 );
369      SHPDestroyObject(new_obj);
370    }
371    SHPDestroyObject(obj);
372  }
373}
374
375/* The first two params are currently unused, but if people ever want to get
376 * access to the completed bits of coastline, this is where to change it */
377static void SplitCoastlines( SHPHandle shp_poly UNUSED, DBFHandle dbf_poly UNUSED, SHPHandle shp_arc, DBFHandle dbf_arc, char *out_filename )
378{
379  SHPHandle shp_arc_out = SHPCreate( out_filename, SHPT_ARC );
380  if( !shp_arc_out )
381  {
382    fprintf( stderr, "Couldn't create shapefile '%s': %s\n", out_filename, strerror(errno));
383    return;
384  }
385 
386  DBFHandle dbf_arc_out = DBFCreate( out_filename );
387  if( !dbf_arc_out )
388  {
389    fprintf( stderr, "Couldn't create DBF '%s': %s\n", out_filename, strerror(errno));
390    return;
391  }
392  DBFAddField( dbf_arc_out, "way_id",   FTInteger, 11, 0 );
393  DBFAddField( dbf_arc_out, "complete", FTInteger, 11, 0 );
394  DBFAddField( dbf_arc_out, "error",    FTInteger, 11, 0 );
395
396//  SplitCoastlines2( shp_poly, dbf_poly, shp_arc_out, dbf_arc_out );
397  SplitCoastlines2( 0, shp_arc,  dbf_arc,  shp_arc_out, dbf_arc_out );
398 
399  SHPClose( shp_arc_out );
400  DBFClose( dbf_arc_out );
401}
402
403static double CalcArea( const SHPObject *obj )
404{
405  int i;
406  double base_x = obj->dfXMin;
407  double base_y = obj->dfYMin;
408  double area = 0;
409  int n = (obj->nParts <= 1) ? obj->nVertices : obj->panPartStart[1];
410//  if(VERBOSE) printf( "CalcArea: n=%d\n", n );
411  for( i=0; i<n; i++ )
412  {
413    int p = ((i==0)?n:i)-1;
414    double x1 = obj->padfX[p] - base_x;
415    double x2 = obj->padfX[i] - base_x;
416    double y1 = obj->padfY[p] - base_y;
417    double y2 = obj->padfY[i] - base_y;
418   
419//    if(VERBOSE) printf( "i=%d (%f,%f)-(%f,%f) %f\n", i, x1,y1,x2,y2,x1*y2 - y1*x2);
420    area += x1*y2 - y1*x2;
421  }
422  return area;
423}
424
425static const int table[3][3] = { {6, 5, 4}, {7, -1, 3}, { 0, 1, 2 } };
426
427/* Determines the quadrant the given point is in. -1 means it's inside the box */
428/* The rule is: left and bottom edges are in, right and top edges are out */
429static inline int GetPosition( double *lb, double *rt, double X, double Y )
430{
431  int x, y;
432
433  x = (X >= rt[0]) - ( X < lb[0] ) + 1;
434  y = (Y >= rt[1]) - ( Y < lb[1] ) + 1;
435
436  return table[x][y];
437}
438
439struct intersect
440{
441  double x, y;  /* Coordinates */
442  double t;     /* 0..1, where the intersection happens */
443};
444
445/* The rule is: left and bottom edges are in, right and top edges are out */
446/* All corners are out, except the lower left one */
447static int CalculateIntersections( double x1, double y1, double x2, double y2, 
448                                   double *lb, double *rt, struct intersect *intersections )
449{
450  int count = 0;
451  double x, y, t;
452 
453  /* Left side */
454  if( (x1 < lb[0]) != (x2 < lb[0]) )
455  {
456    /* Determine intersection */
457    x = lb[0];
458    t = (x1-lb[0]) / (x1-x2);
459    y = y1 + t*(y2-y1);
460   
461    if( lb[1] <= y && y < rt[1] ) /* Include only if in range */
462      intersections[count++] = (struct intersect){ x: x, y: y, t: t };
463  }
464  /* Right side */
465  if( (x1 >= rt[0]) != (x2 >= rt[0]) )
466  {
467    /* Determine intersection */
468    x = rt[0];
469    t = (x1-rt[0]) / (x1-x2);
470    y = y1 + t*(y2-y1);
471   
472    if( lb[1] <= y && y < rt[1] ) /* Include only if in range */
473      intersections[count++] = (struct intersect){ x: x, y: y, t: t };
474  }
475  /* Top side */
476  if( (y1 >= rt[1]) != (y2 >= rt[1]) )
477  {
478    /* Determine intersection */
479    y = rt[1];
480    t = (y1-rt[1]) / (y1-y2);
481    x = x1 + t*(x2-x1);
482   
483    if( lb[0] <= x && x < rt[0] ) /* Include only if in range */
484      intersections[count++] = (struct intersect){ x: x, y: y, t: t };
485  }
486  /* Bottom side */
487  if( (y1 < lb[1]) != (y2 < lb[1]) )
488  {
489    /* Determine intersection */
490    y = lb[1];
491    t = (y1-lb[1]) / (y1-y2);
492    x = x1 + t*(x2-x1);
493   
494    if( lb[0] <= x && x < rt[0] ) /* Include only if in range */
495      intersections[count++] = (struct intersect){ x: x, y: y, t: t };
496  }
497
498  /* Check the count, if we went over we killed the caller's stack... */
499  if( count > 2 )
500  {
501    fprintf( stderr, "Too many intersections (%d)\n", count );
502    exit(1);
503  }
504 
505  if( count == 2 )
506  {
507    /* If there's two intersections, reorder them to match the intersection order */
508    if( intersections[0].t > intersections[1].t )
509    {
510      struct intersect i;
511      i = intersections[1];
512      intersections[1] = intersections[0];
513      intersections[0] = i;
514    }
515  }
516  if(count == 0)
517  {
518    printf( "\nCalculate intersections: (%.2f,%.2f)-(%.2f,%.2f) hit %d\n", x1, y1, x2, y2, count );
519    for( int i=0; i<count; i++ )
520      printf( "   (%.2f,%.2f) t=%.6f\n", intersections[i].x, intersections[i].y, intersections[i].t );
521  } 
522  return count;
523}
524
525static int seg_compare( const void *a, const void *b )
526{
527  const struct segment *aa = (struct segment*)a;
528  const struct segment *bb = (struct segment*)b;
529 
530  if( aa->sa < bb->sa )
531    return -1;
532  if( aa->sa > bb->sa )
533    return 1;
534  return 0;
535}
536
537/* We currently don't use anything from the source DBF file, but the cabability is there */
538void Process( struct state *state, SHPHandle shp, DBFHandle dbf UNUSED, SHPTree *shx, int polygon, SHPObject **objects)
539{
540  int count;
541  int *list = SHPTreeFindLikelyShapes( shx, state->lb, state->rt, &count );
542  int poly_start;
543
544  for( int poly = 0; poly < count; poly++ )
545  {
546    /* Here we track parts that have gone across the box */
547    int vertex;
548    int intersected = 0;
549   
550    int id = list[poly];
551    /* If this is a small polygon, it's been printed, so skip it */
552    if( polygon && small_poly_list[ id>>3 ] & (1<<(id&7)) )
553      continue;
554
555    /* Now we have a candidate object, we need to process it */
556//    SHPObject *obj = SHPReadObject( shp, id );
557    SHPObject *obj = objects[id];
558   if (!obj)
559      continue;
560
561    /* If it's got less than 4 vertices it's not a real object */
562    /* No need to mark it as error here, done in SplitCoastlines */
563    if( obj->nVertices < 4 )
564    {
565    #if 0
566      if( polygon )
567      {
568        int new_id = SHPWriteObject( shp_out, -1, obj );
569        if( new_id < 0 ) { fprintf( stderr, "Output failure: %m\n"); exit(1); }
570        DBFWriteIntegerAttribute( dbf_out, new_id, DBF_OUT_ERROR, 1 );
571        DBFWriteIntegerAttribute( dbf_out, new_id, DBF_OUT_TILE_X, state->x );
572        DBFWriteIntegerAttribute( dbf_out, new_id, DBF_OUT_TILE_Y, state->y );
573      }
574    #endif
575      //SHPDestroyObject( obj );
576      continue;
577    }
578
579    // check for being outside (search results include false positives)
580        if (!SHPCheckBoundsOverlap( state->lb, state->rt, &(obj->dfXMin), &(obj->dfXMax), 2))
581        {
582            //printf( "Shape %d: not in area of interest, but fetched.\n",  id);
583            //SHPDestroyObject( obj );
584            continue;
585        }
586
587
588    if( polygon &&
589        state->lb[0] < obj->dfXMin && obj->dfXMax < state->rt[0] &&
590        state->lb[1] < obj->dfYMin && obj->dfYMax < state->rt[1] )
591    {
592      int sa = state->subarea_count;
593      if( sa > MAX_SUBAREAS-5 )
594        ResizeSubareas( state, 2*MAX_SUBAREAS );
595
596      state->sub_areas[sa].shp = shp;
597      state->sub_areas[sa].index = id;
598      state->sub_areas[sa].x = obj->padfX[0];
599      state->sub_areas[sa].y = obj->padfY[0];
600      state->sub_areas[sa].areasize = CalcArea( obj );
601      state->sub_areas[sa].used = 0;
602
603      state->subarea_count++;
604      state->subarea_nodecount += obj->nVertices;
605      //SHPDestroyObject( obj );
606      continue;
607    }
608 
609    if(VERBOSE) fprintf( stderr, "\nProcessing object %d (%d vertices)\n", poly, obj->nVertices );
610    /* First we need to move along the object until we leave the box. For polygons we
611     * know it will eventually, since we determined already this object is larger
612     * than the box */
613
614    for( vertex=0; vertex < obj->nVertices && GetPosition( state->lb, state->rt, obj->padfX[vertex], obj->padfY[vertex] ) == -1; vertex++ )
615      ;
616
617    if( vertex == obj->nVertices )
618    {
619      if( polygon )
620        fprintf( stderr, "Object %d did not leave box (%d vertices, polygon:%d) (%.2f,%.2f-%.2f,%.2f)\n", id, 
621                obj->nVertices, polygon, obj->dfXMin, obj->dfYMin, obj->dfXMax, obj->dfYMax );
622      //SHPDestroyObject( obj );
623      continue;
624    }
625    /* We need to mark this point, so when we loop back we know where to stop */
626    poly_start = vertex+1;
627   
628    /* This tracks the current sector. */
629    int curr_sect = GetPosition( state->lb, state->rt, obj->padfX[vertex], obj->padfY[vertex] );
630   
631    int winding = 0, max_winding = 0, min_winding = 0;
632   
633    /* We need this flag so when we loop round a polygon we know if we're at the first node or the last */
634    int started = 0;
635
636    /* The basic trick is to analyse each new point and determine the
637     * sector. As long as the sector doesn't change, we're cool. They're
638     * numbered in such a way that any line has to go through the sectors in
639     * numerical order (modulo 8) *or* it intersects the box. The
640     * reason we do this is that if the line never intersects the box, we need
641     * to determine the "winding number" determine if we're inside or outside. */
642
643    if( poly_start == obj->nVertices )
644      poly_start = 1;
645//    fprintf( stderr, "poly_start=%d, vertex=%d\n", poly_start, vertex );
646    for(;;)
647    {
648      vertex++;
649      /* First we need to handle the step to the next node */
650      if( polygon )
651      {
652        /* For polygons we loop around to the start point */
653        if( vertex == obj->nVertices )
654        {
655//          fprintf( stderr, "\nLooping..." );
656//          sleep(1);
657          vertex = 1;
658        }
659        if( vertex == poly_start && started )
660          break;
661      }
662      else
663      {
664        /* Else we just stop at the end */
665        if( vertex == obj->nVertices )
666          break;
667      }
668      started = 1;
669     
670      if( vertex >= obj->nVertices )  /* Shouldn't happen */
671      {
672        fprintf( stderr, "Somehow %d >= %d\n", vertex, obj->nVertices );
673        break;
674      }
675     
676      int sect = GetPosition( state->lb, state->rt, obj->padfX[vertex], obj->padfY[vertex] );
677      if( sect == curr_sect )
678        continue;
679      if(VERBOSE)
680      printf("Moved from %d to %d\n", curr_sect, sect );
681      /* Now we know we've moved to another sector, so we need to know the intersection points */
682      struct intersect intersections[2];
683      int int_count;
684
685//      if( id == 133912 ) fprintf(stderr, "\nVertex %d: Moved from sector %d to %d", vertex, curr_sect, sect );
686     
687      /* If we moved to adjacent positive sector, we don't need to check for intersections... */
688      if( sect >= 0 && curr_sect >= 0 && ( ((curr_sect-sect)&7) == 1 || ((sect-curr_sect)&7) == 1 ) )
689        int_count = 0;
690      else
691        int_count = CalculateIntersections( obj->padfX[vertex-1], obj->padfY[vertex-1], 
692                                            obj->padfX[vertex],   obj->padfY[vertex],
693                                            state->lb, state->rt,
694                                            intersections );
695
696      /* There are corner cases with the calculations of intersections, if
697       * you move exactly on to the edge. With floating point numbers you'd
698       * think it was possibly, but, well, it is. Espescially around the 0
699       * meridian. In this case we get zero intersections even though we
700       * changed inside/outside. What we do is basically pretend we havn't
701       * changed sector at all and check the next point */
702      if( int_count == 0 && (sect == -1 || curr_sect == -1) )
703      {
704        fprintf( stderr, "Went from %d to %d without intersection.\n"
705                      "line (%f,%f)-(%f,%f), box (%f,%f)-(%f,%f)\n",
706                      curr_sect, sect,
707                      obj->padfX[vertex-1], obj->padfY[vertex-1],
708                      obj->padfX[vertex],   obj->padfY[vertex],
709                      state->lb[0], state->lb[1], state->rt[0], state->rt[1] );
710#if TEST
711        continue;
712#endif
713      }
714      /* Another possibility is that we went from positive to positive with
715       * only one intersection. Not good. Have not yet thought of a
716       * heuristic to deal with this case */
717      if( int_count == 1 && ( (sect != -1) == (curr_sect != -1) ) )
718      {
719        fprintf( stderr, "Went from %d to %d with 1 intersection.\n"
720                      "line (%f,%f)-(%f,%f), box (%f,%f)-(%f,%f)\n",
721                      curr_sect, sect,
722                      obj->padfX[vertex-1], obj->padfY[vertex-1],
723                      obj->padfX[vertex],   obj->padfY[vertex],
724                      state->lb[0], state->lb[1], state->rt[0], state->rt[1] );
725      }
726     
727      /* finally, if we got two intersections, we mave have passed straight through the box */
728      if( int_count == 2 && (sect == -1 || curr_sect == -1) )
729      {
730        fprintf( stderr, "Went from %d to %d with 2 intersections.\n"
731                      "line (%f,%f)-(%f,%f), box (%f,%f)-(%f,%f)\n",
732                      curr_sect, sect,
733                      obj->padfX[vertex-1], obj->padfY[vertex-1],
734                      obj->padfX[vertex],   obj->padfY[vertex],
735                      state->lb[0], state->lb[1], state->rt[0], state->rt[1] );
736      }
737     
738      struct segment *seg_ptr = &state->seg_list[state->seg_count];
739     
740      /* Now we know the number of intersections. */
741      if( int_count == 0 )
742      {
743        /* If we have no intersections, we were outside and we're still
744         * outside. Then we need to track the winding number. With no
745         * intersections we can move a maximum of three sections, so we can
746         * tell if we moved clockwise or anti-clockwise. */
747        if( !intersected )
748        {
749          int diff = (sect - curr_sect) & 7;  /* Like mod 8, but always positive */
750          if( diff > 4 )
751            diff -= 8;
752           
753          winding += diff;
754         
755          if( winding < min_winding )
756            min_winding = winding;
757          if( winding > max_winding )
758            max_winding = winding;
759        }
760      }
761      else if( int_count == 1 )
762      {
763        /* If we have one intersection, we went from inside to out, or vice
764         * versa. So we need to add a section or finish an old one. */
765        intersected = 1;
766        if( curr_sect != -1 )
767        {
768          /* Going in... */
769          seg_ptr->shp = shp;
770          seg_ptr->index = id;
771          seg_ptr->sx = intersections[0].x;
772          seg_ptr->sy = intersections[0].y;
773          seg_ptr->snode = vertex;
774        }
775        else
776        {
777          /* Going out... */
778          seg_ptr->ex = intersections[0].x;
779          seg_ptr->ey = intersections[0].y;
780          seg_ptr->enode = vertex-1;
781         
782          (state->seg_count)++;
783        }
784      }
785      else if( int_count == 2 )
786      {
787        /* If we have two intersections, we went straight through. So we need
788         * to make a segment with no internal nodes */
789        intersected = 1;
790       
791        /* Going in and out in one go... */
792        seg_ptr->shp = NULL;
793        seg_ptr->index = -2;
794        seg_ptr->sx = intersections[0].x;
795        seg_ptr->sy = intersections[0].y;
796        seg_ptr->snode = -1;
797        seg_ptr->ex = intersections[1].x;
798        seg_ptr->ey = intersections[1].y;
799        seg_ptr->enode = -1;
800       
801        (state->seg_count)++;
802      }
803      curr_sect = sect;
804    }
805    /* Generally, if we have a high winding number we consider ourselves
806     * inside. However, in the case of bits that stick out, while the ends
807     * may show to not have a high winding number, the min/max winding
808     * numbers can show that the shape made a loop around the box, even
809     * though the endpoints may not look like it. So we also trigger if
810     * the end winding number >3 and the greatest difference at least 6 */
811   
812    /* However, it's more complicated than that. We could be enclosed by a
813     * continental landmass, but also by a large lake. So we must take the
814     * *smallest* enclosing area to determine whether we're enclosed or not.
815     * To facilitate this the enclosed fields tracks the area of the
816     * smallest encloser and the sign determines te positive or negativeness
817     * of the area. */
818   
819    if( !intersected )
820    {
821      double enclosed = 0;
822      if ((winding >= 3) || (winding >= 2 && (max_winding-min_winding) > 6))
823      {
824        if(VERBOSE)
825          printf( "Decided enclosed: intersected=%d, winding=%d, max_winding=%d, min_winding=%d\n", intersected, winding, max_winding, min_winding );
826        enclosed = +1;
827      }
828      if (winding < -4) 
829      {
830        if(VERBOSE)
831          printf( "Decided unenclosed: intersected=%d, winding=%d, max_winding=%d, min_winding=%d\n", intersected, winding, max_winding, min_winding );
832        enclosed = -1;
833      }
834
835      if( enclosed )
836      {
837        // Here we don't need to worry too much about small areas, since
838        // they're unlikely to enclose anything
839        double size = CalcArea(obj);
840        int intsize = copysign( ceil(10*log2(1+fabs(size))), enclosed ); // Scale the size down to a number that will fit in an integer
841        if( state->enclosed == 0 )
842          state->enclosed = intsize;
843        else if( abs(state->enclosed) > abs(intsize) )
844          state->enclosed = intsize;
845         
846        if(VERBOSE)
847          printf( "(%d,%d) New state->enclosed: %d, size=%f\n", state->x, state->y, state->enclosed, size );
848      }
849    }
850   
851    //SHPDestroyObject(obj);
852  }
853  free(list);
854//  printf( "segcount: %d\n", state->seg_count );
855}
856
857/* Compare function to sort subareas in *descending* order */
858static int subarea_compare( const void *a, const void *b )
859{
860  struct subarea *sa = (struct subarea *)a;
861  struct subarea *sb = (struct subarea *)b;
862 
863  double diff = fabs(sb->areasize) - fabs(sa->areasize);
864 
865  if( diff > 0 ) return +1;
866  if( diff < 0 ) return -1;
867  return 0;
868}
869
870void OutputSegs( struct state *state )
871{
872  /* At this point we've processed the whole object and have a list of
873   * segments which cross our box. We now compute the angles they make to
874   * the centre, sort them on that basis and try to make closed areas. There
875   * are as many starts as finishes so it will always finish, but a senseble
876   * end result depends on the original polygon having been sane (i.e.
877   * simple) */
878   
879  double centre_x = (state->lb[0] + state->rt[0])/2;
880  double centre_y = (state->lb[1] + state->rt[1])/2;
881
882  struct segment *seg_list = state->seg_list;
883  int seg_count = state->seg_count;
884 
885//  fprintf( stderr, "%d sub %d  ", state->subarea_count, state->subarea_nodecount );
886  // First we must sort the subareas by decreasing size
887  qsort( state->sub_areas, state->subarea_count, sizeof( state->sub_areas[0] ), subarea_compare );
888 
889  if( seg_count == 0 )
890  {
891    /* No intersections at all, so we check the winding number. With the
892     * water-on-the-right rule we're looking for a positive winding. */
893//      if( abs(winding) > 4 )
894//        printf( "\nNot intersected, winding = %d, min_winding = %d, max_winding = %d\n", winding, min_winding, max_winding );
895     
896    if( state->enclosed > 0 )  // Enclosed by a negative area does not count
897    {
898      seg_list[0].sx = seg_list[0].ex = state->lb[0];
899      seg_list[0].sy = seg_list[0].ey = centre_y;
900      seg_list[0].snode = seg_list[0].enode = -1;
901      seg_list[0].sa = +M_PI;
902      seg_list[0].ea = -M_PI;
903      seg_list[0].next = -1;
904      seg_count++;
905    }
906  }
907  else
908  {
909    /* Work out the angles */
910    for( int k=0; k<seg_count; k++ )
911    {
912      seg_list[k].sa = atan2( seg_list[k].sy - centre_y, seg_list[k].sx - centre_x );
913      seg_list[k].ea = atan2( seg_list[k].ey - centre_y, seg_list[k].ex - centre_x );
914      seg_list[k].next = -1;
915    }
916  }
917  if( seg_count > 0 )
918  {
919    /* Create the helper nodes for the corners */
920    for( int k=0; k<4; k++, seg_count++ )
921    {
922      seg_list[seg_count].sx = seg_list[seg_count].ex = (k<2) ? state->lb[0] : state->rt[0];
923      seg_list[seg_count].sy = seg_list[seg_count].ey = (((k+1)&3)<2) ? state->lb[1] : state->rt[1];
924      seg_list[seg_count].sa = seg_list[seg_count].ea = atan2( seg_list[seg_count].sy - centre_y, seg_list[seg_count].sx - centre_x );
925      seg_list[seg_count].snode = seg_list[seg_count].enode = -2;
926      seg_list[seg_count].next = -2;
927    }
928   
929    /* Sort the nodes by increasing angle */
930    qsort( seg_list, seg_count, sizeof(seg_list[0]), seg_compare );
931   
932    for(;;)
933    {
934      int part_count = 1;
935      Parts[0] = 0;
936      Parttypes[0] = SHPP_RING;
937     
938      /* First we need to find an unused segment */
939      int curr;
940      for( curr=0; curr < seg_count && seg_list[curr].next != -1; curr++ )
941        ;
942      if( curr == seg_count )
943        break;
944      int node_count = 0;
945      for(;;)
946      {
947        if(VERBOSE) printf( "Part %d: ndc=%d, curr=%d (%d-%d)\n", part_count, node_count, curr, seg_list[curr].snode, seg_list[curr].enode );
948        v_x[node_count] = seg_list[curr].sx;
949        v_y[node_count] = seg_list[curr].sy;
950        node_count++;
951
952        if( seg_list[curr].snode >= 0 )
953        {
954          node_count = CopyShapeToArray( seg_list[curr].shp, seg_list[curr].index, seg_list[curr].snode, seg_list[curr].enode, node_count );
955          v_x[node_count] = seg_list[curr].ex;
956          v_y[node_count] = seg_list[curr].ey;
957          node_count++;
958        }
959       
960        double angle = seg_list[curr].ea;
961        /* Determine the first unused segment with a start angle greater than the current angle, with wrapping */
962        int next;
963        for( next = 0; next < seg_count && seg_list[next].sa <= angle; next++ )
964          ;
965        if( next == seg_count )
966          next = 0;
967         
968        seg_list[curr].next = next;
969        /* If we come to an already used segment, we're done */
970        if( seg_list[next].next >= 0 )
971          break;
972         
973        curr = next;
974      }
975      v_x[node_count] = v_x[0];
976      v_y[node_count] = v_y[0];
977      node_count++;
978     
979      ProcessSubareas( state, &part_count, &node_count );
980     
981      if( part_count > MAX_SUBAREAS - 2 )
982        fprintf( stderr, "(%d,%d) Subarea overflow: %d > %d\n", state->x, state->y, part_count, MAX_SUBAREAS-2 );
983      if( node_count > MAX_NODES - 100 )
984        fprintf( stderr, "(%d,%d) Node overflow: %d > %d\n", state->x, state->y, node_count, MAX_NODES - 100 );
985//        fprintf( stderr, "Created object: %d verticies\n", node_count );
986//      SHPObject *shape = SHPCreateSimpleObject( SHPT_POLYGON, node_count, v_x, v_y, v_z );
987      SHPObject *shape = SHPCreateObject( SHPT_POLYGON, -1, 
988                                            part_count, Parts, Parttypes, 
989                                            node_count, v_x, v_y, v_z, NULL );
990      // If a wrongly oriented shape crosses a boundary, sometimes we can see that...
991      // Must do this prior to rewinding object
992      int inverted = (CalcArea( shape ) < 0);
993      if(VERBOSE) printf( "Created shape orientation: %f,%s\n", CalcArea( shape ), !inverted?"good":"bad" );
994      SHPRewindObject( NULL, shape );
995      int new_id = SHPWriteObject( shp_out, -1, shape );
996      if( new_id < 0 ) { fprintf( stderr, "Output failure: %m\n"); exit(1); }
997      SHPDestroyObject( shape );
998     
999      DBFWriteIntegerAttribute( dbf_out, new_id, DBF_OUT_ERROR, inverted ? 1 : 0 );
1000      DBFWriteIntegerAttribute( dbf_out, new_id, DBF_OUT_TILE_X, state->x );
1001      DBFWriteIntegerAttribute( dbf_out, new_id, DBF_OUT_TILE_Y, state->y );
1002    }
1003    /* Check for any remaining sub areas... we just output them, but mark negative ones as errors */
1004    for( int k=0; k<state->subarea_count; k++ )
1005    {
1006      struct subarea *sub_area = &state->sub_areas[k];
1007      if( sub_area->used )
1008        continue;
1009       
1010      if(VERBOSE) printf( "Remaining subarea %d (area=%f)\n", k, sub_area->areasize );
1011      int node_count = CopyShapeToArray( sub_area->shp, sub_area->index, 0, -1, 0 );
1012      int part_count = 1;
1013      Parttypes[0] = SHPP_RING;
1014     
1015      sub_area->used = 1;  // Need to mark as used first, or it will match itself
1016
1017      ProcessSubareas( state, &part_count, &node_count );
1018     
1019//      SHPObject *obj = SHPReadObject( sub_area->shp, sub_area->index );
1020      SHPObject *obj = SHPCreateObject( SHPT_POLYGON, -1, 
1021                                        part_count, Parts, Parttypes, 
1022                                        node_count, v_x, v_y, v_z, NULL );
1023      int new_id = SHPWriteObject( shp_out, -1, obj );
1024      if( new_id < 0 ) { fprintf( stderr, "Output failure: %m\n"); exit(1); }
1025      SHPDestroyObject( obj );
1026      DBFWriteIntegerAttribute( dbf_out, new_id, DBF_OUT_ERROR, (sub_area->areasize > 0)?0:1 );
1027      DBFWriteIntegerAttribute( dbf_out, new_id, DBF_OUT_TILE_X, state->x );
1028      DBFWriteIntegerAttribute( dbf_out, new_id, DBF_OUT_TILE_Y, state->y );
1029    }
1030//    printf( "\nMatching %s with %d vertices, id=%d\n", polygon ? "polygon" : "arc", obj->nVertices, id );
1031    if(VERBOSE)
1032    for( int k=0; k<seg_count; k++ )
1033    {
1034      printf( "%2d %6d (%11.2f,%11.2f)-(%11.2f,%11.2f) (%5d-%5d) (%4.2f-%4.2f) => %2d\n", 
1035          k, seg_list[k].index, seg_list[k].sx, seg_list[k].sy, seg_list[k].ex, seg_list[k].ey, 
1036          seg_list[k].snode, seg_list[k].enode, seg_list[k].sa, seg_list[k].ea,
1037          seg_list[k].next );
1038    }
1039  }
1040}
1041
1042static int contains( double x, double y, double *v_x, double *v_y, int vertices )
1043{
1044#if 1
1045  /* Referenced from Wikipedia */
1046  /* http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html */
1047    {
1048      int i, j, c = 0;
1049      for (i = 0, j = vertices-1; i < vertices; j = i++) {
1050        if ((((v_y[i]<=y) && (y<v_y[j])) ||
1051             ((v_y[j]<=y) && (y<v_y[i]))) &&
1052            (x < (v_x[j] - v_x[i]) * (y - v_y[i]) / (v_y[j] - v_y[i]) + v_x[i]))
1053
1054          c = !c;
1055      }
1056      return c;
1057    }
1058#else
1059  /* Home cooked algorithm, with a bug */
1060  double last_angle = atan2( v_y[0] - y, v_x[0] - x );
1061  double angle = 0;
1062 
1063#if TEST
1064  int debug = (round(x) == round(6538707.962709)) && (round(y) == round(2365505.004066));
1065#endif
1066
1067  double last_x = v_x[0];
1068  double last_y = v_y[0];
1069 
1070  for( int i=0; i<vertices; i++ )
1071  {
1072    /* Only bother checking if we've changed quadrants */
1073    if( i != vertices-1 &&
1074        (v_x[i] < x) == (last_x < x) &&
1075        (v_y[i] < y) == (last_y < y) )
1076      continue;
1077    last_x = v_x[i];
1078    last_y = v_y[i];
1079   
1080    double new_angle = atan2( v_y[i] - y, v_x[i] - x );
1081   
1082#if TEST
1083    if(debug)
1084      printf( "i=%d (%f,%f)-(%f,%f) last_angle=%.2f new_angle=%.2f angle=%.2f offset=%.2f\n", i, v_x[i-1], v_y[i-1], last_x, last_y, last_angle, new_angle, angle, new_angle-angle );
1085#endif
1086    /* What we want to do is set angle to the same angle as new_angle module
1087     * 2pi in such a way that it is within pi of the current angle. */
1088    double diff = fmod(fmod(new_angle-last_angle,2*M_PI) + 3*M_PI,2*M_PI) - M_PI;
1089    last_angle = new_angle;
1090    angle += diff;
1091   
1092  }
1093  double windings = round( angle / (2*M_PI) );
1094
1095#if TEST
1096  if(debug)
1097    printf( "Determined for (%.2f,%.2f) windings=%f\n", x, y, windings );
1098#endif
1099  return windings > 0;
1100#endif
1101}
1102
1103static int CopyShapeToArray( SHPHandle shp, int index, int snode, int enode, int node_count )
1104{
1105  SHPObject *obj = SHPReadObject( shp, index );
1106  int k;
1107 
1108  if( enode == -1 )   /* This mean to copy the whole object */
1109    enode = obj->nVertices - 1;
1110  //          printf( "Copying %d - %d (max %d)\n", seg_list[curr].snode, seg_list[curr].enode, obj->nVertices );
1111 
1112  /* The slightly odd construction of the loop is because it need to handle sections of polygon which may wrap around the end of the shape */
1113  if( snode != enode && enode == 0 )
1114    enode = obj->nVertices-1;
1115  // nodes[0] == node[nVertices-1]
1116  // nodes[1] == node[nVertices]
1117  for( k = snode; k != enode; k++ )
1118  {
1119    if( node_count > MAX_NODES - 100 )
1120    {
1121//      fprintf( stderr, "Node overflow...\n" );
1122      break;
1123    }
1124    if( k == obj->nVertices )
1125    {
1126      k = 1;
1127      if( k == enode )
1128        break;
1129    }
1130     
1131    v_x[node_count] = obj->padfX[k];
1132    v_y[node_count] = obj->padfY[k];
1133    node_count++;
1134  }
1135  /* k = enode now... */
1136  v_x[node_count] = obj->padfX[k];
1137  v_y[node_count] = obj->padfY[k];
1138  node_count++;
1139  SHPDestroyObject(obj);
1140 
1141  return node_count;
1142}
1143
1144void ProcessSubareas(struct state *state, int *pc, int *nc)
1145{
1146  int part_count = *pc;
1147  int node_count = *nc;
1148 
1149  if(VERBOSE) printf( "ProcessSubareas(pc=%d,nc=%d)\n", *pc, *nc );
1150 
1151  // Setting this simplifies the later code as it can be assumed that
1152  // every ring ends at the beginning of the next part
1153  Parts[part_count] = node_count;
1154  // Now we have the outer ring we need to find if any of the subareas are contained in it
1155  for( int k=0; k<state->subarea_count; k++ )
1156  {
1157    struct subarea *sub_area = &state->sub_areas[k];
1158    if( sub_area->used )  // Already used, skip
1159      continue;
1160    if( !contains( sub_area->x, sub_area->y, v_x, v_y, Parts[1] ) ) // If not contained skip
1161      continue;
1162    // Now we have to verify that this new subarea is not contained in any of the existing ones
1163    if(VERBOSE) printf( "subarea %d is contained (size=%f,%s)\n", k, sub_area->areasize, (sub_area->areasize < 0)?"good":"bad");
1164    int contained = 1;
1165    for( int m=1; m < part_count; m++ )
1166    {
1167      if( contains( sub_area->x, sub_area->y, v_x+Parts[m], v_y+Parts[m], Parts[m+1]-Parts[m] ) )
1168      {
1169        contained = 0;
1170        if(VERBOSE) printf( "subarea %d is excluded\n", k);
1171        break;
1172      }
1173    }
1174    if( !contained )
1175      continue;
1176     
1177    // We know the surrounding object has positive area
1178    if( sub_area->areasize < 0 )
1179    {
1180//          printf( "Appending subshape\n" );
1181      Parttypes[0] = SHPP_OUTERRING; /* Multipart object now */
1182      // Append to object
1183      if( part_count >= MAX_SUBAREAS-2 )
1184      {
1185       // ResizeSubareas( state, MAX_SUBAREAS*2 )
1186        fprintf( stderr, "(%d,%d) Parts array overflow (%d)\n", state->x, state->y, MAX_SUBAREAS );
1187      }
1188      else
1189      {
1190        Parttypes[part_count] = SHPP_INNERRING;
1191        Parts[part_count] = node_count;
1192        node_count = CopyShapeToArray( sub_area->shp, sub_area->index, 0, -1, node_count );
1193        Parts[part_count+1] = node_count;
1194      }
1195     
1196      part_count++;
1197    }
1198    else
1199    {
1200      // Error, copy object
1201      SHPObject *obj = SHPReadObject( sub_area->shp, sub_area->index );
1202      if(VERBOSE) printf( "Outputting error shape (x,y)=(%f,%f), n=%d\n", obj->padfX[0], obj->padfY[0], obj->nVertices );
1203      int new_id = SHPWriteObject( shp_out, -1, obj );
1204      if( new_id < 0 ) { fprintf( stderr, "Output failure: %m\n"); exit(1); }
1205      SHPDestroyObject( obj );
1206     
1207      DBFWriteIntegerAttribute( dbf_out, new_id, DBF_OUT_ERROR, 1 );
1208      DBFWriteIntegerAttribute( dbf_out, new_id, DBF_OUT_TILE_X, state->x );
1209      DBFWriteIntegerAttribute( dbf_out, new_id, DBF_OUT_TILE_Y, state->y );
1210    }
1211     
1212    sub_area->used = 1;
1213  }
1214 
1215  *nc = node_count;
1216  *pc = part_count;
1217}
Note: See TracBrowser for help on using the repository browser.