← Back to branch summary

~ubuntu-branches/ubuntu/vivid/grass/vivid-proposed

« back to all changes in this revision

Viewing changes to lib/vector/Vlib/build_sfa.c

  • Committer: Package Import Robot
  • Author(s): Bas Couwenberg
  • Date: 2015-02-20 23:12:08 UTC
  • mfrom: (8.2.6 experimental)
  • Revision ID: package-import@ubuntu.com-20150220231208-1u6qvqm84v430b10
Tags: 7.0.0-1~exp1
* New upstream release.
* Update python-ctypes-ternary.patch to use if/else instead of and/or.
* Drop check4dev patch, rely on upstream check.
* Add build dependency on libpq-dev to grass-dev for libpq-fe.h.
* Drop patches applied upstream, refresh remaining patches.
* Update symlinks for images switched from jpg to png.

Show diffs side-by-side

added added

removed removed

Lines of Context:
lib/vector/Vlib/build_sfa.c
 
1
/*!
 
2
   \file lib/vector/Vlib/build_sfa.c
 
3
 
 
4
   \brief Vector library - Building pseudo-topology for simple feature access
 
5
 
 
6
   Higher level functions for reading/writing/manipulating vectors.
 
7
 
 
8
   Line offset is
 
9
    - centroids   : FID
 
10
    - other types : index of the first record (which is FID) in offset array.
 
11
 
 
12
   (C) 2001-2012 by the GRASS Development Team
 
13
 
 
14
   This program is free software under the GNU General Public License
 
15
   (>=v2). Read the file COPYING that comes with GRASS for details.
 
16
 
 
17
   \author Radim Blazek
 
18
   \author Piero Cavalieri
 
19
   \author Various updates for GRASS 7 by Martin Landa <landa.martin gmail.com>
 
20
*/
 
21
 
 
22
#include <stdlib.h>
 
23
 
 
24
#include <grass/gis.h>
 
25
#include <grass/vector.h>
 
26
#include <grass/glocale.h>
 
27
 
 
28
/*!  
 
29
  \brief This structure keeps info about geometry parts above current
 
30
  geometry, path to curent geometry in the feature. First 'part' number
 
31
  however is feature id */
 
32
struct geom_parts
 
33
{
 
34
    int *part;
 
35
    int a_parts;
 
36
    int n_parts;
 
37
};
 
38
 
 
39
static void init_parts(struct geom_parts *);
 
40
static void reset_parts(struct geom_parts *);
 
41
static void free_parts(struct geom_parts *);
 
42
static void add_part(struct geom_parts *, int);
 
43
static void del_part(struct geom_parts *);
 
44
static void add_parts_to_offset(struct Format_info_offset *,
 
45
                                struct geom_parts *);
 
46
static int add_line(struct Plus_head *, struct Format_info_offset *,
 
47
                    int, struct line_pnts *,
 
48
                    int, struct geom_parts *);
 
49
 
 
50
#ifdef HAVE_POSTGRES
 
51
#include "pg_local_proto.h"
 
52
 
 
53
static int add_geometry_pg(struct Plus_head *,
 
54
                           struct Format_info_pg *,
 
55
                           struct feat_parts *, int, 
 
56
                           int, int,
 
57
 struct geom_parts *);
 
58
static void build_pg(struct Map_info *, int);
 
59
#endif
 
60
 
 
61
#ifdef HAVE_OGR
 
62
#include <ogr_api.h>
 
63
 
 
64
static int add_geometry_ogr(struct Plus_head *,
 
65
                            struct Format_info_ogr *,
 
66
                            OGRGeometryH, int, int,
 
67
                            struct geom_parts *);
 
68
 
 
69
static void build_ogr(struct Map_info *, int);
 
70
#endif
 
71
 
 
72
/*!
 
73
  \brief Init parts
 
74
*/
 
75
void init_parts(struct geom_parts * parts)
 
76
{
 
77
    G_zero(parts, sizeof(struct geom_parts));
 
78
}
 
79
 
 
80
/*!
 
81
  \brief Reset parts
 
82
*/
 
83
void reset_parts(struct geom_parts * parts)
 
84
{
 
85
    parts->n_parts = 0;
 
86
}
 
87
 
 
88
/*!
 
89
  \brief Free parts
 
90
*/
 
