Changeset 66601


Ignore:
Timestamp:
Oct 25, 2015, 9:18:55 PM (9 years ago)
Author:
wenzeslaus
Message:

r.in.lidar: read multiple LAS files specified in a text file [news]

Printing and projection check is simply repeated for multiple files.

guisection Required and part of Optional moved to Input and Output.
input option preserved as single file input and made binary.

Related code in main moved together. Files open (and closed) as needed.
Files iterated on the level of points, so still writting the raster
just once.

Custom resolution and data-based extent supported. Base raster
segment reading fixed for the case of using current region.
Data-based extent is still limited by the current region
when base raster is used unless base raster resolution flag
is used.

Follows refactoring in r66593, r66594, r66595 and r66596.

Individual file extents could be preserved to avoid reading the
files for some row groups. Alternativelly, segement-based
reading and writting could be used.

Location:
grass/trunk/raster/r.in.lidar
Files:
1 added
4 edited

Legend:

Unmodified
Added
Removed

  • grass/trunk/raster/r.in.lidar/info.c

    r66594 r66601  
    9191
    9292
    93 int scan_bounds(LASReaderH LAS_reader, int shell_style, int extents,
     93int scan_bounds(LASReaderH LAS_reader, int shell_style, int extents, int update,
    9494                double zscale, struct Cell_head *region)
    9595{
    … …  
    162162                max_y, min_y, max_x, min_x);
    163163    }
     164    else if (update) {
     165        if (min_x < region->west)
     166            region->west = min_x;
     167        if (max_x > region->east)
     168            region->east = max_x;
     169        if (min_y < region->south)
     170            region->south = min_y;
     171        if (max_y > region->north)
     172            region->north = max_y;
     173    }
    164174    else {
    165175        region->east = max_x;

  • grass/trunk/raster/r.in.lidar/local_proto.h

    r66596 r66601  
    5050/* info.c */
    5151void print_lasinfo(LASHeaderH, LASSRSH);
    52 int scan_bounds(LASReaderH, int, int, double, struct Cell_head *);
     52int scan_bounds(LASReaderH, int, int, int, double, struct Cell_head *);
    5353
    5454/* support.c */
    … …  
    7272                          const char *projstr,
    7373                          int override,
    74                           int shellstyle);
     74                          int verbose);
    7575/* raster reading */
    7676int row_array_get_value_row_col(void *array, int arr_row, int arr_col,
    7777                                int cols, RASTER_MAP_TYPE rtype, double *value);
    7878
     79/* multiple files */
     80
     81struct StringList
     82{
     83    int num_items;
     84    int max_items;
     85    char **items;
     86};
     87
     88void string_list_from_file(struct StringList *string_list, char *filename);
     89void string_list_from_one_item(struct StringList *string_list, char *item);
     90void string_list_free(struct StringList *string_list);
     91
    7992#endif /* __LOCAL_PROTO_H__ */

  • grass/trunk/raster/r.in.lidar/main.c

    r66596 r66601  
    7979    struct Option *method_opt, *base_raster_opt, *zrange_opt, *zscale_opt;
    8080    struct Option *trim_opt, *pth_opt, *res_opt;
     81    struct Option *file_list_opt;
    8182    struct Flag *print_flag, *scan_flag, *shell_style, *over_flag, *extents_flag, *intens_flag;
    8283    struct Flag *base_rast_res_flag;
    … …  
    103104        _("Creates a raster map from LAS LiDAR points using univariate statistics.");
    104105
    105     input_opt = G_define_standard_option(G_OPT_F_INPUT);
     106    input_opt = G_define_standard_option(G_OPT_F_BIN_INPUT);
     107    input_opt->required = NO;
    106108    input_opt->label = _("LAS input file");
    107109    input_opt->description = _("LiDAR input files in LAS format (*.las or *.laz)");
     110    input_opt->guisection = _("Input");
    108111
    109112    output_opt = G_define_standard_option(G_OPT_R_OUTPUT);
     113    output_opt->guisection = _("Output");
     114
     115    file_list_opt = G_define_standard_option(G_OPT_F_INPUT);
     116    file_list_opt->key = "file";
     117    file_list_opt->label = _("File containing names of LAS input files");
     118    file_list_opt->description = _("LiDAR input files in LAS format (*.las or *.laz)");
     119    file_list_opt->required = NO;
     120    file_list_opt->guisection = _("Input");
    110121
    111122    method_opt = G_define_option();
    … …  
    184195    res_opt->description =
    185196        _("Output raster resolution");
     197    res_opt->guisection = _("Output");
    186198
    187199    filter_opt = G_define_option();
    … …  
    214226    extents_flag->description =
    215227        _("Extend region extents based on new dataset");
     228    extents_flag->guisection = _("Output");
    216229
    217230    over_flag = G_define_flag();
    … …  
    245258        exit(EXIT_FAILURE);
    246259
     260    struct StringList infiles;
     261
     262    if (file_list_opt->answer) {
     263        if (access(file_list_opt->answer, F_OK) != 0)
     264            G_fatal_error(_("File <%s> does not exist"), file_list_opt->answer);
     265        string_list_from_file(&infiles, file_list_opt->answer);
     266    }
     267    else {
     268        string_list_from_one_item(&infiles, input_opt->answer);
     269    }
    247270
    248271    /* parse input values */
    249     infile = input_opt->answer;
    250272    outmap = output_opt->answer;
    251273
    … …  
    254276    }
    255277
    256     /* Don't crash on cmd line if file not found */
    257     if (access(infile, F_OK) != 0) {
    258         G_fatal_error(_("Input file <%s> does not exist"), infile);
    259     }
    260     /* Open LAS file*/
    261     LAS_reader = LASReader_Create(infile);
    262     if (LAS_reader == NULL) {
    263         G_fatal_error(_("Unable to open file <%s>"), infile);
    264     }
    265    
    266     LAS_header = LASReader_GetHeader(LAS_reader);
    267     if  (LAS_header == NULL) {
    268         G_fatal_error(_("Input file <%s> is not a LAS LiDAR point cloud"),
    269                         infile);
    270     }
    271 
    272     LAS_srs = LASHeader_GetSRS(LAS_header);
    273 
    274     /* Print LAS header */
    275     if (print_flag->answer) {
    276         /* print... */
    277         print_lasinfo(LAS_header, LAS_srs);
    278 
    279         LASSRS_Destroy(LAS_srs);
    280         LASHeader_Destroy(LAS_header);
    281         LASReader_Destroy(LAS_reader);
    282 
    283         exit(EXIT_SUCCESS);
    284     }
     278    /* check zrange and extent relation */
     279    if (scan_flag->answer || extents_flag->answer) {
     280        if (zrange_opt->answer)
     281            G_warning(_("zrange will not be taken into account during scan"));
     282    }
     283
     284    Rast_get_window(&region);
     285    /* G_get_window seems to be unreliable if the location has been changed */
     286    G_get_set_window(&loc_wind);        /* TODO: v.in.lidar uses G_get_default_window() */
     287
     288    estimated_lines = 0;
     289    int i;
     290    for (i = 0; i < infiles.num_items; i++) {
     291        infile = infiles.items[i];
     292        /* don't if file not found */
     293        if (access(infile, F_OK) != 0)
     294            G_fatal_error(_("Input file <%s> does not exist"), infile);
     295        /* Open LAS file*/
     296        LAS_reader = LASReader_Create(infile);
     297        if (LAS_reader == NULL)
     298            G_fatal_error(_("Unable to open file <%s>"), infile);
     299        LAS_header = LASReader_GetHeader(LAS_reader);
     300        if  (LAS_header == NULL) {
     301            G_fatal_error(_("Input file <%s> is not a LAS LiDAR point cloud"),
     302                            infile);
     303        }
     304
     305        LAS_srs = LASHeader_GetSRS(LAS_header);
     306
     307        /* print info or check projection if we are actually importing */
     308        if (print_flag->answer) {
     309            /* print filename when there is more than one file */
     310            if (infiles.num_items > 1)
     311                fprintf(stdout, "File: %s\n", infile);
     312            /* Print LAS header */
     313            print_lasinfo(LAS_header, LAS_srs);
     314        }
     315        else {
     316            /* report that we are checking more files */
     317            if (i == 1)
     318                G_message(_("First file's projection checked,"
     319                            " checking projection of the other files..."));
     320            /* Fetch input map projection in GRASS form. */
     321            projstr = LASSRS_GetWKT_CompoundOK(LAS_srs);
     322            /* we are printing the non-warning messages only for first file */
     323            projection_check_wkt(cellhd, loc_wind, projstr, over_flag->answer,
     324                                 shell_style->answer || i);
     325            /* if there is a problem in some other file, first OK message
     326             * is printed but than a warning, this is not ideal but hopefully
     327             * not so confusing when importing multiple files */
     328        }
     329        if (scan_flag->answer || extents_flag->answer) {
     330            /* we assign to the first one (i==0) but update for the rest */
     331            scan_bounds(LAS_reader, shell_style->answer, extents_flag->answer, i,
     332                        zscale, &region);
     333        }
     334        /* number of estimated point across all files */
     335        /* TODO: this should be ull which won't work with percent report */
     336        estimated_lines += LASHeader_GetPointRecordsCount(LAS_header);
     337        /* We are closing all again and we will be opening them later,
     338         * so we don't have to worry about limit for open files. */
     339        LASSRS_Destroy(LAS_srs);
     340        LASHeader_Destroy(LAS_header);
     341        LASReader_Destroy(LAS_reader);
     342    }
     343    /* if we are not importing, end */
     344    if (print_flag->answer || scan_flag->answer)
     345        exit(EXIT_SUCCESS);
    285346
    286347    return_filter = LAS_ALL;
    … …  
    300361    class_filter_create_from_strings(&class_filter, class_opt->answers);
    301362
    302     /* Fetch input map projection in GRASS form. */
    303     projstr = LASSRS_GetWKT_CompoundOK(LAS_srs);
    304 
    305     if (TRUE) {
    306         projection_check_wkt(cellhd, loc_wind, projstr, over_flag->answer, shell_style->answer);
    307     }
    308 
    309363    percent = atoi(percent_opt->answer);
    310364    zscale = atof(zscale_opt->answer);
    … …  
    339393        rtype = CELL_TYPE;
    340394
    341 
    342     Rast_get_window(&region);
    343 
    344 
    345     if (scan_flag->answer || extents_flag->answer) {
    346         if (zrange_opt->answer)
    347             G_warning(_("zrange will not be taken into account during scan"));
    348 
    349         scan_bounds(LAS_reader, shell_style->answer, extents_flag->answer,
    350                     zscale, &region);
    351 
    352         if (!extents_flag->answer) {
    353             LASHeader_Destroy(LAS_header);
    354             LASReader_Destroy(LAS_reader);
    355 
    356             exit(EXIT_SUCCESS);
    357         }
    358     }
    359    
    360395    if (res_opt->answer) {
    361396        /* align to resolution */
    … …  
    397432    int use_segment = 0;
    398433    int use_base_raster_res = 0;
     434    /* TODO: see if the input region extent is smaller than the raster
     435     * if yes, the we need to load the whole base raster if the -e
     436     * flag was defined (alternatively clip the regions) */
    399437    if (base_rast_res_flag->answer)
    400438        use_base_raster_res = 1;
    401     if (base_raster_opt->answer && (res_opt->answer || use_base_raster_res))
     439    if (base_raster_opt->answer && (res_opt->answer || use_base_raster_res
     440                                    || extents_flag->answer))
    402441        use_segment = 1;
    403442    if (base_raster_opt->answer && !use_segment) {
    … …  
    409448    if (base_raster_opt->answer && use_segment) {
    410449        if (use_base_raster_res) {
    411             /* TODO: how to get cellhd already stored in the open map? */
     450            /* read raster actual extent and resolution */
    412451            Rast_get_cellhd(base_raster_opt->answer, "", &input_region);
    413452            /* TODO: make it only as small as the output is or points are */
    414453            Rast_set_input_window(&input_region);  /* we have split window */
     454        } else {
     455            Rast_get_input_window(&input_region);
    415456        }
    416457        rast_segment_open(&base_segment, base_raster_opt->answer, &base_raster_data_type);
    … …  
    428469    out_fd = Rast_open_new(outmap, rtype);
    429470
    430     estimated_lines = LASHeader_GetPointRecordsCount(LAS_header);
    431 
    432471    /* allocate memory for a single row of output data */
    433472    raster_row = Rast_allocate_output_buf(rtype);
    … …  
    442481    /* main binning loop(s) */
    443482    for (pass = 1; pass <= npasses; pass++) {
    444         LASError LAS_error;
    445483
    446484        if (npasses > 1)
    447485            G_message(_("Pass #%d (of %d) ..."), pass, npasses);
    448 
    449         LAS_error = LASReader_Seek(LAS_reader, 0);
    450        
    451         if (LAS_error != LE_None)
    452             G_fatal_error(_("Could not rewind input file"));
    453486
    454487        /* figure out segmentation */
    … …  
    477510        G_percent_reset();
    478511
    479         while ((LAS_point = LASReader_GetNextPoint(LAS_reader)) != NULL) {
    480             line++;
    481             counter++;
    482 
    483             if (counter == 100000) {    /* speed */
    484                 if (line < estimated_lines)
    485                     G_percent(line, estimated_lines, 3);
    486                 counter = 0;
    487             }
    488 
    489             if (!LASPoint_IsValid(LAS_point)) {
    490                 continue;
    491             }
    492 
    493             x = LASPoint_GetX(LAS_point);
    494             y = LASPoint_GetY(LAS_point);
    495             if (intens_flag->answer)
    496                 /* use z variable here to allow for scaling of intensity below */
    497                 z = LASPoint_GetIntensity(LAS_point);
    498             else
    499                 z = LASPoint_GetZ(LAS_point);
    500 
    501         int return_n = LASPoint_GetReturnNumber(LAS_point);
    502         int n_returns = LASPoint_GetNumberOfReturns(LAS_point);
    503         if (return_filter_is_out(&return_filter_struct, return_n, n_returns)) {
    504             n_filtered++;
    505             continue;
    506         }
    507         point_class = (int) LASPoint_GetClassification(LAS_point);
    508         if (class_filter_is_out(&class_filter, point_class))
    509             continue;
    510 
    511             if (y <= pass_south || y > pass_north) {
    512                 continue;
    513             }
    514             if (x < region.west || x >= region.east) {
    515                 continue;
    516             }
    517 
    518             z = z * zscale;
    519 
    520             /* find the bin in the current array box */
    521             arr_row = (int)((pass_north - y) / region.ns_res);
    522             arr_col = (int)((x - region.west) / region.ew_res);
    523 
    524             if (base_array) {
    525                 double base_z;
    526                 if (row_array_get_value_row_col(base_array, arr_row, arr_col,
    527                                                 cols, base_raster_data_type,
    528                                                 &base_z))
    529                     z -= base_z;
     512        /* loop of input files */
     513        for (i = 0; i < infiles.num_items; i++) {
     514            infile = infiles.items[i];
     515            /* we already know file is there, so just do basic checks */
     516            LAS_reader = LASReader_Create(infile);
     517            if (LAS_reader == NULL)
     518                G_fatal_error(_("Unable to open file <%s>"), infile);
     519
     520            while ((LAS_point = LASReader_GetNextPoint(LAS_reader)) != NULL) {
     521                line++;
     522                counter++;
     523
     524                if (counter == 100000) {        /* speed */
     525                    if (line < estimated_lines)
     526                        G_percent(line, estimated_lines, 3);
     527                    counter = 0;
     528                }
     529
     530                if (!LASPoint_IsValid(LAS_point)) {
     531                    continue;
     532                }
     533
     534                x = LASPoint_GetX(LAS_point);
     535                y = LASPoint_GetY(LAS_point);
     536                if (intens_flag->answer)
     537                    /* use z variable here to allow for scaling of intensity below */
     538                    z = LASPoint_GetIntensity(LAS_point);
    530539                else
     540                    z = LASPoint_GetZ(LAS_point);
     541
     542                int return_n = LASPoint_GetReturnNumber(LAS_point);
     543                int n_returns = LASPoint_GetNumberOfReturns(LAS_point);
     544                if (return_filter_is_out(&return_filter_struct, return_n, n_returns)) {
     545                    n_filtered++;
    531546                    continue;
    532             }
    533             else if (use_segment) {
    534                 double base_z;
    535                 if (rast_segment_get_value_xy(&base_segment, &input_region,
    536                                               base_raster_data_type, x, y,
    537                                               &base_z))
    538                     z -= base_z;
    539                 else
     547                }
     548                point_class = (int) LASPoint_GetClassification(LAS_point);
     549                if (class_filter_is_out(&class_filter, point_class))
    540550                    continue;
    541             }
    542 
    543             if (zrange_opt->answer) {
    544                 if (z < zrange_min || z > zrange_max) {
     551
     552                if (y <= pass_south || y > pass_north) {
    545553                    continue;
    546554                }
    547             }
    548 
    549             count++;
    550             /*          G_debug(5, "x: %f, y: %f, z: %f", x, y, z); */
    551 
    552             update_value(&point_binning, &bin_index_nodes, cols, arr_row, arr_col, rtype, z);
    553         }                       /* while !EOF */
     555                if (x < region.west || x >= region.east) {
     556                    continue;
     557                }
     558
     559                z = z * zscale;
     560
     561                /* find the bin in the current array box */
     562                arr_row = (int)((pass_north - y) / region.ns_res);
     563                arr_col = (int)((x - region.west) / region.ew_res);
     564
     565                if (base_array) {
     566                    double base_z;
     567                    if (row_array_get_value_row_col(base_array, arr_row, arr_col,
     568                                                    cols, base_raster_data_type,
     569                                                    &base_z))
     570                        z -= base_z;
     571                    else
     572                        continue;
     573                }
     574                else if (use_segment) {
     575                    double base_z;
     576                    if (rast_segment_get_value_xy(&base_segment, &input_region,
     577                                                  base_raster_data_type, x, y,
     578                                                  &base_z))
     579                        z -= base_z;
     580                    else
     581                        continue;
     582                }
     583
     584                if (zrange_opt->answer) {
     585                    if (z < zrange_min || z > zrange_max) {
     586                        continue;
     587                    }
     588                }
     589
     590                count++;
     591                /*          G_debug(5, "x: %f, y: %f, z: %f", x, y, z); */
     592
     593                update_value(&point_binning, &bin_index_nodes, cols, arr_row, arr_col, rtype, z);
     594            }                        /* while !EOF of one input file */
     595            /* close input LAS file */
     596            LASReader_Destroy(LAS_reader);
     597        }           /* end of loop for all input files files */
    554598
    555599        G_percent(1, 1, 1);     /* flush */
    … …  
    568612        /* free memory */
    569613        point_binning_free(&point_binning, &bin_index_nodes);
     614    string_list_free(&infiles);
    570615    }                           /* passes loop */
    571616    if (base_array)
    … …  
    577622    G_free(raster_row);
    578623
    579     /* close input LAS file */
    580     LASHeader_Destroy(LAS_header);
    581     LASReader_Destroy(LAS_reader);
    582 
    583624    /* close raster file & write history */
    584625    Rast_close(out_fd);

  • grass/trunk/raster/r.in.lidar/projection.c

    r66594 r66601  
    100100void projection_check_wkt(struct Cell_head cellhd,
    101101                          struct Cell_head loc_wind,
    102                           const char *projstr, int override, int shellstyle)
     102                          const char *projstr, int override, int verbose)
    103103{
    104104    struct Key_Value *loc_proj_info = NULL, *loc_proj_units = NULL;
    … …  
    115115
    116116    /* Does the projection of the current location match the dataset? */
    117     /* G_get_window seems to be unreliable if the location has been changed */
    118     G_get_set_window(&loc_wind);        /* TODO: v.in.lidar uses G_get_default_window() */
     117
    119118    /* fetch LOCATION PROJ info */
    120119    if (loc_wind.proj != PROJECTION_XY) {
    … …  
    126125        cellhd.proj = loc_wind.proj;
    127126        cellhd.zone = loc_wind.zone;
    128         G_message(_("Over-riding projection check"));
     127        if (verbose)
     128            G_message(_("Over-riding projection check"));
    129129    }
    130130    else if (loc_wind.proj != cellhd.proj
    … …  
    136136                                   proj_info, proj_units, err);
    137137    }
    138     else if (!shellstyle) {
     138    else if (verbose) {
    139139        G_message(_("Projection of input dataset and current location "
    140140                    "appear to match"));
Note: See TracChangeset for help on using the changeset viewer.