← Back to branch summary

~ubuntu-branches/ubuntu/wily/grass/wily

« back to all changes in this revision

Viewing changes to vector/v.in.lidar/main.c

Tags: 7.0.0~rc1+ds1-1~exp1
http://bugs.debian.org/776349
* New upstream release candidate.
* Repack upstream tarball, remove precompiled Python objects.
* Add upstream metadata.
* Update gbp.conf and Vcs-Git URL to use the experimental branch.
* Update watch file for GRASS 7.0.
* Drop build dependencies for Tcl/Tk, add build dependencies:
  python-numpy, libnetcdf-dev, netcdf-bin, libblas-dev, liblapack-dev
* Update Vcs-Browser URL to use cgit instead of gitweb.
* Update paths to use grass70.
* Add configure options: --with-netcdf, --with-blas, --with-lapack,
  remove --with-tcltk-includes.
* Update patches for GRASS 7.
* Update copyright file, changes:
  - Update copyright years
  - Group files by license
  - Remove unused license sections
* Add patches for various typos.
* Fix desktop file with patch instead of d/rules.
* Use minimal dh rules.
* Bump Standards-Version to 3.9.6, no changes.
* Use dpkg-maintscript-helper to replace directories with symlinks.
  (closes: #776349)
* Update my email to use @debian.org address.

Show diffs side-by-side

added added

removed removed

Lines of Context:
vector/v.in.lidar/main.c
 
1
 
 
2
/****************************************************************
 
3
 *
 
4
 * MODULE:       v.in.lidar
 
5
 *
 
6
 * AUTHOR(S):    Markus Metz
 
7
 *               based on v.in.ogr
 
8
 *
 
9
 * PURPOSE:      Import LiDAR LAS points
 
10
 *
 
11
 * COPYRIGHT:    (C) 2011-2014 by the GRASS Development Team
 
12
 *
 
13
 *               This program is free software under the
 
14
 *               GNU General Public License (>=v2).
 
15
 *               Read the file COPYING that comes with GRASS
 
16
 *               for details.
 
17
 *
 
18
 * TODO: - make fixed field length of OFTIntegerList dynamic
 
19
 *       - several other TODOs below
 
20
**************************************************************/
 
21
 
 
22
#include <stdlib.h>
 
23
#include <string.h>
 
24
#include <unistd.h>
 
25
#include <grass/gis.h>
 
26
#include <grass/dbmi.h>
 
27
#include <grass/vector.h>
 
28
#include <grass/gprojects.h>
 
29
#include <grass/glocale.h>
 
30
#include <liblas/capi/liblas.h>
 
31
 
 
32
#ifndef MAX
 
33
#  define MIN(a,b)      ((a<b) ? a : b)
 
34
#  define MAX(a,b)      ((a>b) ? a : b)
 
35
#endif
 
36
 
 
37
#define LAS_ALL 0
 
38
#define LAS_FIRST 1
 
39
#define LAS_LAST 2
 
40
#define LAS_MID 3
 
41
 
 
42
/*
 
43
 * ASPRS Standard LIDAR Point Classes
 
44
 * Classification Value (bits 0:4) : Meaning
 
45
 *      0 : Created, never classified
 
46
 *      1 : Unclassified
 
47
 *      2 : Ground
 
48
 *      3 : Low Vegetation
 
49
 *      4 : Medium Vegetation
 
50
 *      5 : High Vegetation
 
51
 *      6 : Building
 
52
 *      7 : Low Point (noise)
 
53
 *      8 : Model Key-point (mass point)
 
54
 *      9 : Water
 
55
 *     10 : Reserved for ASPRS Definition
 
56
 *     11 : Reserved for ASPRS Definition
 
57
 *     12 : Overlap Points
 
58
 *  13-31 : Reserved for ASPRS Definition
 
59
 */
 
60
 
 
61
/* Classification Bit Field Encoding
 
62
 * Bits | Field Name     | Description
 
63
 *  0-4 | Classification | Standard ASPRS classification as defined in the
 
64
 *                         above classification table.
 
65
 *    5 | Synthetic      | If set then this point was created by a technique
 
66
 *                         other than LIDAR collection such as digitized from
 
67
 *                         a photogrammetric stereo model or by traversing
 
68
 *                         a waveform.
 
69
 *    6 | Key-point      | If set, this point is considered to be a model 
 
70
 *                         key-point and thus generally should not be withheld
 
71
 *                         in a thinning algorithm.
 
72
 *    7 | Withheld       | If set, this point should not be included in
 
73
 *                         processing (synonymous with Deleted).
 
74
*/
 
75
 
 
76
struct class_table
 
77
{
 
78
    int code;
 
79
    char *name;
 
80
};
 
81
 
 
82
static struct class_table class_val[] = {
 
83
    {0, "Created, never classified"},
 
84
    {1, "Unclassified"},
 
85
    {2, "Ground"},
 
86
    {3, "Low Vegetation"},
 
87
    {4, "Medium Vegetation"},
 
88
    {5, "High Vegetation"},
 
89
    {6, "Building"},
 
90
    {7, "Low Point (noise)"},
 
91
    {8, "Model Key-point (mass point)"},
 
92
    {9, "Water"},
 
93
    {10, "Reserved for ASPRS Definition"},
 
94
    {11, "Reserved for ASPRS Definition"},
 
95
    {12, "Overlap Points"},
 
96
    {13 /* 13 - 31 */, "Reserved for ASPRS Definition"},
 
97
    {0, 0}
 
98
};
 
99
 
 
100
static struct class_table class_type[] = {
 
101
    {5, "Synthetic"},
 
102
    {6, "Key-point"},
 
103
    {7, "Withheld"},
 
104
    {0, 0}
 
105
};
 
106
 
 
107
void print_lasinfo(LASHeaderH LAS_header, LASSRSH LAS_srs);
 
108
 
 
109
int main(int argc, char *argv[])
 
110
{
 
111
    int i;
 
112
    float xmin = 0., ymin = 0., xmax = 0., ymax = 0.;
 
113
    struct GModule *module;
 
114
    struct Option *in_opt, *out_opt, *spat_opt, *filter_opt, *class_opt;
 
115
    struct Option *outloc_opt;
 
116
    struct Flag *print_flag, *notab_flag, *region_flag, *notopo_flag;
 
117
    struct Flag *over_flag, *extend_flag, *no_import_flag;
 
118
    char buf[2000];
 
119
    struct Key_Value *loc_proj_info = NULL, *loc_proj_units = NULL;
 
120
    struct Key_Value *proj_info, *proj_units;
 
121
    const char *projstr;
 
122
    struct Cell_head cellhd, loc_wind, cur_wind;
 
123
    char error_msg[8192];
 
124
 
 
125
    /* Vector */
 
126
    struct Map_info Map;
 
127
    int cat;
 
128
 
 
129
    /* Attributes */
 
130
    struct field_info *Fi;
 
131
    dbDriver *driver;
 
132
    dbString sql, strval;
 
133
    
 
134
    /* LAS */
 
135
    LASReaderH LAS_reader;
 
136
    LASHeaderH LAS_header;
 
137
    LASSRSH LAS_srs;
 
138
    LASPointH LAS_point;
 
139
    double scale_x, scale_y, scale_z, offset_x, offset_y, offset_z;
 
140
    int las_point_format;
 
141
    int have_time, have_color;
 
142
    int return_filter;
 
143
    int skipme;
 
144
    int point_class;
 
145
    unsigned int not_valid;     
 
146
 
 
147
    struct line_pnts *Points;
 
148
    struct line_cats *Cats;
 
149
 
 
150
    unsigned int n_features, feature_count, n_outside, n_filtered, n_class_filtered;
 
151
    int overwrite;
 
152
 
 
153
    G_gisinit(argv[0]);
 
154
 
 
155
    module = G_define_module();
 
156
    G_add_keyword(_("vector"));
 
157
    G_add_keyword(_("import"));
 
158
    G_add_keyword(_("LIDAR"));
 
159
    module->description = _("Converts LAS LiDAR point clouds to a GRASS vector map with libLAS.");
 
160
 
 
161
    in_opt = G_define_standard_option(G_OPT_F_INPUT);
 
162
    in_opt->label = _("LAS input file");
 
163
    in_opt->description = _("LiDAR input files in LAS format (*.las or *.laz)");
 
164
 
 
165
    out_opt = G_define_standard_option(G_OPT_V_OUTPUT);
 
166
    
 
167
    spat_opt = G_define_option();
 
168
    spat_opt->key = "spatial";
 
169
    spat_opt->type = TYPE_DOUBLE;
 
170
    spat_opt->multiple = YES;
 
171
    spat_opt->required = NO;
 
172
    spat_opt->key_desc = "xmin,ymin,xmax,ymax";
 
173
    spat_opt->label = _("Import subregion only");
 
174
    spat_opt->guisection = _("Subregion");
 
175
    spat_opt->description =
 
176
        _("Format: xmin,ymin,xmax,ymax - usually W,S,E,N");
 
177
 
 
178
    outloc_opt = G_define_option();
 
179
    outloc_opt->key = "location";
 
180
    outloc_opt->type = TYPE_STRING;
 
181
    outloc_opt->required = NO;
 
182
    outloc_opt->description = _("Name for new location to create");
 
183
    outloc_opt->key_desc = "name";
 
184
    
 
185
    filter_opt = G_define_option();
 
186
    filter_opt->key = "return_filter";
 
187
    filter_opt->type = TYPE_STRING;
 
188
    filter_opt->required = NO;
 
189
    filter_opt->label = _("Only import points of selected return type");
 
190
    filter_opt->description = _("If not specified, all points are imported");
 
191
    filter_opt->options = "first,last,mid";
 
192
 
 
193
    class_opt = G_define_option();
 
194
    class_opt->key = "class_filter";
 
195
    class_opt->type = TYPE_INTEGER;
 
196
    class_opt->multiple = YES;
 
197
    class_opt->required = NO;
 
198
    class_opt->label = _("Only import points of selected class(es)");
 
199
    class_opt->description = _("Input is comma separated integers. "
 
200
                               "If not specified, all points are imported.");
 
201
 
 
202
    print_flag = G_define_flag();
 
203
    print_flag->key = 'p';
 
204
    print_flag->description =
 
205
        _("Print LAS file info and exit");
 
206
    print_flag->suppress_required = YES;
 
207
    
 
208
    notab_flag = G_define_standard_flag(G_FLG_V_TABLE);
 
209
    notab_flag->guisection = _("Attributes");
 
210
 
 
211
    over_flag = G_define_flag();
 
212
    over_flag->key = 'o';
 
213
    over_flag->description =
 
214
        _("Override dataset projection (use location's projection)");
 
215
 
 
216
    region_flag = G_define_flag();
 
217
    region_flag->key = 'r';
 
218
    region_flag->guisection = _("Subregion");
 
219
    region_flag->description = _("Limit import to the current region");
 
220
 
 
221
    extend_flag = G_define_flag();
 
222
    extend_flag->key = 'e';
 
223
    extend_flag->description =
 
224
        _("Extend region extents based on new dataset");
 
225
 
 
226
    no_import_flag = G_define_flag();
 
227
    no_import_flag->key = 'i';
 
228
    no_import_flag->description =
 
229
        _("Create the location specified by the \"location\" parameter and exit."
 
230
          " Do not import the vector data.");
 
231
    no_import_flag->suppress_required = YES;
 
232
 
 
233
    notopo_flag = G_define_standard_flag(G_FLG_V_TOPO);
 
234
 
 
235
    /* The parser checks if the map already exists in current mapset, this is
 
236
     * wrong if location options is used, so we switch out the check and do it
 
237
     * in the module after the parser */
 
238
    overwrite = G_check_overwrite(argc, argv);
 
239
 
 
240
    if (G_parser(argc, argv))
 
241
        exit(EXIT_FAILURE);
 
242
 
 
243
    /* Don't crash on cmd line if file not found */
 
244
    if (access(in_opt->answer, F_OK) != 0) {
 
245
        G_fatal_error(_("Input file <%s> does not exist"), in_opt->answer);
 
246
    }
 
247
    /* Open LAS file*/
 
248
    LAS_reader = LASReader_Create(in_opt->answer);
 
249
    LAS_header = LASReader_GetHeader(LAS_reader);
 
250
 
 
251
    if  (LAS_header == NULL) {
 
252
        G_fatal_error(_("Input file <%s> is not a LAS LiDAR point cloud"),
 
253
                        in_opt->answer);
 
254
    }
 
255
 
 
256
    LAS_srs = LASHeader_GetSRS(LAS_header);
 
257
 
 
258
    scale_x = LASHeader_GetScaleX(LAS_header);
 
259
    scale_y = LASHeader_GetScaleY(LAS_header);
 
260
    scale_z = LASHeader_GetScaleZ(LAS_header);
 
261
 
 
262
    offset_x = LASHeader_GetOffsetX(LAS_header);
 
263
    offset_y = LASHeader_GetOffsetY(LAS_header);
 
264
    offset_z = LASHeader_GetOffsetZ(LAS_header);
 
265
 
 
266
    xmin = LASHeader_GetMinX(LAS_header);
 
267
    xmax = LASHeader_GetMaxX(LAS_header);
 
268
    ymin = LASHeader_GetMinY(LAS_header);
 
269
    ymax = LASHeader_GetMaxY(LAS_header);
 
270
 
 
271
    /* Print LAS header */
 
272
    if (print_flag->answer) {
 
273
        /* print... */
 
274
        print_lasinfo(LAS_header, LAS_srs);
 
275
 
 
276
        LASSRS_Destroy(LAS_srs);
 
277
        LASHeader_Destroy(LAS_header);
 
278
        LASReader_Destroy(LAS_reader);
 
279
 
 
280
        exit(EXIT_SUCCESS);
 
281
    }
 
282
 
 
283
    return_filter = LAS_ALL;
 
284
    if (filter_opt->answer) {
 
285
        if (strcmp(filter_opt->answer, "first") == 0)
 
286
            return_filter = LAS_FIRST;
 
287
        else if (strcmp(filter_opt->answer, "last") == 0)
 
288
            return_filter = LAS_LAST;
 
289
        else if (strcmp(filter_opt->answer, "mid") == 0)
 
290
            return_filter = LAS_MID;
 
291
        else
 
292
            G_fatal_error(_("Unknown filter option <%s>"), filter_opt->answer);
 
293
    }
 
294
 
 
295
    if (region_flag->answer) {
 
296
        if (spat_opt->answer)
 
297
            G_fatal_error(_("Select either the current region flag or the spatial option, not both"));
 
298
 
 
299
        G_get_window(&cur_wind);
 
300
        xmin = cur_wind.west;
 
301
        xmax = cur_wind.east;
 
302
        ymin = cur_wind.south;
 
303
        ymax = cur_wind.north;
 
304
    }
 
305
    if (spat_opt->answer) {
 
306
        /* See as reference: gdal/ogr/ogr_capi_test.c */
 
307
 
 
308
        /* cut out a piece of the map */
 
309
        /* order: xmin,ymin,xmax,ymax */
 
310
        int arg_s_num = 0;
 
311
        i = 0;
 
312
        while (spat_opt->answers[i]) {
 
313
            if (i == 0)
 
314
                xmin = atof(spat_opt->answers[i]);
 
315
            if (i == 1)
 
316
                ymin = atof(spat_opt->answers[i]);
 
317
            if (i == 2)
 
318
                xmax = atof(spat_opt->answers[i]);
 
319
            if (i == 3)
 
320
                ymax = atof(spat_opt->answers[i]);
 
321
            arg_s_num++;
 
322
            i++;
 
323
        }
 
324
        if (arg_s_num != 4)
 
325
            G_fatal_error(_("4 parameters required for 'spatial' parameter"));
 
326
    }
 
327
    if (spat_opt->answer || region_flag->answer) {
 
328
        G_debug(2, "cut out with boundaries: xmin:%f ymin:%f xmax:%f ymax:%f",
 
329
                xmin, ymin, xmax, ymax);
 
330
    }
 
331
 
 
332
    /* fetch boundaries */
 
333
    G_get_window(&cellhd);
 
334
    cellhd.north = ymax;
 
335
    cellhd.south = ymin;
 
336
    cellhd.west = xmin;
 
337
    cellhd.east = xmax;
 
338
    cellhd.rows = 20;   /* TODO - calculate useful values */
 
339
    cellhd.cols = 20;
 
340
    cellhd.ns_res = (cellhd.north - cellhd.south) / cellhd.rows;
 
341
    cellhd.ew_res = (cellhd.east - cellhd.west) / cellhd.cols;
 
342
 
 
343
    /* Fetch input map projection in GRASS form. */
 
344
    proj_info = NULL;
 
345
    proj_units = NULL;
 
346
    projstr = LASSRS_GetWKT_CompoundOK(LAS_srs);
 
347
 
 
348
    /* Do we need to create a new location? */
 
349
    if (outloc_opt->answer != NULL) {
 
350
        /* Convert projection information non-interactively as we can't
 
351
         * assume the user has a terminal open */
 
352
        if (GPJ_wkt_to_grass(&cellhd, &proj_info,
 
353
                             &proj_units, projstr, 0) < 0) {
 
354
            G_fatal_error(_("Unable to convert input map projection to GRASS "
 
355
                            "format; cannot create new location."));
 
356
        }
 
357
        else {
 
358
            if (0 != G_make_location(outloc_opt->answer, &cellhd,
 
359
                                     proj_info, proj_units)) {
 
360
                G_fatal_error(_("Unable to create new location <%s>"),
 
361
                              outloc_opt->answer);
 
362
            }
 
363
            G_message(_("Location <%s> created"), outloc_opt->answer);
 
364
        }
 
365
 
 
366
        /* If the i flag is set, clean up? and exit here */
 
367
        if(no_import_flag->answer)
 
368
            exit(EXIT_SUCCESS);
 
369
 
 
370
        /*  TODO: */
 
371
        G_warning("Import into new location not yet implemented");
 
372
        /* at this point the module should be using G_create_alt_env()
 
373
            to change context to the newly created location; once done
 
374
            it should switch back with G_switch_env(). See r.in.gdal */
 
375
    }
 
376
    else {
 
377
        int err = 0;
 
378
 
 
379
        /* Projection only required for checking so convert non-interactively */
 
380
        if (GPJ_wkt_to_grass(&cellhd, &proj_info,
 
381
                             &proj_units, projstr, 0) < 0)
 
382
            G_warning(_("Unable to convert input map projection information to "
 
383
                       "GRASS format for checking"));
 
384
 
 
385
        /* Does the projection of the current location match the dataset? */
 
386
        /* G_get_window seems to be unreliable if the location has been changed */
 
387
        G_get_default_window(&loc_wind);
 
388
        /* fetch LOCATION PROJ info */
 
389
        if (loc_wind.proj != PROJECTION_XY) {
 
390
            loc_proj_info = G_get_projinfo();
 
391
            loc_proj_units = G_get_projunits();
 
392
        }
 
393
 
 
394
        if (over_flag->answer) {
 
395
            cellhd.proj = loc_wind.proj;
 
396
            cellhd.zone = loc_wind.zone;
 
397
            G_message(_("Over-riding projection check"));
 
398
        }
 
399
        else if (loc_wind.proj != cellhd.proj
 
400
                 || (err =
 
401
                     G_compare_projections(loc_proj_info, loc_proj_units,
 
402
                                           proj_info, proj_units)) != TRUE) {
 
403
            int i_value;
 
404
 
 
405
            strcpy(error_msg,
 
406
                   _("Projection of dataset does not"
 
407
                     " appear to match current location.\n\n"));
 
408
 
 
409
            /* TODO: output this info sorted by key: */
 
410
            if (loc_wind.proj != cellhd.proj || err != -2) {
 
411
                if (loc_proj_info != NULL) {
 
412
                    strcat(error_msg, _("GRASS LOCATION PROJ_INFO is:\n"));
 
413
                    for (i_value = 0; i_value < loc_proj_info->nitems;
 
414
                         i_value++)
 
415
                        sprintf(error_msg + strlen(error_msg), "%s: %s\n",
 
416
                                loc_proj_info->key[i_value],
 
417
                                loc_proj_info->value[i_value]);
 
418
                    strcat(error_msg, "\n");
 
419
                }
 
420
 
 
421
                if (proj_info != NULL) {
 
422
                    strcat(error_msg, _("Import dataset PROJ_INFO is:\n"));
 
423
                    for (i_value = 0; i_value < proj_info->nitems; i_value++)
 
424
                        sprintf(error_msg + strlen(error_msg), "%s: %s\n",
 
425
                                proj_info->key[i_value],
 
426
                                proj_info->value[i_value]);
 
427
                }
 
428
                else {
 
429
                    strcat(error_msg, _("Import dataset PROJ_INFO is:\n"));
 
430
                    if (cellhd.proj == PROJECTION_XY)
 
431
                        sprintf(error_msg + strlen(error_msg),
 
432
                                "Dataset proj = %d (unreferenced/unknown)\n",
 
433
                                cellhd.proj);
 
434
                    else if (cellhd.proj == PROJECTION_LL)
 
435
                        sprintf(error_msg + strlen(error_msg),
 
436
                                "Dataset proj = %d (lat/long)\n",
 
437
                                cellhd.proj);
 
438
                    else if (cellhd.proj == PROJECTION_UTM)
 
439
                        sprintf(error_msg + strlen(error_msg),
 
440
                                "Dataset proj = %d (UTM), zone = %d\n",
 
441
                                cellhd.proj, cellhd.zone);
 
442
                    else
 
443
                        sprintf(error_msg + strlen(error_msg),
 
444
                                "Dataset proj = %d (unknown), zone = %d\n",
 
445
                                cellhd.proj, cellhd.zone);
 
446
                }
 
447
            }
 
448
            else {
 
449
                if (loc_proj_units != NULL) {
 
450
                    strcat(error_msg, "GRASS LOCATION PROJ_UNITS is:\n");
 
451
                    for (i_value = 0; i_value < loc_proj_units->nitems;
 
452
                         i_value++)
 
453
                        sprintf(error_msg + strlen(error_msg), "%s: %s\n",
 
454
                                loc_proj_units->key[i_value],
 
455
                                loc_proj_units->value[i_value]);
 
456
                    strcat(error_msg, "\n");
 
457
                }
 
458
 
 
459
                if (proj_units != NULL) {
 
460
                    strcat(error_msg, "Import dataset PROJ_UNITS is:\n");
 
461
                    for (i_value = 0; i_value < proj_units->nitems; i_value++)
 
462
                        sprintf(error_msg + strlen(error_msg), "%s: %s\n",
 
463
                                proj_units->key[i_value],
 
464
                                proj_units->value[i_value]);
 
465
                }
 
466
            }
 
467
            sprintf(error_msg + strlen(error_msg),
 
468
                    _("\nYou can use the -o flag to %s to override this projection check.\n"),
 
469
                    G_program_name());
 
470
            strcat(error_msg,
 
471
                   _("Consider generating a new location with 'location' parameter"
 
472
                    " from input data set.\n"));
 
473
            G_fatal_error(error_msg);
 
474
        }
 
475
        else {
 
476
            G_verbose_message(_("Projection of input dataset and current "
 
477
                                "location appear to match"));
 
478
        }
 
479
    }
 
480
 
 
481
    db_init_string(&sql);
 
482
    db_init_string(&strval);
 
483
 
 
484
    if (!outloc_opt->answer) {  /* Check if the map exists */
 
485
        if (G_find_vector2(out_opt->answer, G_mapset())) {
 
486
            if (overwrite)
 
487
                G_warning(_("Vector map <%s> already exists and will be overwritten"),
 
488
                          out_opt->answer);
 
489
            else
 
490
                G_fatal_error(_("Vector map <%s> already exists"),
 
491
                              out_opt->answer);
 
492
        }
 
493
    }
 
494
 
 
495
    /* open output vector */
 
496
    sprintf(buf, "%s", out_opt->answer);
 
497
    /* strip any @mapset from vector output name */
 
498
    G_find_vector(buf, G_mapset());
 
499
    if (Vect_open_new(&Map, out_opt->answer, 1) < 0)
 
500
        G_fatal_error(_("Unable to create vector map <%s>"), out_opt->answer);
 
501
 
 
502
    Vect_hist_command(&Map);
 
503
    
 
504
    n_features = LASHeader_GetPointRecordsCount(LAS_header);
 
505
    las_point_format = LASHeader_GetDataFormatId(LAS_header);
 
506
 
 
507
    have_time = (las_point_format == 1 || las_point_format == 3 || 
 
508
                 las_point_format == 4 || las_point_format == 5);
 
509
 
 
510
    have_color = (las_point_format == 2 || las_point_format == 3 || 
 
511
                   las_point_format == 5);
 
512
 
 
513
    /* Add DB link */
 
514
    if (!notab_flag->answer) {
 
515
        char *cat_col_name = GV_KEY_COLUMN;
 
516
 
 
517
        Fi = Vect_default_field_info(&Map, 1, NULL, GV_1TABLE);
 
518
 
 
519
        Vect_map_add_dblink(&Map, 1, out_opt->answer, Fi->table,
 
520
                            cat_col_name, Fi->database, Fi->driver);
 
521
 
 
522
        /* check available LAS info, depends on POINT DATA RECORD FORMAT [0-5] */
 
523
        /* X (double),
 
524
         * Y (double), 
 
525
         * Z (double), 
 
526
         * intensity (double), 
 
527
         * return number (int), 
 
528
         * number of returns (int),
 
529
         * scan direction (int),
 
530
         * flight line edge (int),
 
531
         * classification type (char),
 
532
         * class (char),
 
533
         * time (double) (FORMAT 1, 3, 4, 5),
 
534
         * scan angle rank (int),
 
535
         * source ID (int),
 
536
         * user data (char), ???
 
537
         * red (int)  (FORMAT 2, 3, 5),
 
538
         * green (int) (FORMAT 2, 3, 5),
 
539
         * blue (int) (FORMAT 2, 3, 5)*/
 
540
         
 
541
        /* Create table */
 
542
        sprintf(buf, "create table %s (%s integer", Fi->table,
 
543
                cat_col_name);
 
544
        db_set_string(&sql, buf);
 
545
        
 
546
        /* x, y, z */
 
547
        sprintf(buf, ", x_coord double precision");
 
548
        db_append_string(&sql, buf);
 
549
        sprintf(buf, ", y_coord double precision");
 
550
        db_append_string(&sql, buf);
 
551
        sprintf(buf, ", z_coord double precision");
 
552
        db_append_string(&sql, buf);
 
553
        /* intensity */
 
554
        sprintf(buf, ", intensity integer");
 
555
        db_append_string(&sql, buf);
 
556
        /* return number */
 
557
        sprintf(buf, ", return integer");
 
558
        db_append_string(&sql, buf);
 
559
        /* number of returns */
 
560
        sprintf(buf, ", n_returns integer");
 
561
        db_append_string(&sql, buf);
 
562
        /* scan direction */
 
563
        sprintf(buf, ", scan_dir integer");
 
564
        db_append_string(&sql, buf);
 
565
        /* flight line edge */
 
566
        sprintf(buf, ", edge integer");
 
567
        db_append_string(&sql, buf);
 
568
        /* classification type */
 
569
        sprintf(buf, ", cl_type varchar(20)");
 
570
        db_append_string(&sql, buf);
 
571
        /* classification class */
 
572
        sprintf(buf, ", class varchar(40)");
 
573
        db_append_string(&sql, buf);
 
574
        /* GPS time */
 
575
        if (have_time) {
 
576
            sprintf(buf, ", gps_time double precision");
 
577
            db_append_string(&sql, buf);
 
578
        }
 
579
        /* scan angle */
 
580
        sprintf(buf, ", angle integer");
 
581
        db_append_string(&sql, buf);
 
582
        /* source id */
 
583
        sprintf(buf, ", src_id integer");
 
584
        db_append_string(&sql, buf);
 
585
        /* user data */
 
586
        sprintf(buf, ", usr_data integer");
 
587
        db_append_string(&sql, buf);
 
588
        /* colors */
 
589
        if (have_color) {
 
590
            sprintf(buf, ", red integer, green integer, blue integer");
 
591
            db_append_string(&sql, buf);
 
592
            sprintf(buf, ", GRASSRGB varchar(11)");
 
593
            db_append_string(&sql, buf);
 
594
        }
 
595
 
 
596
        db_append_string(&sql, ")");
 
597
        G_debug(3, db_get_string(&sql));
 
598
 
 
599
        driver =
 
600
            db_start_driver_open_database(Fi->driver,
 
601
                                          Vect_subst_var(Fi->database,
 
602
                                                         &Map));
 
603
        if (driver == NULL) {
 
604
            G_fatal_error(_("Unable open database <%s> by driver <%s>"),
 
605
                          Vect_subst_var(Fi->database, &Map), Fi->driver);
 
606
        }
 
607
        db_set_error_handler_driver(driver);
 
608
 
 
609
        if (db_execute_immediate(driver, &sql) != DB_OK) {
 
610
            G_fatal_error(_("Unable to create table: '%s'"),
 
611
                          db_get_string(&sql));
 
612
        }
 
613
 
 
614
        if (db_create_index2(driver, Fi->table, cat_col_name) != DB_OK)
 
615
            G_warning(_("Unable to create index for table <%s>, key <%s>"),
 
616
                      Fi->table, cat_col_name);
 
617
 
 
618
        if (db_grant_on_table
 
619
            (driver, Fi->table, DB_PRIV_SELECT,
 
620
             DB_GROUP | DB_PUBLIC) != DB_OK)
 
621
            G_fatal_error(_("Unable to grant privileges on table <%s>"),
 
622
                          Fi->table);
 
623
 
 
624
        db_begin_transaction(driver);
 
625
    }
 
626
 
 
627
    /* Import feature */
 
628
    cat = 1;
 
629
    not_valid = 0;
 
630
    feature_count = 0;
 
631
    n_outside = 0;
 
632
    n_filtered = 0;
 
633
    n_class_filtered = 0;
 
634
 
 
635
    Points = Vect_new_line_struct();
 
636
    Cats = Vect_new_cats_struct();
 
637
    
 
638
    G_important_message(_("Scanning %d points..."), n_features);
 
639
    while ((LAS_point = LASReader_GetNextPoint(LAS_reader)) != NULL) {
 
640
        double x, y, z;
 
641
 
 
642
        G_percent(feature_count++, n_features, 1);      /* show something happens */
 
643
        
 
644
        if (!LASPoint_IsValid(LAS_point)) {
 
645
            not_valid++;
 
646
            continue;
 
647
        }
 
648
 
 
649
        Vect_reset_line(Points);
 
650
        Vect_reset_cats(Cats);
 
651
 
 
652
        x = LASPoint_GetX(LAS_point);
 
653
        y = LASPoint_GetY(LAS_point);
 
654
        z = LASPoint_GetZ(LAS_point);
 
655
 
 
656
        if (spat_opt->answer || region_flag->answer) {
 
657
            if (x < xmin || x > xmax || y < ymin || y > ymax) {
 
658
                n_outside++;
 
659
                continue;
 
660
            }
 
661
        }
 
662
        if (return_filter != LAS_ALL) {
 
663
            int return_no = LASPoint_GetReturnNumber(LAS_point);
 
664
            int n_returns = LASPoint_GetNumberOfReturns(LAS_point);
 
665
            skipme = 1;
 
666
 
 
667
            switch (return_filter) {
 
668
            case LAS_FIRST:
 
669
                if (return_no == 1)
 
670
                    skipme = 0;
 
671
                break;
 
672
            case LAS_MID:
 
673
                if (return_no > 1 && return_no < n_returns)
 
674
                    skipme = 0;
 
675
                break;
 
676
            case LAS_LAST:
 
677
                if (n_returns > 1 && return_no == n_returns)
 
678
                    skipme = 0;
 
679
                break;
 
680
            }
 
681
            
 
682
            if (skipme) {
 
683
                n_filtered++;
 
684
                continue;
 
685
            }
 
686
        }
 
687
        if (class_opt->answer) {
 
688
            point_class = (int) LASPoint_GetClassification(LAS_point);
 
689
            i = 0;
 
690
            skipme = TRUE;
 
691
            while (class_opt->answers[i]) {
 
692
                if (point_class == atoi(class_opt->answers[i])) {
 
693
                    skipme = FALSE;
 
694
                    break;
 
695
                }
 
696
                i++;
 
697
            }
 
698
            if (skipme) {
 
699
                n_class_filtered++;
 
700
                continue;
 
701
            }
 
702
        }
 
703
 
 
704
        Vect_append_point(Points, x, y, z);
 
705
        Vect_cat_set(Cats, 1, cat);
 
706
        Vect_write_line(&Map, GV_POINT, Points, Cats);
 
707
 
 
708
        /* Attributes */
 
709
        if (!notab_flag->answer) {
 
710
            char class_flag;
 
711
            int las_class_type, las_class;
 
712
 
 
713
             /* use LASPoint_Validate (LASPointH hPoint) to check for
 
714
              * return number, number of returns, scan direction, flight line edge,
 
715
              * classification, scan angle rank */
 
716
            sprintf(buf, "insert into %s values ( %d", Fi->table, cat);
 
717
            db_set_string(&sql, buf);
 
718
 
 
719
            /* x, y, z */
 
720
            sprintf(buf, ", %f", x);
 
721
            db_append_string(&sql, buf);
 
722
            sprintf(buf, ", %f", y);
 
723
            db_append_string(&sql, buf);
 
724
            sprintf(buf, ", %f", z);
 
725
            db_append_string(&sql, buf);
 
726
            /* intensity */
 
727
            sprintf(buf, ", %d", LASPoint_GetIntensity(LAS_point));
 
728
            db_append_string(&sql, buf);
 
729
            /* return number */
 
730
            sprintf(buf, ", %d", LASPoint_GetReturnNumber(LAS_point));
 
731
            db_append_string(&sql, buf);
 
732
            /* number of returns */
 
733
            sprintf(buf, ",  %d", LASPoint_GetNumberOfReturns(LAS_point));
 
734
            db_append_string(&sql, buf);
 
735
            /* scan direction */
 
736
            sprintf(buf, ", %d",  LASPoint_GetScanDirection(LAS_point));
 
737
            db_append_string(&sql, buf);
 
738
            /* flight line edge */
 
739
            sprintf(buf, ",  %d", LASPoint_GetFlightLineEdge(LAS_point));
 
740
            db_append_string(&sql, buf);
 
741
            class_flag = LASPoint_GetClassification(LAS_point);
 
742
            /* classification type int or char ? */
 
743
            las_class_type = class_flag / 32;
 
744
            sprintf(buf, ", \"%s\"", class_type[las_class_type].name);
 
745
            db_append_string(&sql, buf);
 
746
            /* classification class int or char ? */
 
747
            las_class = class_flag % 32;
 
748
            if (las_class > 13)
 
749
                las_class = 13;
 
750
            sprintf(buf, ", \"%s\"", class_val[las_class].name);
 
751
            db_append_string(&sql, buf);
 
752
            /* GPS time */
 
753
            if (have_time) {
 
754
                sprintf(buf, ", %f", LASPoint_GetTime(LAS_point));
 
755
                db_append_string(&sql, buf);
 
756
            }
 
757
            /* scan angle */
 
758
            sprintf(buf, ", %d", LASPoint_GetScanAngleRank(LAS_point));
 
759
            db_append_string(&sql, buf);
 
760
            /* source id */
 
761
            sprintf(buf, ", %d", LASPoint_GetPointSourceId(LAS_point));
 
762
            db_append_string(&sql, buf);
 
763
            /* user data */
 
764
            sprintf(buf, ", %d", LASPoint_GetUserData(LAS_point));
 
765
            db_append_string(&sql, buf);
 
766
            /* colors */
 
767
            if (have_color) {
 
768
                LASColorH LAS_color = LASPoint_GetColor(LAS_point);
 
769
                int red = LASColor_GetRed(LAS_color);
 
770
                int green = LASColor_GetGreen(LAS_color);
 
771
                int blue = LASColor_GetBlue(LAS_color);
 
772
 
 
773
                sprintf(buf, ", %d, %d, %d", red, green, blue);
 
774
                db_append_string(&sql, buf);
 
775
                sprintf(buf, ", \"%03d:%03d:%03d\"", red, green, blue);
 
776
                db_append_string(&sql, buf);
 
777
            }
 
778
            db_append_string(&sql, " )");
 
779
            G_debug(3, db_get_string(&sql));
 
780
 
 
781
            if (db_execute_immediate(driver, &sql) != DB_OK) {
 
782
                G_fatal_error(_("Cannot insert new row: %s"),
 
783
                              db_get_string(&sql));
 
784
            }
 
785
        }
 
786
 
 
787
        cat++;
 
788
    }
 
789
    G_percent(n_features, n_features, 1);       /* finish it */
 
790
 
 
791
    if (!notab_flag->answer) {
 
792
        db_commit_transaction(driver);
 
793
        db_close_database_shutdown_driver(driver);
 
794
    }
 
795
    
 
796
    LASSRS_Destroy(LAS_srs);
 
797
    LASHeader_Destroy(LAS_header);
 
798
    LASReader_Destroy(LAS_reader);
 
799
 
 
800
    /* close map */
 
801
    if (!notopo_flag->answer)
 
802
        Vect_build(&Map);
 
803
    Vect_close(&Map);
 
804
    
 
805
    G_message(_("%d points imported"),
 
806
              n_features - not_valid - n_outside - n_filtered - n_class_filtered);
 
807
    if (not_valid)
 
808
        G_message(_("%d input points were not valid"), not_valid);
 
809
    if (n_outside)
 
810
        G_message(_("%d input points were outside of the selected area"), n_outside);
 
811
    if (n_filtered)
 
812
        G_message(_("%d input points were filtered out by return number"), n_filtered);
 
813
    if (n_class_filtered)
 
814
        G_message(_("%d input points were filtered out by class number"), n_class_filtered);
 
815
 
 
816
    /* -------------------------------------------------------------------- */
 
817
    /*      Extend current window based on dataset.                         */
 
818
    /* -------------------------------------------------------------------- */
 
819
    if (extend_flag->answer) {
 
820
        G_get_set_window(&loc_wind);
 
821
 
 
822
        loc_wind.north = MAX(loc_wind.north, cellhd.north);
 
823
        loc_wind.south = MIN(loc_wind.south, cellhd.south);
 
824
        loc_wind.west = MIN(loc_wind.west, cellhd.west);
 
825
        loc_wind.east = MAX(loc_wind.east, cellhd.east);
 
826
 
 
827
        loc_wind.rows = (int)ceil((loc_wind.north - loc_wind.south)
 
828
                                  / loc_wind.ns_res);
 
829
        loc_wind.south = loc_wind.north - loc_wind.rows * loc_wind.ns_res;
 
830
 
 
831
        loc_wind.cols = (int)ceil((loc_wind.east - loc_wind.west)
 
832
                                  / loc_wind.ew_res);
 
833
        loc_wind.east = loc_wind.west + loc_wind.cols * loc_wind.ew_res;
 
834
 
 
835
        G_put_window(&loc_wind);
 
836
    }
 
837
 
 
838
    exit(EXIT_SUCCESS);
 
839
}
 
840
 
 
841
void print_lasinfo(LASHeaderH LAS_header, LASSRSH LAS_srs)
 
842
{
 
843
    char *las_srs_proj4 = LASSRS_GetProj4(LAS_srs);
 
844
    int las_point_format = LASHeader_GetDataFormatId(LAS_header);
 
845
 
 
846
    fprintf(stdout, "\nUsing LAS Library Version '%s'\n\n",
 
847
                    LAS_GetFullVersion());
 
848
    fprintf(stdout, "LAS File Version:                  %d.%d\n",
 
849
                    LASHeader_GetVersionMajor(LAS_header),
 
850
                    LASHeader_GetVersionMinor(LAS_header));
 
851
    fprintf(stdout, "System ID:                         '%s'\n",
 
852
                    LASHeader_GetSystemId(LAS_header));
 
853
    fprintf(stdout, "Generating Software:               '%s'\n",
 
854
                    LASHeader_GetSoftwareId(LAS_header));
 
855
    fprintf(stdout, "File Creation Day/Year:            %d/%d\n",
 
856
                    LASHeader_GetCreationDOY(LAS_header),
 
857
                    LASHeader_GetCreationYear(LAS_header));
 
858
    fprintf(stdout, "Point Data Format:                 %d\n",
 
859
                    las_point_format);
 
860
    fprintf(stdout, "Number of Point Records:           %d\n",
 
861
                    LASHeader_GetPointRecordsCount(LAS_header));
 
862
    fprintf(stdout, "Number of Points by Return:        %d %d %d %d %d\n",
 
863
                    LASHeader_GetPointRecordsByReturnCount(LAS_header, 0),
 
864
                    LASHeader_GetPointRecordsByReturnCount(LAS_header, 1),
 
865
                    LASHeader_GetPointRecordsByReturnCount(LAS_header, 2),
 
866
                    LASHeader_GetPointRecordsByReturnCount(LAS_header, 3),
 
867
                    LASHeader_GetPointRecordsByReturnCount(LAS_header, 4));
 
868
    fprintf(stdout, "Scale Factor X Y Z:                %g %g %g\n",
 
869
                    LASHeader_GetScaleX(LAS_header),
 
870
                    LASHeader_GetScaleY(LAS_header),
 
871
                    LASHeader_GetScaleZ(LAS_header));
 
872
    fprintf(stdout, "Offset X Y Z:                      %g %g %g\n",
 
873
                    LASHeader_GetOffsetX(LAS_header),
 
874
                    LASHeader_GetOffsetY(LAS_header),
 
875
                    LASHeader_GetOffsetZ(LAS_header));
 
876
    fprintf(stdout, "Min X Y Z:                         %g %g %g\n",
 
877
                    LASHeader_GetMinX(LAS_header),
 
878
                    LASHeader_GetMinY(LAS_header),
 
879
                    LASHeader_GetMinZ(LAS_header));
 
880
    fprintf(stdout, "Max X Y Z:                         %g %g %g\n",
 
881
                    LASHeader_GetMaxX(LAS_header),
 
882
                    LASHeader_GetMaxY(LAS_header),
 
883
                    LASHeader_GetMaxZ(LAS_header));
 
884
    if (las_srs_proj4 && strlen(las_srs_proj4) > 0) {
 
885
        fprintf(stdout, "Spatial Reference:\n");
 
886
        fprintf(stdout, "%s\n", las_srs_proj4);
 
887
    }
 
888
    else {
 
889
        fprintf(stdout, "Spatial Reference:                 None\n");
 
890
    }
 
891
    
 
892
    fprintf(stdout, "\nData Fields:\n");
 
893
    fprintf(stdout, "  'X'\n  'Y'\n  'Z'\n  'Intensity'\n  'Return Number'\n");
 
894
    fprintf(stdout, "  'Number of Returns'\n  'Scan Direction'\n");
 
895
    fprintf(stdout, "  'Flighline Edge'\n  'Classification'\n  'Scan Angle Rank'\n");
 
896
    fprintf(stdout, "  'User Data'\n  'Point Source ID'\n");
 
897
    if (las_point_format == 1 || las_point_format == 3 || las_point_format == 4 || las_point_format == 5) {
 
898
        fprintf(stdout, "  'GPS Time'\n");
 
899
    }
 
900
    if (las_point_format == 2 || las_point_format == 3 || las_point_format == 5) {
 
901
        fprintf(stdout, "  'Red'\n  'Green'\n  'Blue'\n");
 
902
    }
 
903
    fprintf(stdout, "\n");
 
904
    fflush(stdout);
 
905
 
 
906
    return;
 
907
}