91
void free_parts(struct geom_parts * parts)
 
92
{
 
93
    G_free(parts->part);
 
94
    G_zero(parts, sizeof(struct geom_parts));
 
95
}
 
96
 
 
97
/*!
 
98
  \brief Add new part number to parts
 
99
*/
 
100
void add_part(struct geom_parts *parts, int part)
 
101
{
 
102
    if (parts->a_parts == parts->n_parts) {
 
103
        parts->a_parts += 10;
 
104
        parts->part = (int *) G_realloc((void *)parts->part,
 
105
                                        parts->a_parts * sizeof(int));
 
106
    }
 
107
    parts->part[parts->n_parts] = part;
 
108
    parts->n_parts++;
 
109
}
 
110
 
 
111
/*!
 
112
  \brief Remove last part
 
113
*/
 
114
void del_part(struct geom_parts *parts)
 
115
{
 
116
    parts->n_parts--;
 
117
}
 
118
 
 
119
/*!
 
120
  \brief Add parts to offset
 
121
*/
 
122
void add_parts_to_offset(struct Format_info_offset *offset,
 
123
                         struct geom_parts *parts)
 
124
{
 
125
    int i, j;
 
126
 
 
127
    if (offset->array_num + parts->n_parts >= offset->array_alloc) {
 
128
        offset->array_alloc += parts->n_parts + 1000;
 
129
        offset->array = (int *)G_realloc(offset->array,
 
130
                                         offset->array_alloc * sizeof(int));
 
131
    }
 
132
    j = offset->array_num;
 
133
    for (i = 0; i < parts->n_parts; i++) {
 
134
        G_debug(4, "add offset %d", parts->part[i]);
 
135
        offset->array[j] = parts->part[i];
 
136
        j++;
 
137
    }
 
138
    offset->array_num += parts->n_parts;
 
139
}
 
140
 
 
141
/*!
 
142
  \brief Add line to support structures
 
143
*/
 
144
int add_line(struct Plus_head *plus, struct Format_info_offset *offset,
 
145
             int type, struct line_pnts *Points,
 
146
             int FID, struct geom_parts *parts)
 
147
{
 
148
    int line;
 
149
    long offset_value;
 
150
    struct bound_box box;
 
151
 
 
152
    if (type != GV_CENTROID) {
 
153
        /* beginning in the offset array */
 
154
        offset_value = offset->array_num;
 
155
    }
 
156
    else {
 
157
        /* TODO : could be used to statore category ? */
 
158
        /* because centroids are read from topology, not from layer */
 
159
        offset_value = FID; 
 
160
    }
 
161
 
 
162
    G_debug(4, "Register line: FID = %d offset = %ld", FID, offset_value);
 
163
    dig_line_box(Points, &box);
 
164
    line = dig_add_line(plus, type, Points, &box, offset_value);
 
165
    G_debug(4, "Line registered with line = %d", line);
 
166
 
 
167
    /* Set box */
 
168
    if (line == 1)
 
169
        Vect_box_copy(&(plus->box), &box);
 
170
    else
 
171
        Vect_box_extend(&(plus->box), &box);
 
172
 
 
173
    if (type != GV_BOUNDARY) {
 
174
        dig_cidx_add_cat(plus, 1, (int)FID, line, type);
 
175
    }
 
176
    else {
 
177
        dig_cidx_add_cat(plus, 0, 0, line, type);
 
178
    }
 
179
 
 
180
    /* because centroids are read from topology, not from layer */
 
181
    if (type != GV_CENTROID)
 
182
        add_parts_to_offset(offset, parts);
 
183
    
 
184
    return line;
 
185
}
 
186
 
 
187
#ifdef HAVE_POSTGRES
 
188
/*!
 
189
  \brief Recursively add geometry (PostGIS) to topology
 
190
*/
 
191
int add_geometry_pg(struct Plus_head *plus,
 
192
                    struct Format_info_pg *pg_info,
 
193
                    struct feat_parts *fparts, int ipart,
 
194
                    int FID, int build, struct geom_parts *parts)
 
