source: subversion/applications/utils/coastcheck/coast2shp.c @ 9743

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

Make coast2shp limit lattitudes to +-85.0511 degrees to avoid problems with points at infinity

File size: 22.1 KB
Line 
1/*
2#-----------------------------------------------------------------------------
3# coast2shp - converts list of polygons into shapefilesfile into PostgreSQL
4# compatible output suitable to be rendered by mapnik
5# Use: coast2shp polygonlist coastline.osm.gz
6#-----------------------------------------------------------------------------
7# Based upon osm2pgsql:
8#   Original Python implementation by Artem Pavlenko
9#   Re-implementation by Jon Burgess, Copyright 2006
10# Reused for coast2shp by Martijn van Oosterhout 2007-2008
11#
12# This program is free software; you can redistribute it and/or
13# modify it under the terms of the GNU General Public License
14# as published by the Free Software Foundation; either version 2
15# of the License, or (at your option) any later version.
16#
17# This program is distributed in the hope that it will be useful,
18# but WITHOUT ANY WARRANTY; without even the implied warranty of
19# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
20# GNU General Public License for more details.
21#
22# You should have received a copy of the GNU General Public License
23# along with this program; if not, write to the Free Software
24# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
25#-----------------------------------------------------------------------------
26*/
27
28#define _GNU_SOURCE
29#include <stdio.h>
30#include <unistd.h>
31#include <stdlib.h>
32#include <string.h>
33#include <assert.h>
34#include <getopt.h>
35#include <errno.h>
36#include <math.h>
37
38#include <libxml/xmlstring.h>
39#include <libxml/xmlreader.h>
40
41#include <shapefil.h>
42#include <proj_api.h>
43
44#include "osmtypes.h"
45#include "keyvals.h"
46#include "input.h"
47#include "rb.h"
48
49/* Mercator projection limits for a square map.
50 * Points at 90 degrees are off at infinity
51 */
52const double merc_lat_min = -85.0511;
53const double merc_lat_max = +85.0511;
54
55int MAX_VERTICES = 1024*1024;
56
57static int count_node,    max_node;
58static int count_way,     max_way;
59static int count_rel,     max_rel;
60
61/* Since {node,way} elements are not nested we can guarantee the
62   values in an end tag must match those of the corresponding
63   start tag and can therefore be cached.
64*/
65static double node_lon, node_lat;
66static struct keyval nds;
67static int osm_id;
68struct rb_table *nodes_table;
69struct rb_table *ways_table;
70
71int verbose;
72int latlong;
73
74static projPJ pj_ll, pj_merc;
75
76static void project_init(void)
77{
78        pj_ll   = pj_init_plus("+proj=latlong +ellps=GRS80 +no_defs +over");
79        pj_merc = pj_init_plus("+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs +over");
80
81        if (!pj_ll || !pj_merc) {
82                fprintf(stderr, "Projection code failed to initialise\n");
83                exit(1);
84        }
85}
86
87static void project_exit(void)
88{
89        pj_free(pj_ll);
90        pj_ll = NULL;
91        pj_free(pj_merc);
92        pj_merc = NULL;
93}
94
95static void printStatus(void)
96{
97    if( isatty(STDERR_FILENO) )
98      fprintf(stderr, "\rProcessing: Node(%dk) Way(%dk) Relation(%dk)",
99              count_node/1000, count_way/1000, count_rel/1000);
100}
101
102
103void StartElement(xmlTextReaderPtr reader, const xmlChar *name)
104{
105    xmlChar *xid, *xlat, *xlon /**xk, *xv, *xrole, *xtype*/;
106//    char *k;
107
108    if (xmlStrEqual(name, BAD_CAST "node")) {
109        xid  = xmlTextReaderGetAttribute(reader, BAD_CAST "id");
110        xlon = xmlTextReaderGetAttribute(reader, BAD_CAST "lon");
111        xlat = xmlTextReaderGetAttribute(reader, BAD_CAST "lat");
112        assert(xid); assert(xlon); assert(xlat);
113
114        osm_id  = strtol((char *)xid, NULL, 10);
115        node_lon = strtod((char *)xlon, NULL);
116        node_lat = strtod((char *)xlat, NULL);
117
118        if (osm_id > max_node)
119            max_node = osm_id;
120
121        count_node++;
122        if (count_node%10000 == 0)
123            printStatus();
124
125        xmlFree(xid);
126        xmlFree(xlon);
127        xmlFree(xlat);
128    } else if (xmlStrEqual(name, BAD_CAST "tag")) {
129#if 0
130        xk = xmlTextReaderGetAttribute(reader, BAD_CAST "k");
131        assert(xk);
132
133        /* 'created_by' and 'source' are common and not interesting to mapnik renderer */
134        if (strcmp((char *)xk, "created_by") && strcmp((char *)xk, "source")) {
135            char *p;
136            xv = xmlTextReaderGetAttribute(reader, BAD_CAST "v");
137            assert(xv);
138            k  = (char *)xmlStrdup(xk);
139            while ((p = strchr(k, ' ')))
140                *p = '_';
141
142            addItem(&tags, k, (char *)xv, 0);
143            xmlFree(k);
144            xmlFree(xv);
145        }
146        xmlFree(xk);
147#endif
148    } else if (xmlStrEqual(name, BAD_CAST "way")) {
149        xid  = xmlTextReaderGetAttribute(reader, BAD_CAST "id");
150        assert(xid);
151        osm_id   = strtol((char *)xid, NULL, 10);
152
153        if (osm_id > max_way)
154            max_way = osm_id;
155
156        count_way++;
157        if (count_way%1000 == 0)
158            printStatus();
159
160        xmlFree(xid);
161    } else if (xmlStrEqual(name, BAD_CAST "nd")) {
162        xid  = xmlTextReaderGetAttribute(reader, BAD_CAST "ref");
163        assert(xid);
164
165        addItem(&nds, "id", (char *)xid, 0);
166
167        xmlFree(xid);
168    } else if (xmlStrEqual(name, BAD_CAST "relation")) {
169#if 0
170        xid  = xmlTextReaderGetAttribute(reader, BAD_CAST "id");
171        assert(xid);
172        osm_id   = strtol((char *)xid, NULL, 10);
173
174        if (osm_id > max_rel)
175            max_rel = osm_id;
176
177        xmlFree(xid);
178#endif
179        count_rel++;
180        if (count_rel%1000 == 0)
181            printStatus();
182
183    } else if (xmlStrEqual(name, BAD_CAST "member")) {
184#if 0
185        xrole = xmlTextReaderGetAttribute(reader, BAD_CAST "role");
186        assert(xrole);
187
188        xtype = xmlTextReaderGetAttribute(reader, BAD_CAST "type");
189        assert(xtype);
190
191        xid  = xmlTextReaderGetAttribute(reader, BAD_CAST "ref");
192        assert(xid);
193
194        /* Currently we are only interested in 'way' members since these form polygons with holes */
195        if (xmlStrEqual(xtype, BAD_CAST "way"))
196            addItem(&members, (char *)xrole, (char *)xid, 0);
197
198        xmlFree(xid);
199        xmlFree(xrole);
200        xmlFree(xtype);
201#endif
202    } else if (xmlStrEqual(name, BAD_CAST "osm")) {
203        /* ignore */
204    } else if (xmlStrEqual(name, BAD_CAST "bound")) {
205        /* ignore */
206    } else {
207        fprintf(stderr, "%s: Unknown element name: %s\n", __FUNCTION__, name);
208    }
209}
210
211void EndElement(const xmlChar *name)
212{
213    if (xmlStrEqual(name, BAD_CAST "node")) {
214        struct osmNode  * storenode;
215        storenode = (struct osmNode *) calloc(1,sizeof(struct osmNode));
216        if (storenode==NULL)
217        {
218                fprintf(stderr,"out of memory\n");
219                exit(1);
220        }
221
222      if (node_lat < merc_lat_min)
223            node_lat = merc_lat_min;
224      else if (node_lat > merc_lat_max)
225            node_lat = merc_lat_max;
226
227        storenode->id = osm_id;
228        storenode->lat = node_lat * DEG_TO_RAD;
229        storenode->lon = node_lon * DEG_TO_RAD;
230
231        rb_insert(nodes_table, storenode);
232    } else if (xmlStrEqual(name, BAD_CAST "way")) {
233        struct osmWay  * storeway;
234        struct keyval *p;
235        int i;
236        storeway = (struct osmWay *) calloc(1,sizeof(struct osmWay)+sizeof(int)*countList(&nds));
237        if (storeway==NULL)
238        {
239                fprintf(stderr,"out of memory\n");
240                exit(1);
241        }
242        storeway->id = osm_id;
243        for( i=0, p = popItem(&nds); p; i++, p = popItem(&nds) )
244                storeway->nds[i] = strtol(p->value, NULL, 10);
245        storeway->nds[i] = 0;
246        if( i == 1 )
247          fprintf(stderr, "Wierd: way %d only has %d nodes\n", osm_id, i);
248        if( i >= 2 )
249                rb_insert( ways_table, storeway );
250        resetList(&nds);
251    } else if (xmlStrEqual(name, BAD_CAST "relation")) {
252        /* ignore */
253    } else if (xmlStrEqual(name, BAD_CAST "tag")) {
254        /* ignore */
255    } else if (xmlStrEqual(name, BAD_CAST "nd")) {
256        /* ignore */
257    } else if (xmlStrEqual(name, BAD_CAST "member")) {
258        /* ignore */
259    } else if (xmlStrEqual(name, BAD_CAST "osm")) {
260        printStatus();
261    } else if (xmlStrEqual(name, BAD_CAST "bound")) {
262        /* ignore */
263    } else {
264        fprintf(stderr, "%s: Unknown element name: %s\n", __FUNCTION__, name);
265    }
266}
267
268static void processNode(xmlTextReaderPtr reader) {
269    xmlChar *name;
270    name = xmlTextReaderName(reader);
271    if (name == NULL)
272        name = xmlStrdup(BAD_CAST "--");
273       
274    switch(xmlTextReaderNodeType(reader)) {
275        case XML_READER_TYPE_ELEMENT:
276            StartElement(reader, name);
277            if (xmlTextReaderIsEmptyElement(reader))
278                EndElement(name); /* No end_element for self closing tags! */
279            break;
280        case XML_READER_TYPE_END_ELEMENT:
281            EndElement(name);
282            break;
283        case XML_READER_TYPE_SIGNIFICANT_WHITESPACE:
284            /* Ignore */
285            break;
286        default:
287            fprintf(stderr, "Unknown node type %d\n", xmlTextReaderNodeType(reader));
288            break;
289    }
290
291    xmlFree(name);
292}
293
294static int streamFile(char *filename) {
295    xmlTextReaderPtr reader;
296    int ret = 0;
297
298    reader = inputUTF8(filename);
299
300    if (reader != NULL) {
301        ret = xmlTextReaderRead(reader);
302        while (ret == 1) {
303            processNode(reader);
304            ret = xmlTextReaderRead(reader);
305        }
306
307        if (ret != 0) {
308            fprintf(stderr, "%s : failed to parse\n", filename);
309            return ret;
310        }
311
312        xmlFreeTextReader(reader);
313    } else {
314        fprintf(stderr, "Unable to open %s\n", filename);
315        return 1;
316    }
317    return 0;
318}
319
320void exit_nicely(void)
321{
322    fprintf(stderr, "Error occurred, cleaning up\n");
323    exit(1);
324}
325 
326static void usage(const char *arg0)
327{
328    const char *name = basename(arg0);
329
330    fprintf(stderr, "Usage:\n");
331    fprintf(stderr, "\t%s coastline.txt coastline.osm.gz outputprefix\n", name);
332    fprintf(stderr, "\n");
333}
334
335/* Can be used to compare any osm object, cause it only uses the first two fields */
336static int osm_compare (const void *pa, const void *pb, void *param){
337        const struct osmNode * na=pa;
338        const struct osmNode * nb=pb;
339        if (na->id < nb->id) return -1;
340        if (na->id > nb->id) return 1;
341        return 0;
342        param=param;
343}
344
345static int processList(char *filename, char *output_prefix)
346{
347  char out_c[64], out_p[64], out_i[64];
348  FILE *in;
349  SHPHandle shp_poly, shp_arc, shp_point;
350  DBFHandle dbf_poly, dbf_arc, dbf_point;
351 
352  in = fopen( filename, "rt" );
353  if( !in )
354  {
355      fprintf(stderr, "Error opening file %s: %s\n", filename, strerror(errno) );
356      exit_nicely();
357  }
358  if( strlen( output_prefix ) > sizeof(out_c)-10 )
359  {
360      fprintf( stderr, "Output prefix too long\n" );
361      exit_nicely();
362  }
363  sprintf( out_c, "%s_c", output_prefix );
364  sprintf( out_p, "%s_p", output_prefix );
365  sprintf( out_i, "%s_i", output_prefix );
366 
367  int max_vertex_count = 0;
368 
369  shp_poly = SHPCreate( out_c,  SHPT_POLYGON );
370  shp_arc = SHPCreate( out_i, SHPT_ARC );
371  shp_point = SHPCreate( out_p, SHPT_POINT );
372  dbf_poly = DBFCreate( out_c );
373  dbf_arc = DBFCreate( out_i );
374  dbf_point = DBFCreate( out_p );
375  DBFAddField( dbf_poly, "type", FTInteger, 5, 0 );
376  DBFAddField( dbf_poly, "length", FTInteger, 10, 0 );
377  DBFAddField( dbf_poly, "way_id", FTInteger, 10, 0 );
378  DBFAddField( dbf_arc, "type", FTInteger, 5, 0 );
379  DBFAddField( dbf_arc, "length", FTInteger, 10, 0 );
380  DBFAddField( dbf_arc, "way_id", FTInteger, 10, 0 );
381  DBFAddField( dbf_point, "type", FTInteger, 5, 0 );
382 
383  int shp_poly_count = 0, shp_arc_count = 0;
384 
385  double *v_x = malloc( MAX_VERTICES*sizeof(double) );
386  double *v_y = malloc( MAX_VERTICES*sizeof(double) );
387  double *v_z = malloc( MAX_VERTICES*sizeof(double) );
388  if( !v_x || !v_y || !v_z )
389  {
390    fprintf( stderr, "Out of memory allocating vertex buffers\n" );
391    exit_nicely();
392  }
393  int line = 0;
394  while( !feof(in) )
395  {
396    char type;
397    int length;
398    int i;
399    line++;
400    fscanf( in, "%c", &type );
401    if( type == 'P' )
402    {
403      int t;
404      double x, y;
405      double z = 0;
406      fscanf( in, "%d %lf %lf\n", &t, &y, &x );
407
408      // Coordinates very close to +/-180 sometimes get wrapped to the wrong side. Hack to fix...
409//      double sign = (fabs(x)>179.9) ? copysign( 1, x ) : 0;
410//      if( sign != 0 )
411//        fprintf( stderr, "Processing point (%f,%f)\n", x,y );
412      x *= DEG_TO_RAD;
413      y *= DEG_TO_RAD;
414//      printf( "Before transform (%f,%f)\n", x, y );
415      pj_transform(pj_ll, pj_merc, 1, 1, &x, &y, &z);
416//      if( sign != 0 )
417//        fprintf( stderr, "Projected to (%f,%f)\n", x, y );
418//      if( sign < 0 && x > 0 )
419//        sign -= 2*20037508.34;
420//      if( sign > 0 && x < 0 )
421//        sign += 2*20037508.34;
422     
423//      if( sign != 0 )
424//        fprintf( stderr, "Corrected to (%f,%f)\n", x, y );
425     
426//      printf( "After transform (%f,%f)\n", x, y );
427     
428//      if( x > 0.0 && y > 5527259.12027 && x < 78271.516964 && y < 5605266.15512 )
429//      {
430//        fprintf(stderr, "Found at line %d\n", line );
431//      }
432      SHPObject *p;
433      p = SHPCreateSimpleObject( SHPT_POINT, 1, &x, &y, NULL );
434      int idx = SHPWriteObject( shp_point, -1, p );
435      if( idx < 0 ) { fprintf(stderr, "Write failure: %m\n"); exit(1); }
436      DBFWriteIntegerAttribute( dbf_point, idx, 0, t );
437      SHPDestroyObject(p);
438     
439      continue;
440    }
441    if( type != 'C' && type != 'I' )
442    {
443      fprintf( stderr, "Got bad type at offset %ld\n", ftell(in) );
444      exit_nicely();
445    }
446    fscanf( in, "%d ", &length );
447    printf("Generating type=%c, length=%d, poly_count=%d, arc_count=%d\n", type, length, shp_poly_count, shp_arc_count );
448    struct osmNode *last_match = NULL;
449    struct osmNode *first_match = NULL;
450    int vertex_count = 0;
451    int way_id = 0;
452    for( i=0; i<length; i++ )
453    {
454      struct osmWay key_way;
455      struct osmNode key_node;
456      int j;
457      if( fscanf( in, "%d ", &key_way.id ) != 1 )
458      {
459        fprintf( stderr, "Failed to read\n");
460        exit_nicely();
461      } 
462      struct osmWay *way_match = rb_find( ways_table, &key_way );
463      if( !way_match )
464      {
465        fprintf( stderr, "Failed to find way %d\n", key_way.id );
466        continue;
467      }
468      if( way_id == 0 )
469        way_id = key_way.id;
470       
471      for( j=0; way_match->nds[j]; j++ )
472      {
473//        printf("Adding way %d, node %d (%d)\n", key_way.id, j, way_match->nds[j] );
474        key_node.id = way_match->nds[j];
475        struct osmNode *node_match = rb_find( nodes_table, &key_node );
476        if( !node_match )
477        {
478          fprintf( stderr, "Failed to find node %d\n", key_node.id );
479          continue;
480        }
481        if( !first_match )
482          first_match = node_match;
483        if( last_match == node_match )
484          continue;
485        if( last_match && fabs( last_match->lat - node_match->lat ) < 3e-6*DEG_TO_RAD
486                       && fabs( last_match->lon - node_match->lon ) < 3e-6*DEG_TO_RAD )
487        {
488#if 0
489          if( way_id == 22650041 )
490            fprintf(stderr, "Skipped (too close) {%.6f,%.6f} ", node_match->lat, node_match->lon );
491#endif
492          continue;
493        }
494        /* Extend array as necessary */
495        if( vertex_count > MAX_VERTICES-10 )
496        {
497          MAX_VERTICES <<= 1;
498          v_x = realloc( v_x, MAX_VERTICES*sizeof(double) );
499          v_y = realloc( v_y, MAX_VERTICES*sizeof(double) );
500          v_z = realloc( v_z, MAX_VERTICES*sizeof(double) );
501          memset( v_z, 0, MAX_VERTICES*sizeof(double) );
502          fprintf( stderr, "Resized vertex arrays to %d\n", MAX_VERTICES );
503        }
504        v_y[vertex_count] = node_match->lat;
505        v_x[vertex_count] = node_match->lon;
506       
507        if( vertex_count > 1 )
508        {
509          if( fabs( v_y[vertex_count] - v_y[vertex_count-1] ) > 1*DEG_TO_RAD ||
510              fabs( v_x[vertex_count] - v_x[vertex_count-1] ) > 1*DEG_TO_RAD )
511          {
512            fprintf( stderr, "Problem found at way %d, node %d (%d) (%.3f,%.3f %.3f,%.3f)\n", 
513                      way_id, j, way_match->nds[j],
514                     v_x[vertex_count-1]/DEG_TO_RAD, v_y[vertex_count-1]/DEG_TO_RAD, 
515                     v_x[vertex_count]  /DEG_TO_RAD, v_y[vertex_count]  /DEG_TO_RAD );
516          }
517        }
518        printf( "(%.6f,%.6f) ", node_match->lat, node_match->lon );
519        vertex_count++;
520        last_match = node_match;
521      }
522      printf("|");
523    }
524    printf("\n");
525    if( type == 'C' && vertex_count < MAX_VERTICES-1 ) // Make sure polygon is closed
526    {
527      if( last_match != first_match )
528      {
529        if( fabs( last_match->lat - first_match->lat ) < 1e-6 &&
530            fabs( last_match->lon - first_match->lon ) < 1e-6 )
531        {
532          // If it's close, make sure it's equal
533          v_y[vertex_count-1] = first_match->lat;
534          v_x[vertex_count-1] = first_match->lon;
535        }
536        else
537        {
538          if( fabs( last_match->lat - first_match->lat ) > 0.5 ||
539              fabs( last_match->lon - first_match->lon ) > 0.5 )
540          {
541            fprintf( stderr, "Warning: Overly long last leg: %d\n", way_id );
542          }
543          // Otherwise, link to the first node
544          v_y[vertex_count] = first_match->lat;
545          v_x[vertex_count] = first_match->lon;
546          vertex_count++;
547        }
548      }
549    }
550    if( vertex_count > max_vertex_count )
551      max_vertex_count = vertex_count;
552
553#if 0
554    if( way_id == 20806670 )
555    {
556      fprintf( stderr, "Before:\n");
557      for(int y=0; y<vertex_count; y++)
558        fprintf( stderr, "(%f,%f) ", v_x[y]/DEG_TO_RAD,v_y[y]/DEG_TO_RAD);
559      fprintf(stderr, "\n");
560    }
561#endif
562    memset( v_z, 0, sizeof(double) * ((vertex_count < MAX_VERTICES) ? vertex_count : MAX_VERTICES) );
563    pj_transform(pj_ll, pj_merc, (vertex_count < MAX_VERTICES) ? vertex_count : MAX_VERTICES, 1, v_x, v_y, v_z);
564     
565    SHPObject *tmp = SHPCreateSimpleObject( (type == 'C') ? SHPT_POLYGON : SHPT_ARC, 
566                                            (vertex_count < MAX_VERTICES) ? vertex_count : MAX_VERTICES, 
567                                            v_x, v_y, NULL );
568    if( (tmp->dfXMax - tmp->dfXMin) > 30000000 )
569    {
570      fprintf( stderr, "Oversize: way_id=%d\n", way_id );
571      int scoreleft = 0, scoreright = 0;
572      for(int y=0; y<vertex_count; y++)
573      {
574        if( v_x[y] < -20000000 ) scoreleft++;
575        if( v_x[y] > +20000000 ) scoreright++;
576      }
577      for(int y=0; y<vertex_count; y++)
578      {
579        if( scoreleft < scoreright && v_x[y] < -20000000 ) v_x[y] += 2*20037508.34f;
580        if( scoreleft > scoreright && v_x[y] > +20000000 ) v_x[y] -= 2*20037508.34f;
581      }
582      // Recreate object
583      SHPDestroyObject(tmp);
584      tmp = SHPCreateSimpleObject( (type == 'C') ? SHPT_POLYGON : SHPT_ARC, 
585                                   (vertex_count < MAX_VERTICES) ? vertex_count : MAX_VERTICES, 
586                                   v_x, v_y, NULL );
587      if( (tmp->dfXMax - tmp->dfXMin) > 30000000 )
588      {
589        fprintf( stderr, "Still Oversize: way_id=%d\n", way_id );
590        for( int y=0; y<vertex_count; y++)
591          fprintf( stderr, "(%f,%f), ", v_x[y], v_y[y] );
592        fprintf(stderr, "\n");
593      }
594    }
595#if 0
596    if( way_id == 20806670 )
597    {
598      fprintf( stderr, "After:\n");
599      for(int y=0; y<vertex_count; y++)
600        fprintf( stderr, "(%f,%f) ", v_x[y],v_y[y]);
601      fprintf(stderr, "\n");
602    }
603#endif
604    if( vertex_count < 2 )
605    {
606      fprintf(stderr, "Way %d only has %d node(s)\n", way_id, vertex_count);
607    }
608    if( type == 'C' && vertex_count >= 4 )
609    {
610      int idx = SHPWriteObject( shp_poly, -1, tmp );
611      if( idx < 0 ) { fprintf(stderr, "Write failure: %m\n"); exit(1); }
612      shp_poly_count++;
613      DBFWriteIntegerAttribute( dbf_poly, idx, 0, 0 );
614      DBFWriteIntegerAttribute( dbf_poly, idx, 1, vertex_count );
615      DBFWriteIntegerAttribute( dbf_poly, idx, 2, way_id );
616    }
617    else if( vertex_count > 0 )
618    {
619      tmp->nSHPType = SHPT_ARC;
620      int idx = SHPWriteObject( shp_arc, -1, tmp );
621      if( idx < 0 ) { fprintf(stderr, "Write failure: %m\n"); exit(1); }
622      shp_arc_count++;
623      DBFWriteIntegerAttribute( dbf_arc, idx, 0, 1 );
624      DBFWriteIntegerAttribute( dbf_arc, idx, 1, vertex_count );
625      DBFWriteIntegerAttribute( dbf_arc, idx, 2, way_id );
626     
627      SHPObject *p;
628      p = SHPCreateSimpleObject( SHPT_POINT, 1, &v_x[0], &v_y[0], NULL );
629      idx = SHPWriteObject( shp_point, -1, p );
630      if( idx < 0 ) { fprintf(stderr, "Write failure: %m\n"); exit(1); }
631      DBFWriteIntegerAttribute( dbf_point, idx, 0, 2 );
632      SHPDestroyObject(p);
633     
634      if( type != 'C' )
635      {
636        p = SHPCreateSimpleObject( SHPT_POINT, 1, &v_x[vertex_count-1], &v_y[vertex_count-1], NULL );
637        idx = SHPWriteObject( shp_point, -1, p );
638        if( idx < 0 ) { fprintf(stderr, "Write failure: %m\n"); exit(1); }
639        DBFWriteIntegerAttribute( dbf_point, idx, 0, 2 );
640        SHPDestroyObject(p);
641      }
642    }
643     
644    SHPDestroyObject(tmp);
645  }
646  fprintf( stderr, "Max vertex count: %d\n", max_vertex_count );
647  fprintf( stderr, "Polygons: %d, Arcs: %d\n", shp_poly_count, shp_arc_count );
648  if( max_vertex_count >= MAX_VERTICES )
649    fprintf( stderr, "objects cropped to %d vertices\n", MAX_VERTICES );
650  SHPCreateTree( shp_poly, 2, 10, NULL, NULL );
651  SHPCreateTree( shp_arc, 2, 10, NULL, NULL );
652  SHPCreateTree( shp_point, 2, 10, NULL, NULL );
653  SHPClose( shp_poly );
654  SHPClose( shp_arc );
655  SHPClose( shp_point );
656  DBFClose( dbf_poly );
657  DBFClose( dbf_arc );
658  DBFClose( dbf_point );
659  return 0;
660}
661
662int main(int argc, char *argv[])
663{
664    fprintf(stderr, "coast2shp SVN version %s $Rev: 4895 $ \n", VERSION);
665
666    if (argc != 4) {
667        usage(argv[0]);
668        exit(EXIT_FAILURE);
669    }
670
671    initList(&nds);
672    nodes_table = rb_create (osm_compare, NULL, NULL);
673    ways_table = rb_create (osm_compare, NULL, NULL);
674
675    count_node = max_node = 0;
676    count_way = max_way = 0;
677    count_rel = max_rel = 0;
678
679    project_init();
680    LIBXML_TEST_VERSION
681
682    fprintf(stderr, "\nReading in file: %s\n", argv[2]);
683    if (streamFile(argv[2]) != 0)
684        exit_nicely();
685
686    xmlCleanupParser();
687    xmlMemoryDump();
688
689    fprintf(stderr, "\nReading in file: %s\n", argv[1]);
690    processList(argv[1], argv[3]);
691   
692    if (count_node || count_way || count_rel) {
693        fprintf(stderr, "\n");
694        fprintf(stderr, "Node stats: total(%d), max(%d)\n", count_node, max_node);
695        fprintf(stderr, "Way stats: total(%d), max(%d)\n", count_way, max_way);
696        fprintf(stderr, "Relation stats: total(%d), max(%d)\n", count_rel, max_rel);
697    }
698
699    project_exit();
700
701    return 0;
702}
703
Note: See TracBrowser for help on using the repository browser.