source: subversion/applications/utils/coastcheck/closeshp.c @ 9744

Last change on this file since 9744 was 7969, checked in by martinvoosterhout, 12 years ago

No new features, just plenty of error checking, growing the nodes array and
plugging some memory leaks.

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