195
{
 
196
    int line, i, idx, area, isle, outer_area, ret;
 
197
    int lines[1];
 
198
    double area_size, x, y;
 
199
    SF_FeatureType ftype;
 
200
    struct bound_box box;
 
201
    struct Format_info_offset *offset;
 
202
    struct line_pnts *line_i;
 
203
    
 
204
    ftype = fparts->ftype[ipart];
 
205
    
 
206
    G_debug(4, "add_geometry_pg() FID = %d ftype = %d", FID, ftype);
 
207
    
 
208
    offset = &(pg_info->offset);
 
209
    
 
210
    outer_area = 0;
 
211
    
 
212
    switch(ftype) {     
 
213
    case SF_POINT:
 
214
        G_debug(4, "Point");
 
215
        line_i = pg_info->cache.lines[fparts->idx[ipart]];
 
216
        add_line(plus, offset, GV_POINT, line_i,
 
217
                 FID, parts);
 
218
        break;
 
219
    case SF_LINESTRING:
 
220
        G_debug(4, "LineString");
 
221
        line_i = pg_info->cache.lines[fparts->idx[ipart]];
 
222
        add_line(plus, offset, GV_LINE, line_i,
 
223
                 FID, parts);
 
224
        break;
 
225
    case SF_POLYGON:
 
226
        G_debug(4, "Polygon");
 
227
 
 
228
        /* register boundaries */
 
229
        idx = fparts->idx[ipart];
 
230
        for (i = 0; i < fparts->nlines[ipart]; i++) {
 
231
            line_i = pg_info->cache.lines[idx++];
 
232
            G_debug(4, "part %d", i);
 
233
            add_part(parts, i);
 
234
            line = add_line(plus, offset, GV_BOUNDARY,
 
235
                            line_i, FID, parts);
 
236
            del_part(parts);
 
237
 
 
238
            if (build < GV_BUILD_AREAS)
 
239
                continue;
 
240
            
 
241
            /* add area (each inner ring is also area) */
 
242
            dig_line_box(line_i, &box);
 
243
            dig_find_area_poly(line_i, &area_size);
 
244
 
 
245
            if (area_size > 0)          /* area clockwise */
 
246
                lines[0] = line;
 
247
            else
 
248
                lines[0] = -line;
 
249
 
 
250
            area = dig_add_area(plus, 1, lines, &box);
 
251
 
 
252
            /* each area is also isle */
 
253
            lines[0] = -lines[0];       /* island is counter clockwise */
 
254
 
 
255
            isle = dig_add_isle(plus, 1, lines, &box);
 
256
 
 
257
            if (build < GV_BUILD_ATTACH_ISLES)
 
258
                continue;
 
259
            
 
260
            if (i == 0) {       /* outer ring */
 
261
                outer_area = area;
 
262
            }
 
263
            else {              /* inner ring */
 
264
                struct P_isle *Isle;
 
265
 
 
266
                Isle = plus->Isle[isle];
 
267
                Isle->area = outer_area;
 
268
 
 
269
                dig_area_add_isle(plus, outer_area, isle);
 
270
            }
 
271
        }
 
272
 
 
273
        if (build >= GV_BUILD_CENTROIDS) {
 
274
            /* create virtual centroid */
 
275
            ret = Vect_get_point_in_poly_isl((const struct line_pnts *) pg_info->cache.lines[fparts->idx[ipart]],
 
276
                                             (const struct line_pnts **) &pg_info->cache.lines[fparts->idx[ipart]] + 1,
 
277
                                             fparts->nlines[ipart] - 1, &x, &y);
 
278
            if (ret < -1) {
 
279
                G_warning(_("Unable to calculate centroid for area %d"),
 
280
                          outer_area);
 
281
            }
 
282
            else {
 
283
                struct P_area *Area;
 
284
                struct P_topo_c *topo;
 
285
                struct P_line *Line;
 
286
                struct line_pnts *line_c;
 
287
 
 
288
                G_debug(4, "  Centroid: %f, %f", x, y);
 
289
                line_c = Vect_new_line_struct();
 
290
                Vect_append_point(line_c, x, y, 0.0);
 
291
                line = add_line(plus, offset, GV_CENTROID, line_c, FID, parts);
 
292
                
 
293
                Line = plus->Line[line];
 
294
                topo = (struct P_topo_c *)Line->topo;
 
295
                topo->area = outer_area;
 
296
                
 
297
                /* register centroid to area */
 
298
                Area = plus->Area[outer_area];
 
299
                Area->centroid = line;
 
300
                Vect_destroy_line_struct(line_c);
 
301
            }
 
302
        }
 
303
        break;
 
304
    default:
 
305
        G_warning(_("Feature type %d not supported"), ftype);
 
306
        break;
 
307
    }
 
308
    
 
309
    return 0;
 
310
}
 
311
 
 
312
/*!
 
313
  \brief Build pseudo-topology for PostGIS layers
 
314
*/
 
315
void build_pg(struct Map_info *Map, int build)
 
316
{
 
317
    int iFeature, ipart, fid, nrecords, npoints;
 
318
    char *wkb_data;
 
319
    
 
320
    struct Format_info_pg *pg_info;
 
321
    
 
322
    struct feat_parts fparts;
 
323
    struct geom_parts parts;
 
324
 
 
325
    pg_info = &(Map->fInfo.pg);
 
326
 
 
327
    /* initialize data structures */
 
328
    init_parts(&parts);
 
329
    G_zero(&fparts, sizeof(struct feat_parts));
 
330
    
 
331
    /* get all features */
 
332
    if (Vect__open_cursor_next_line_pg(pg_info, TRUE, Map->plus.built) != 0)
 
333
        return;
 
334
    
 
335
    /* scan records */
 
336
    npoints = 0;
 
337
    nrecords = PQntuples(pg_info->res);
 
338
    G_debug(4, "build_pg(): nrecords = %d", nrecords);
 
339
    G_message(_("Registering primitives..."));
 
340
    for (iFeature = 0; iFeature < nrecords; iFeature++) {
 
341
        /* get feature id */
 
342
        fid  = atoi(PQgetvalue(pg_info->res, iFeature, 1));
 
343
        if (fid < 1)
 
344
            continue; /* PostGIS Topology: skip features with negative
 
345
                       * fid (isles, universal face, ...) */
 
346
 
 
347
        wkb_data = PQgetvalue(pg_info->res, iFeature, 0);
 
348
        
 
349
        G_progress(iFeature + 1, 1e4);
 
350
 
 
351
        /* cache feature (lines) */
 
352
        if (SF_NONE == Vect__cache_feature_pg(wkb_data, FALSE, FALSE,
 
353
                                              &(pg_info->cache), &fparts)) {
 
354
            G_warning(_("Feature %d without geometry skipped"),
 
355
                      iFeature + 1);
 
356
            continue;
 
357
        }
 
358
        
 
359
        /* register all parts */
 
360
        reset_parts(&parts);
 
361
        add_part(&parts, fid);
 
362
        for (ipart = 0; ipart < fparts.n_parts; ipart++) {
 
363
            if (fparts.nlines[ipart] < 1) {
 
364
                G_warning(_("Feature %d without geometry skipped"), fid);
 
365
                continue;
 
366
            }
 
367
            
 
368
            npoints += pg_info->cache.lines[ipart]->n_points;
 
369
 
 
370
            G_debug(4, "Feature: fid = %d part = %d", fid, ipart);
 
371
            
 
372
            if (fparts.n_parts > 1)
 
373
                add_part(&parts, ipart);
 
374
            add_geometry_pg(&(Map->plus), pg_info, &fparts, ipart,
 
375
                            fid, build, &parts);
 
376
            if (fparts.n_parts > 1)
 
377
                del_part(&parts);
 
378
        }
 
379
 
 
380
        /* read next feature from cache */
 
381
        pg_info->cache.lines_next = 0;
 
382
    }
 
383
    G_progress(1, 1);
 
384
 
 
385
    G_message(_n("One primitive registered", "%d primitives registered", Map->plus.n_lines), Map->plus.n_lines);
 
386
    G_message(_n("One vertex registered", "%d vertices registered", npoints), npoints);
 
387
 
 
388
    Map->plus.built = GV_BUILD_BASE;
 
389
 
 
390
    PQclear(pg_info->res);
 
391
    pg_info->res = NULL;
 
392
    
 
393
    /* free allocated space */
 
394
    free_parts(&parts);
 
395
}
 
396
#endif /* HAVE_POSTGRES */
 
397
 
 
398
#ifdef HAVE_OGR
 
399
/*!
 
400
  \brief Recursively add geometry (OGR) to topology
 
401
*/
 
402
int add_geometry_ogr(struct Plus_head *plus,
 
403
                     struct Format_info_ogr *ogr_info,
 
404
                     OGRGeometryH hGeom, int FID, int build,
 
405
                     struct geom_parts *parts)
 
406
{
 
407
    int i, ret, npoints, line;
 
408
    int area, isle, outer_area;
 
409
    int lines[1];
 
410
    double area_size, x, y;
 
411
    int eType, nRings, iPart, nParts, nPoints;
 
412
 
 
413
    struct bound_box box;
 
414
    struct P_line *Line;
 
415
    struct Format_info_offset *offset;
 
416
    
 
417
    OGRGeometryH hGeom2, hRing;
 
418
 
 
419
    G_debug(4, "add_geometry_ogr() FID = %d", FID);
 
420
 
 
421
    offset = &(ogr_info->offset);
 
422
    
 
423
    /* allocate space in cache */
 
424
    if (!ogr_info->cache.lines) {
 
425
        ogr_info->cache.lines_alloc = 1;
 
426
        ogr_info->cache.lines = (struct line_pnts **) G_malloc(sizeof(struct line_pnts *));
 
427
        
 
428
        ogr_info->cache.lines_types = (int *) G_malloc(sizeof(int));
 
429
        ogr_info->cache.lines[0] = Vect_new_line_struct();
 
430
        ogr_info->cache.lines_types[0] = -1;
 
431
    }
 
432
 
 
433
    npoints = outer_area = 0;
 
434
    eType = wkbFlatten(OGR_G_GetGeometryType(hGeom));
 
435
    G_debug(4, "OGR type = %d", eType);
 
436
 
 
437
    switch (eType) {
 
438
    case wkbPoint:
 
439
        G_debug(4, "Point");
 
440
        
 
441
        ogr_info->cache.lines_types[0] = GV_POINT;
 
442
        Vect_reset_line(ogr_info->cache.lines[0]);
 
443
        Vect_append_point(ogr_info->cache.lines[0], OGR_G_GetX(hGeom, 0),
 
444
                          OGR_G_GetY(hGeom, 0), OGR_G_GetZ(hGeom, 0));
 
445
        add_line(plus, offset, GV_POINT, ogr_info->cache.lines[0], FID, parts);
 
446
        npoints += ogr_info->cache.lines[0]->n_points;
 
447
        break;
 
448
 
 
449
    case wkbLineString:
 
450
        G_debug(4, "LineString");
 
451
        
 
452
        ogr_info->cache.lines_types[0] = GV_LINE;
 
453
        nPoints = OGR_G_GetPointCount(hGeom);
 
454
        Vect_reset_line(ogr_info->cache.lines[0]);
 
455
        for (i = 0; i < nPoints; i++) {
 
456
            Vect_append_point(ogr_info->cache.lines[0],
 
457
                              OGR_G_GetX(hGeom, i), OGR_G_GetY(hGeom, i),
 
458
                              OGR_G_GetZ(hGeom, i));
 
459
        }
 
460
        add_line(plus, offset, GV_LINE, ogr_info->cache.lines[0], FID, parts);  
 
461
        npoints += ogr_info->cache.lines[0]->n_points;
 
462
        break;
 
463
 
 
464
    case wkbPolygon:
 
465
        G_debug(4, "Polygon");
 
466
 
 
467
        nRings = OGR_G_GetGeometryCount(hGeom);
 
468
        G_debug(4, "Number of rings: %d", nRings);
 
469
 
 
470
        /* alloc space for islands if needed */
 
471
        if (nRings > ogr_info->cache.lines_alloc) {
 
472
            ogr_info->cache.lines_alloc += nRings;
 
473
            ogr_info->cache.lines = (struct line_pnts **) G_realloc(ogr_info->cache.lines,
 
474
                                                                    ogr_info->cache.lines_alloc *
 
475
                                                                    sizeof(struct line_pnts *));
 
476
            ogr_info->cache.lines_types = (int *) G_realloc(ogr_info->cache.lines_types,
 
477
                                                            ogr_info->cache.lines_alloc * sizeof(int));
 
478
                
 
479
            for (i = ogr_info->cache.lines_alloc - nRings; i < ogr_info->cache.lines_alloc; i++) {
 
480
                ogr_info->cache.lines[i] = Vect_new_line_struct();
 
481
                ogr_info->cache.lines_types[i] = -1;
 
482
            }
 
483
        }
 
484
 
 
485
        /* go thru rings */
 
486
        for (iPart = 0; iPart < nRings; iPart++) {
 
487
            ogr_info->cache.lines_types[iPart] = GV_BOUNDARY;
 
488
            hRing = OGR_G_GetGeometryRef(hGeom, iPart);
 
489
            nPoints = OGR_G_GetPointCount(hRing);
 
490
            G_debug(4, "  ring %d : nPoints = %d", iPart, nPoints);
 
491
 
 
492
            Vect_reset_line(ogr_info->cache.lines[iPart]);
 
493
            for (i = 0; i < nPoints; i++) {
 
494
                Vect_append_point(ogr_info->cache.lines[iPart],
 
495
                                  OGR_G_GetX(hRing, i), OGR_G_GetY(hRing, i),
 
496
                                  OGR_G_GetZ(hRing, i));
 
497
            }
 
498
            npoints += ogr_info->cache.lines[iPart]->n_points;
 
499
            
 
500
            /* register boundary */
 
501
            add_part(parts, iPart);
 
502
            line = add_line(plus, offset, GV_BOUNDARY, ogr_info->cache.lines[iPart], FID, parts);
 
503
            del_part(parts);
 
504
 
 
505
            if (build < GV_BUILD_AREAS)
 
506
                continue;
 
507
            
 
508
            /* add area (each inner ring is also area) */
 
509
            dig_line_box(ogr_info->cache.lines[iPart], &box);
 
510
            dig_find_area_poly(ogr_info->cache.lines[iPart], &area_size);
 
511
 
 
512
            if (area_size > 0)          /* area clockwise */
 
513
                lines[0] = line;
 
514
            else
 
515
                lines[0] = -line;
 
516
 
 
517
            area = dig_add_area(plus, 1, lines, &box);
 
518
 
 
519
            /* each area is also isle */
 
520
            lines[0] = -lines[0];       /* island is counter clockwise */
 
521
 
 
522
            isle = dig_add_isle(plus, 1, lines, &box);
 
523
 
 
524
            if (build < GV_BUILD_ATTACH_ISLES)
 
525
                continue;
 
526
            
 
527
            if (iPart == 0) {   /* outer ring */
 
528
                outer_area = area;
 
529
            }
 
530
            else {              /* inner ring */
 
531
                struct P_isle *Isle;
 
532
 
 
533
                Isle = plus->Isle[isle];
 
534
                Isle->area = outer_area;
 
535
 
 
536
                dig_area_add_isle(plus, outer_area, isle);
 
537
            }
 
538
        }
 
539
        
 
540
        if (build >= GV_BUILD_CENTROIDS) {
 
541
            /* create virtual centroid */
 
542
            ret = Vect_get_point_in_poly_isl((const struct line_pnts *) ogr_info->cache.lines[0],
 
543
                                             (const struct line_pnts **) ogr_info->cache.lines + 1,
 
544
                                             nRings - 1, &x, &y);
 
545
            if (ret < -1) {
 
546
                G_warning(_("Unable to calculate centroid for area %d"),
 
547
                          outer_area);
 
548
            }
 
549
            else {
 
550
                struct P_area *Area;
 
551
                struct P_topo_c *topo;
 
552
                
 
553
                G_debug(4, "  Centroid: %f, %f", x, y);
 
554
                Vect_reset_line(ogr_info->cache.lines[0]);
 
555
                Vect_append_point(ogr_info->cache.lines[0], x, y, 0.0);
 
556
                line = add_line(plus, offset, GV_CENTROID, ogr_info->cache.lines[0],
 
557
                                FID, parts);
 
558
                
 
559
                Line = plus->Line[line];
 
560
                topo = (struct P_topo_c *)Line->topo;
 
561
                topo->area = outer_area;
 
562
                
 
563
                /* register centroid to area */
 
564
                Area = plus->Area[outer_area];
 
565
                Area->centroid = line;
 
566
            }
 
567
        }
 
568
        break;
 
569
 
 
570
    case wkbMultiPoint:
 
571
    case wkbMultiLineString:
 
572
    case wkbMultiPolygon:
 
573
    case wkbGeometryCollection:
 
574
        nParts = OGR_G_GetGeometryCount(hGeom);
 
575
        G_debug(4, "%d geoms -> next level", nParts);
 
576
 
 
577
        /* alloc space for parts if needed */
 
578
        if (nParts > ogr_info->cache.lines_alloc) {
 
579
            ogr_info->cache.lines_alloc += nParts;
 
580
            ogr_info->cache.lines = (struct line_pnts **) G_realloc(ogr_info->cache.lines,
 
581
                                                                    ogr_info->cache.lines_alloc *
 
582
                                                                    sizeof(struct line_pnts *));
 
583
            ogr_info->cache.lines_types = (int *) G_realloc(ogr_info->cache.lines_types,
 
584
                                                            ogr_info->cache.lines_alloc * sizeof(int));
 
585
                
 
586
            for (i = ogr_info->cache.lines_alloc - nParts; i < ogr_info->cache.lines_alloc; i++) {
 
587
                ogr_info->cache.lines[i] = Vect_new_line_struct();
 
588
                ogr_info->cache.lines_types[i] = -1;
 
589
            }
 
590
        }
 
591
        
 
592
        /* go thru all parts */
 
593
        for (i = 0; i < nParts; i++) {
 
594
            add_part(parts, i);
 
595
            hGeom2 = OGR_G_GetGeometryRef(hGeom, i);
 
596
            npoints += add_geometry_ogr(plus, ogr_info, hGeom2,
 
597
                                        FID, build, parts);
 
598
            del_part(parts);
 
599
        }
 
600
        break;
 
601
 
 
602
    default:
 
603
        G_warning(_("OGR feature type %d not supported"), eType);
 
604
        break;
 
605
    }
 
606
 
 
607
    return npoints;
 
608
}
 
609
 
 
610
void build_ogr(struct Map_info *Map, int build)
 
611
{
 
612
    int iFeature, FID, npoints, nskipped;
 
613
    
 
614
    struct Format_info_ogr *ogr_info;
 
615
    
 
616
    OGRFeatureH hFeature;
 
617
    OGRGeometryH hGeom;
 
618
                
 
619
    struct geom_parts parts;
 
620
    
 
621
    ogr_info = &(Map->fInfo.ogr);
 
622
    
 
623
    /* initialize data structures */
 
624
    init_parts(&parts);
 
625
 
 
626
    /* Note: Do not use OGR_L_GetFeatureCount (it may scan all features) */
 
627
    OGR_L_ResetReading(ogr_info->layer);
 
628
    npoints = iFeature = nskipped = 0;
 
629
    G_message(_("Registering primitives..."));
 
630
    while ((hFeature = OGR_L_GetNextFeature(ogr_info->layer)) != NULL) {
 
631
        G_debug(3, "   Feature %d", iFeature);
 
632
        
 
633
        G_progress(++iFeature, 1e4);
 
634
        
 
635
        hGeom = OGR_F_GetGeometryRef(hFeature);
 
636
        if (hGeom == NULL) {
 
637
            G_debug(3, "Feature %d without geometry skipped", iFeature);
 
638
            OGR_F_Destroy(hFeature);
 
639
            nskipped++;
 
640
            continue;
 
641
        }
 
642
        
 
643
        FID = (int) OGR_F_GetFID(hFeature);
 
644
        if (FID == OGRNullFID) {
 
645
            G_debug(3, "OGR feature %d without ID skipped", iFeature);
 
646
            OGR_F_Destroy(hFeature);
 
647
            nskipped++;
 
648
            continue;
 
649
        }
 
650
        G_debug(4, "    FID = %d", FID);
 
651
        
 
652
        reset_parts(&parts);
 
653
        add_part(&parts, FID);
 
654
        npoints += add_geometry_ogr(&(Map->plus), ogr_info, hGeom,
 
655
                                    FID, build, &parts);
 
656
        
 
657
        OGR_F_Destroy(hFeature);
 
658
    } /* while */
 
659
    G_progress(1, 1);
 
660
    
 
661
    G_message(_n("One primitive registered", "%d primitives registered", Map->plus.n_lines), Map->plus.n_lines);
 
662
    G_message(_n("One vertex registered", "%d vertices registered", npoints), npoints);
 
663
    
 
664
    if (nskipped > 0)
 
665
        G_warning(_n("One feature without geometry skipped", "%d features without geometry skipped", nskipped), nskipped);
 
666
        
 
667
    Map->plus.built = GV_BUILD_BASE;
 
668
    
 
669
    free_parts(&parts);
 
670
}
 
671
#endif /* HAVE_OGR */
 
672
 
 
673
/*!
 
674
   \brief Build pseudo-topology (for simple features) - internal use only
 
675
 
 
676
   See Vect_build_ogr() and Vect_build_pg() for implemetation issues.
 
677
   
 
678
   Build levels:
 
679
    - GV_BUILD_NONE
 
680
    - GV_BUILD_BASE
 
681
    - GV_BUILD_ATTACH_ISLES
 
682
    - GV_BUILD_CENTROIDS
 
683
    - GV_BUILD_ALL
 
684
   
 
685
   \param Map pointer to Map_info structure
 
686
   \param build build level
 
687
 
 
688
   \return 1 on success
 
689
   \return 0 on error
 
690
*/
 
691
int Vect__build_sfa(struct Map_info *Map, int build)
 
692
{
 
693
    struct Plus_head *plus;
 
694
    
 
695
    plus = &(Map->plus);
 
696
    
 
697
    /* check if upgrade or downgrade */
 
698
    if (build < plus->built) {
 
699
        /* -> downgrade */
 
700
        Vect__build_downgrade(Map, build);
 
701
        return 1;
 
702
    }
 
703
    
 
704
    /* -> upgrade */
 
705
    if (plus->built < GV_BUILD_BASE) {
 
706
        if (Map->format == GV_FORMAT_OGR ||
 
707
            Map->format == GV_FORMAT_OGR_DIRECT) {
 
708
#ifdef HAVE_OGR
 
709
            build_ogr(Map, build);
 
710
#else
 
711
            G_fatal_error(_("GRASS is not compiled with OGR support"));
 
712
#endif
 
713
        }
 
714
        else if (Map->format == GV_FORMAT_POSTGIS) {
 
715
#ifdef HAVE_POSTGRES
 
716
            build_pg(Map, build);
 
717
#else
 
718
            G_fatal_error(_("GRASS is not compiled with PostgreSQL support"));
 
719
#endif
 
720
        }
 
721
        else {
 
722
            G_fatal_error(_("%s: Native format unsupported"),
 
723
                          "Vect__build_sfa()");
 
724
        }
 
725
    }
 
726
    
 
727
    plus->built = build;
 
728
    
 
729
    return 1;
 
730
}
 
731
 
 
732
/*!
 
733
   \brief Dump feature index to file
 
734
 
 
735
   \param Map pointer to Map_info struct
 
736
   \param out file for output (stdout/stderr for example)
 
737
 
 
738
   \return 1 on success
 
739
   \return 0 on error
 
740
 */
 
741
int Vect_fidx_dump(const struct Map_info *Map, FILE *out)
 
742
{
 
743
    int i;
 
744
    const struct Format_info_offset *offset;
 
745
 
 
746
    if (Map->format != GV_FORMAT_OGR &&
 
747
        Map->format != GV_FORMAT_POSTGIS) {
 
748
        G_warning(_("Feature index is built only for non-native formats. "
 
749
                    "Nothing to dump."));
 
750
        return 0;
 
751
    }
 
752
 
 
753
    if (Map->format == GV_FORMAT_OGR)
 
754
        offset = &(Map->fInfo.ogr.offset);
 
755
    else
 
756
        offset = &(Map->fInfo.pg.offset);
 
757
    
 
758
    fprintf(out, "---------- FEATURE INDEX DUMP ----------\n");
 
759
    
 
760
    fprintf(out, "format: %s\n", Vect_maptype_info(Map));
 
761
    if (Vect_maptype(Map) == GV_FORMAT_POSTGIS &&
 
762
        Map->fInfo.pg.toposchema_name)
 
763
        fprintf(out, "topology: PostGIS\n");
 
764
    else
 
765
        fprintf(out, "topology: pseudo\n");
 
766
    fprintf(out, "feature type: %s\n", 
 
767
            Vect_get_finfo_geometry_type(Map));
 
768
    fprintf(out, "number of features: %d\n\noffset : value (fid or part idx):\n",
 
769
            Vect_get_num_lines(Map));
 
770
    for (i = 0; i < offset->array_num; i++) {
 
771
        fprintf(out, "%6d : %d\n", i, offset->array[i]);
 
772
    }
 
773
 
 
774
    return 1;
 
775
}