[GRASS-SVN] r46570 - in grass/trunk/raster: . r.in.lidar
svn_grass at osgeo.org
svn_grass at osgeo.org
Sun Jun 5 05:07:11 EDT 2011
Author: mmetz
Date: 2011-06-05 02:07:11 -0700 (Sun, 05 Jun 2011)
New Revision: 46570
Added:
grass/trunk/raster/r.in.lidar/
grass/trunk/raster/r.in.lidar/r.in.lidar.html
Removed:
grass/trunk/raster/r.in.lidar/r.in.xyz.html
Modified:
grass/trunk/raster/r.in.lidar/Makefile
grass/trunk/raster/r.in.lidar/local_proto.h
grass/trunk/raster/r.in.lidar/main.c
grass/trunk/raster/r.in.lidar/support.c
Log:
add LAS LiDAR version of r.in.xyz
Modified: grass/trunk/raster/r.in.lidar/Makefile
===================================================================
--- grass/trunk/raster/r.in.xyz/Makefile 2011-06-01 06:31:05 UTC (rev 46482)
+++ grass/trunk/raster/r.in.lidar/Makefile 2011-06-05 09:07:11 UTC (rev 46570)
@@ -1,10 +1,13 @@
MODULE_TOPDIR = ../..
-PGM = r.in.xyz
+PGM = r.in.lidar
-LIBES = $(RASTERLIB) $(GISLIB) $(MATHLIB)
+LIBES = $(RASTERLIB) $(GISLIB) $(MATHLIB) $(GPROJLIB) $(LASLIBS)
DEPENDENCIES = $(RASTERDEP) $(GISDEP)
include $(MODULE_TOPDIR)/include/Make/Module.make
+ifneq ($(USE_LIBLAS),)
default: cmd
+endif
+
Modified: grass/trunk/raster/r.in.lidar/local_proto.h
===================================================================
--- grass/trunk/raster/r.in.xyz/local_proto.h 2011-06-01 06:31:05 UTC (rev 46482)
+++ grass/trunk/raster/r.in.lidar/local_proto.h 2011-06-05 09:07:11 UTC (rev 46570)
@@ -1,24 +1,27 @@
-/*
- * r.in.xyz
+/****************************************************************************
*
- * Calculates univariate statistics from the non-null cells of a GRASS
- * raster map
+ * MODULE: r.in.Lidar
+ *
+ * AUTHOR(S): Markus Metz
+ * Based on r.in.xyz by Hamish Bowman, Volker Wichmann
*
- * Copyright 2006 by M. Hamish Bowman, and The GRASS Development Team
- * Author: M. Hamish Bowman, University of Otago, Dunedin, New Zealand
+ * PURPOSE: Imports LAS LiDAR point clouds to a raster map using
+ * aggregate statistics.
*
- * This program is free software licensed under the GPL (>=v2).
- * Read the COPYING file that comes with GRASS for details.
+ * COPYRIGHT: (C) 2011 Markus Metz and the The GRASS Development Team
*
- * This program is intended as a replacement for the GRASS 5
- * s.cellstats module.
- */
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ *****************************************************************************/
#ifndef __LOCAL_PROTO_H__
#define __LOCAL_PROTO_H__
#include <grass/gis.h>
+#include <liblas/capi/liblas.h>
#define BUFFSIZE 256
@@ -38,7 +41,7 @@
#define METHOD_TRIMMEAN 13
/* main.c */
-int scan_bounds(FILE *, int, int, int, char *, int, int, double);
+int scan_bounds(LASReaderH, int, double);
/* support.c */
int blank_array(void *, int, int, RASTER_MAP_TYPE, int);
Modified: grass/trunk/raster/r.in.lidar/main.c
===================================================================
--- grass/trunk/raster/r.in.xyz/main.c 2011-06-01 06:31:05 UTC (rev 46482)
+++ grass/trunk/raster/r.in.lidar/main.c 2011-06-05 09:07:11 UTC (rev 46570)
@@ -1,21 +1,21 @@
-/*
- * r.in.xyz
+/****************************************************************************
*
- * Calculates univariate statistics from the non-null cells of a GRASS raster map
+ * MODULE: r.in.Lidar
+ *
+ * AUTHOR(S): Markus Metz
+ * Based on r.in.xyz by Hamish Bowman, Volker Wichmann
*
- * Copyright 2006 by M. Hamish Bowman, and The GRASS Development Team
- * Author: M. Hamish Bowman, University of Otago, Dunedin, New Zealand
+ * PURPOSE: Imports LAS LiDAR point clouds to a raster map using
+ * aggregate statistics.
*
- * Extended 2007 by Volker Wichmann to support the aggregate functions
- * median, percentile, skewness and trimmed mean
+ * COPYRIGHT: (C) 2011 Markus Metz and the The GRASS Development Team
*
- * This program is free software licensed under the GPL (>=v2).
- * Read the COPYING file that comes with GRASS for details.
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
*
- * This program is intended as a replacement for the GRASS 5 s.cellstats module.
- */
+ *****************************************************************************/
-#include <grass/config.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
@@ -23,6 +23,7 @@
#include <sys/types.h>
#include <grass/gis.h>
#include <grass/raster.h>
+#include <grass/gprojects.h>
#include <grass/glocale.h>
#include "local_proto.h"
@@ -94,21 +95,14 @@
int main(int argc, char *argv[])
{
-
- FILE *in_fp;
int out_fd;
char *infile, *outmap;
- int xcol, ycol, zcol, max_col, percent;
+ int percent;
int do_zfilter;
int method = -1;
int bin_n, bin_min, bin_max, bin_sum, bin_sumsq, bin_index;
double zrange_min, zrange_max, d_tmp;
- char *fs; /* field delim */
- off_t filesize;
- int linesize;
long estimated_lines;
- int from_stdin;
- int can_seek;
RASTER_MAP_TYPE rtype;
struct History history;
@@ -121,11 +115,10 @@
int row, col; /* counters */
int pass, npasses;
- unsigned long line;
+ unsigned long line, line_total;
+ unsigned int counter;
char buff[BUFFSIZE];
double x, y, z;
- char **tokens;
- int ntokens; /* number of tokens */
double pass_north, pass_south;
int arr_row, arr_col;
unsigned long count, count_total;
@@ -146,14 +139,22 @@
int r_low, r_up;
struct GModule *module;
- struct Option *input_opt, *output_opt, *delim_opt, *percent_opt,
- *type_opt;
- struct Option *method_opt, *xcol_opt, *ycol_opt, *zcol_opt, *zrange_opt,
- *zscale_opt;
+ struct Option *input_opt, *output_opt, *percent_opt, *type_opt;
+ struct Option *method_opt, *zrange_opt, *zscale_opt;
struct Option *trim_opt, *pth_opt;
- struct Flag *scan_flag, *shell_style, *skipline;
+ struct Flag *scan_flag, *shell_style, *over_flag;
+ LASReaderH LAS_reader;
+ LASHeaderH LAS_header;
+ LASSRSH LAS_srs;
+ LASPointH LAS_point;
+ struct Key_Value *loc_proj_info = NULL, *loc_proj_units = NULL;
+ struct Key_Value *proj_info, *proj_units;
+ const char *projstr;
+ struct Cell_head cellhd, loc_wind;
+
+
G_gisinit(argv[0]);
module = G_define_module();
@@ -161,11 +162,11 @@
G_add_keyword(_("import"));
G_add_keyword(_("LIDAR"));
module->description =
- _("Create a raster map from an assemblage of many coordinates using univariate statistics.");
+ _("Create a raster map from LAS LiDAR points using univariate statistics.");
input_opt = G_define_standard_option(G_OPT_F_INPUT);
input_opt->description =
- _("ASCII file containing input data (or \"-\" to read from stdin)");
+ _("LiDAR LAS input file");
output_opt = G_define_standard_option(G_OPT_R_OUTPUT);
@@ -187,34 +188,6 @@
type_opt->answer = "FCELL";
type_opt->description = _("Storage type for resultant raster map");
- delim_opt = G_define_standard_option(G_OPT_F_SEP);
- delim_opt->guisection = _("Input");
-
- xcol_opt = G_define_option();
- xcol_opt->key = "x";
- xcol_opt->type = TYPE_INTEGER;
- xcol_opt->required = NO;
- xcol_opt->answer = "1";
- xcol_opt->description =
- _("Column number of x coordinates in input file (first column is 1)");
- xcol_opt->guisection = _("Input");
-
- ycol_opt = G_define_option();
- ycol_opt->key = "y";
- ycol_opt->type = TYPE_INTEGER;
- ycol_opt->required = NO;
- ycol_opt->answer = "2";
- ycol_opt->description = _("Column number of y coordinates in input file");
- ycol_opt->guisection = _("Input");
-
- zcol_opt = G_define_option();
- zcol_opt->key = "z";
- zcol_opt->type = TYPE_INTEGER;
- zcol_opt->required = NO;
- zcol_opt->answer = "3";
- zcol_opt->description = _("Column number of data values in input file");
- zcol_opt->guisection = _("Input");
-
zrange_opt = G_define_option();
zrange_opt->key = "zrange";
zrange_opt->type = TYPE_DOUBLE;
@@ -256,6 +229,11 @@
_("Discard <trim> percent of the smallest and <trim> percent of the largest observations");
trim_opt->guisection = _("Statistic");
+ over_flag = G_define_flag();
+ over_flag->key = 'o';
+ over_flag->description =
+ _("Override dataset projection (use location's projection)");
+
scan_flag = G_define_flag();
scan_flag->key = 's';
scan_flag->description = _("Scan data file for extent then exit");
@@ -265,10 +243,6 @@
shell_style->description =
_("In scan mode, print using shell script style");
- skipline = G_define_flag();
- skipline->key = 'i';
- skipline->description = _("Ignore broken lines");
-
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
@@ -281,22 +255,132 @@
scan_flag->answer = 1; /* pointer not int, so set = shell_style->answer ? */
}
- fs = delim_opt->answer;
- if (strcmp(fs, "\\t") == 0)
- fs = "\t";
- if (strcmp(fs, "tab") == 0)
- fs = "\t";
- if (strcmp(fs, "space") == 0)
- fs = " ";
+ /* Open LAS file*/
+ LAS_reader = LASReader_Create(infile);
+ LAS_header = LASReader_GetHeader(LAS_reader);
- xcol = atoi(xcol_opt->answer);
- ycol = atoi(ycol_opt->answer);
- zcol = atoi(zcol_opt->answer);
- if ((xcol < 0) || (ycol < 0) || (zcol < 0))
- G_fatal_error(_("Please specify a reasonable column number."));
- max_col = (xcol > ycol) ? xcol : ycol;
- max_col = (zcol > max_col) ? zcol : max_col;
+ if (LAS_header == NULL) {
+ G_fatal_error(_("Input file <%s> is not a LAS LiDAR point cloud"),
+ infile);
+ }
+ LAS_srs = LASHeader_GetSRS(LAS_header);
+
+ /* Fetch input map projection in GRASS form. */
+ proj_info = NULL;
+ proj_units = NULL;
+ projstr = LASSRS_GetWKT_CompoundOK(LAS_srs);
+
+ if (1) {
+ int err = 0;
+ char error_msg[8192];
+
+ /* Projection only required for checking so convert non-interactively */
+ if (GPJ_wkt_to_grass(&cellhd, &proj_info,
+ &proj_units, projstr, 0) < 0)
+ G_warning(_("Unable to convert input map projection information to "
+ "GRASS format for checking"));
+
+ /* Does the projection of the current location match the dataset? */
+ /* G_get_window seems to be unreliable if the location has been changed */
+ G_get_set_window(&loc_wind);
+ /* fetch LOCATION PROJ info */
+ if (loc_wind.proj != PROJECTION_XY) {
+ loc_proj_info = G_get_projinfo();
+ loc_proj_units = G_get_projunits();
+ }
+
+ if (over_flag->answer) {
+ cellhd.proj = loc_wind.proj;
+ cellhd.zone = loc_wind.zone;
+ G_message(_("Over-riding projection check"));
+ }
+ else if (loc_wind.proj != cellhd.proj
+ || (err =
+ G_compare_projections(loc_proj_info, loc_proj_units,
+ proj_info, proj_units)) != TRUE) {
+ int i_value;
+
+ strcpy(error_msg,
+ _("Projection of dataset does not"
+ " appear to match current location.\n\n"));
+
+ /* TODO: output this info sorted by key: */
+ if (loc_wind.proj != cellhd.proj || err != -2) {
+ if (loc_proj_info != NULL) {
+ strcat(error_msg, _("GRASS LOCATION PROJ_INFO is:\n"));
+ for (i_value = 0; i_value < loc_proj_info->nitems;
+ i_value++)
+ sprintf(error_msg + strlen(error_msg), "%s: %s\n",
+ loc_proj_info->key[i_value],
+ loc_proj_info->value[i_value]);
+ strcat(error_msg, "\n");
+ }
+
+ if (proj_info != NULL) {
+ strcat(error_msg, _("Import dataset PROJ_INFO is:\n"));
+ for (i_value = 0; i_value < proj_info->nitems; i_value++)
+ sprintf(error_msg + strlen(error_msg), "%s: %s\n",
+ proj_info->key[i_value],
+ proj_info->value[i_value]);
+ }
+ else {
+ strcat(error_msg, _("Import dataset PROJ_INFO is:\n"));
+ if (cellhd.proj == PROJECTION_XY)
+ sprintf(error_msg + strlen(error_msg),
+ "Dataset proj = %d (unreferenced/unknown)\n",
+ cellhd.proj);
+ else if (cellhd.proj == PROJECTION_LL)
+ sprintf(error_msg + strlen(error_msg),
+ "Dataset proj = %d (lat/long)\n",
+ cellhd.proj);
+ else if (cellhd.proj == PROJECTION_UTM)
+ sprintf(error_msg + strlen(error_msg),
+ "Dataset proj = %d (UTM), zone = %d\n",
+ cellhd.proj, cellhd.zone);
+ else if (cellhd.proj == PROJECTION_SP)
+ sprintf(error_msg + strlen(error_msg),
+ "Dataset proj = %d (State Plane), zone = %d\n",
+ cellhd.proj, cellhd.zone);
+ else
+ sprintf(error_msg + strlen(error_msg),
+ "Dataset proj = %d (unknown), zone = %d\n",
+ cellhd.proj, cellhd.zone);
+ }
+ }
+ else {
+ if (loc_proj_units != NULL) {
+ strcat(error_msg, "GRASS LOCATION PROJ_UNITS is:\n");
+ for (i_value = 0; i_value < loc_proj_units->nitems;
+ i_value++)
+ sprintf(error_msg + strlen(error_msg), "%s: %s\n",
+ loc_proj_units->key[i_value],
+ loc_proj_units->value[i_value]);
+ strcat(error_msg, "\n");
+ }
+
+ if (proj_units != NULL) {
+ strcat(error_msg, "Import dataset PROJ_UNITS is:\n");
+ for (i_value = 0; i_value < proj_units->nitems; i_value++)
+ sprintf(error_msg + strlen(error_msg), "%s: %s\n",
+ proj_units->key[i_value],
+ proj_units->value[i_value]);
+ }
+ }
+ sprintf(error_msg + strlen(error_msg),
+ _("\nYou can use the -o flag to %s to override this projection check.\n"),
+ G_program_name());
+ strcat(error_msg,
+ _("Consider generating a new location with 'location' parameter"
+ " from input data set.\n"));
+ G_fatal_error(error_msg);
+ }
+ else {
+ G_message(_("Projection of input dataset and current location "
+ "appear to match"));
+ }
+ }
+
percent = atoi(percent_opt->answer);
zscale = atof(zscale_opt->answer);
@@ -339,6 +423,13 @@
bin_sumsq = FALSE;
bin_index = FALSE;
+ n_array = NULL;
+ min_array = NULL;
+ max_array = NULL;
+ sum_array = NULL;
+ sumsq_array = NULL;
+ index_array = NULL;
+
if (strcmp(method_opt->answer, "n") == 0) {
method = METHOD_N;
bin_n = TRUE;
@@ -419,7 +510,7 @@
rtype = CELL_TYPE;
- G_get_window(®ion);
+ Rast_get_window(®ion);
rows = (int)(region.rows * (percent / 100.0));
cols = region.cols;
@@ -464,34 +555,14 @@
}
- /* open input file */
- if (strcmp("-", infile) == 0) {
- from_stdin = TRUE;
- in_fp = stdin;
- infile = G_store("stdin"); /* filename for history metadata */
- }
- else {
- if ((in_fp = fopen(infile, "r")) == NULL)
- G_fatal_error(_("Unable to open input file <%s>"), infile);
- }
-
- can_seek = fseek(in_fp, 0L, SEEK_SET) == 0;
-
- /* can't rewind() non-files */
- if (!can_seek && npasses != 1) {
- G_warning(_("If input is not from a file it is only possible to perform a single pass."));
- npasses = 1;
- }
-
if (scan_flag->answer) {
if (zrange_opt->answer)
G_warning(_("zrange will not be taken into account during scan"));
- scan_bounds(in_fp, xcol, ycol, zcol, fs, shell_style->answer,
- skipline->answer, zscale);
+ scan_bounds(LAS_reader, shell_style->answer, zscale);
- if (!from_stdin)
- fclose(in_fp);
+ LASHeader_Destroy(LAS_header);
+ LASReader_Destroy(LAS_reader);
exit(EXIT_SUCCESS);
}
@@ -500,44 +571,37 @@
/* open output map */
out_fd = Rast_open_new(outmap, rtype);
- if (can_seek) {
- /* guess at number of lines in the file without actually reading it all in */
- for (line = 0; line < 10; line++) { /* arbitrarily use 10th line for guess */
- if (0 == G_getl2(buff, BUFFSIZE - 1, in_fp))
- break;
- linesize = strlen(buff) + 1;
- }
- G_fseek(in_fp, 0L, SEEK_END);
- filesize = G_ftell(in_fp);
- rewind(in_fp);
- if (linesize < 6) /* min possible: "0,0,0\n" */
- linesize = 6;
- estimated_lines = filesize / linesize;
- G_debug(2, "estimated number of lines in file: %ld", estimated_lines);
- }
- else
- estimated_lines = -1;
+ estimated_lines = LASHeader_GetPointRecordsCount(LAS_header);
/* allocate memory for a single row of output data */
raster_row = Rast_allocate_buf(rtype);
G_message(_("Reading data ..."));
- count_total = 0;
+ count_total = line_total = 0;
+ /* init northern border */
+ pass_south = region.north;
+
/* main binning loop(s) */
for (pass = 1; pass <= npasses; pass++) {
+ LASError LAS_error;
+
if (npasses > 1)
G_message(_("Pass #%d (of %d) ..."), pass, npasses);
- if (can_seek)
- rewind(in_fp);
+ LAS_error = LASReader_Seek(LAS_reader, 0);
+
+ if (LAS_error != LE_None)
+ G_fatal_error(_("Could not rewind input file"));
/* figure out segmentation */
- pass_north = region.north - (pass - 1) * rows * region.ns_res;
- if (pass == npasses)
+ pass_north = pass_south; /* exact copy to avoid fp errors */
+ pass_south = region.north - pass * rows * region.ns_res;
+ if (pass == npasses) {
rows = region.rows - (pass - 1) * rows;
- pass_south = pass_north - rows * region.ns_res;
+ pass_south = region.south; /* exact copy to avoid fp errors */
+ }
G_debug(2, "pass=%d/%d pass_n=%f pass_s=%f rows=%d",
pass, npasses, pass_north, pass_south, rows);
@@ -577,43 +641,24 @@
line = 0;
count = 0;
+ counter = 0;
G_percent_reset();
- while (0 != G_getl2(buff, BUFFSIZE - 1, in_fp)) {
+ while ((LAS_point = LASReader_GetNextPoint(LAS_reader)) != NULL) {
line++;
+ counter++;
- if (line % 10000 == 0) { /* mod for speed */
- if (!can_seek)
- G_clicker();
- else if (line < estimated_lines)
+ if (counter == 10000) { /* speed */
+ if (line < estimated_lines)
G_percent(line, estimated_lines, 3);
+ counter = 0;
}
- if ((buff[0] == '#') || (buff[0] == '\0')) {
- continue; /* line is a comment or blank */
+ if (!LASPoint_IsValid(LAS_point)) {
+ continue;
}
- G_chop(buff); /* remove leading and trailing whitespace from the string. unneded?? */
- tokens = G_tokenize(buff, fs);
- ntokens = G_number_of_tokens(tokens);
- if ((ntokens < 3) || (max_col > ntokens)) {
- if (skipline->answer) {
- G_warning(_("Not enough data columns. "
- "Incorrect delimiter or column number? "
- "Found the following character(s) in row %lu:\n[%s]"),
- line, buff);
- G_warning(_("Line ignored as requested"));
- continue; /* line is garbage */
- }
- else {
- G_fatal_error(_("Not enough data columns. "
- "Incorrect delimiter or column number? "
- "Found the following character(s) in row %lu:\n[%s]"),
- line, buff);
- }
- }
-
/* too slow?
if ( G_projection() == PROJECTION_LL ) {
G_scan_easting( tokens[xcol-1], &x, region.proj);
@@ -621,36 +666,27 @@
}
else {
*/
- if (1 != sscanf(tokens[ycol - 1], "%lf", &y))
- G_fatal_error(_("Bad y-coordinate line %lu column %d. <%s>"),
- line, ycol, tokens[ycol - 1]);
+ x = LASPoint_GetX(LAS_point);
+ y = LASPoint_GetY(LAS_point);
+ z = LASPoint_GetZ(LAS_point);
+
if (y <= pass_south || y > pass_north) {
- G_free_tokens(tokens);
continue;
}
- if (1 != sscanf(tokens[xcol - 1], "%lf", &x))
- G_fatal_error(_("Bad x-coordinate line %lu column %d. <%s>"),
- line, xcol, tokens[xcol - 1]);
if (x < region.west || x > region.east) {
- G_free_tokens(tokens);
continue;
}
- if (1 != sscanf(tokens[zcol - 1], "%lf", &z))
- G_fatal_error(_("Bad z-coordinate line %lu column %d. <%s>"),
- line, zcol, tokens[zcol - 1]);
z = z * zscale;
if (zrange_opt->answer) {
if (z < zrange_min || z > zrange_max) {
- G_free_tokens(tokens);
continue;
}
}
count++;
/* G_debug(5, "x: %f, y: %f, z: %f", x, y, z); */
- G_free_tokens(tokens);
/* find the bin in the current array box */
arr_row = (int)((pass_north - y) / region.ns_res);
@@ -715,8 +751,8 @@
G_percent(1, 1, 1); /* flush */
G_debug(2, "pass %d finished, %lu coordinates in box", pass, count);
count_total += count;
+ line_total += line;
-
/* calc stats and output */
G_message(_("Writing to map ..."));
for (row = 0; row < rows; row++) {
@@ -788,23 +824,31 @@
if (n == 0)
Rast_set_null_value(ptr, 1, rtype);
+ else if (n == 1)
+ Rast_set_d_value(ptr, 0.0, rtype);
else {
variance = (sumsq - sum * sum / n) / n;
if (variance < GRASS_EPSILON)
variance = 0.0;
- if (method == METHOD_STDDEV)
- Rast_set_d_value(ptr, sqrt(variance), rtype);
+ /* nan test */
+ if (variance != variance)
+ Rast_set_null_value(ptr, 1, rtype);
+ else {
- else if (method == METHOD_VARIANCE)
+ if (method == METHOD_STDDEV)
+ variance = sqrt(variance);
+
+ else if (method == METHOD_COEFF_VAR)
+ variance = 100 * sqrt(variance) / (sum / n);
+
+ /* nan test */
+ if (variance != variance)
+ variance = 0.0; /* OK for n > 0 ?*/
+
Rast_set_d_value(ptr, variance, rtype);
+ }
- else if (method == METHOD_COEFF_VAR)
- Rast_set_d_value(ptr,
- 100 * sqrt(variance) / (sum /
- n),
- rtype);
-
}
ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
}
@@ -1039,9 +1083,9 @@
G_percent(1, 1, 1); /* flush */
G_free(raster_row);
- /* close input file */
- if (!from_stdin)
- fclose(in_fp);
+ /* close input LAS file */
+ LASHeader_Destroy(LAS_header);
+ LASReader_Destroy(LAS_reader);
/* close raster file & write history */
Rast_close(out_fd);
@@ -1058,7 +1102,7 @@
sprintf(buff, _("%lu points found in region."), count_total);
G_done_msg(buff);
- G_debug(1, "Processed %lu lines.", line);
+ G_debug(1, "Processed %lu points.", line_total);
exit(EXIT_SUCCESS);
@@ -1066,64 +1110,35 @@
-int scan_bounds(FILE * fp, int xcol, int ycol, int zcol, char *fs,
- int shell_style, int skipline, double zscale)
+int scan_bounds(LASReaderH LAS_reader, int shell_style, double zscale)
{
unsigned long line;
- int first, max_col;
- char buff[BUFFSIZE];
+ int first;
double min_x, max_x, min_y, max_y, min_z, max_z;
- char **tokens;
- int ntokens; /* number of tokens */
double x, y, z;
+ LASPointH LAS_point;
- max_col = (xcol > ycol) ? xcol : ycol;
- max_col = (zcol > max_col) ? zcol : max_col;
-
line = 0;
first = TRUE;
+
+ /* init to nan in case no points are found */
+ min_x = max_x = min_y = max_y = min_z = max_z = 0.0 / 0.0;
G_verbose_message(_("Scanning data ..."));
+
+ LASReader_Seek(LAS_reader, 0);
- while (0 != G_getl2(buff, BUFFSIZE - 1, fp)) {
+ while ((LAS_point = LASReader_GetNextPoint(LAS_reader)) != NULL) {
line++;
- if ((buff[0] == '#') || (buff[0] == '\0')) {
- continue; /* line is a comment or blank */
+ if (!LASPoint_IsValid(LAS_point)) {
+ continue;
}
- G_chop(buff); /* remove leading and trailing whitespace. unneded?? */
- tokens = G_tokenize(buff, fs);
- ntokens = G_number_of_tokens(tokens);
+ x = LASPoint_GetX(LAS_point);
+ y = LASPoint_GetY(LAS_point);
+ z = LASPoint_GetZ(LAS_point);
- if ((ntokens < 3) || (max_col > ntokens)) {
- if (skipline) {
- G_warning(_("Not enough data columns. "
- "Incorrect delimiter or column number? "
- "Found the following character(s) in row %lu:\n[%s]"),
- line, buff);
- G_warning(_("Line ignored as requested"));
- continue; /* line is garbage */
- }
- else {
- G_fatal_error(_("Not enough data columns. "
- "Incorrect delimiter or column number? "
- "Found the following character(s) in row %lu:\n[%s]"),
- line, buff);
- }
- }
-
- /* too slow?
- if ( G_projection() == PROJECTION_LL ) {
- G_scan_easting( tokens[xcol-1], &x, region.proj);
- G_scan_northing( tokens[ycol-1], &y, region.proj);
- }
- else {
- */
- if (1 != sscanf(tokens[xcol - 1], "%lf", &x))
- G_fatal_error(_("Bad x-coordinate line %lu column %d. <%s>"), line,
- xcol, tokens[xcol - 1]);
-
if (first) {
min_x = x;
max_x = x;
@@ -1135,10 +1150,6 @@
max_x = x;
}
- if (1 != sscanf(tokens[ycol - 1], "%lf", &y))
- G_fatal_error(_("Bad y-coordinate line %lu column %d. <%s>"), line,
- ycol, tokens[ycol - 1]);
-
if (first) {
min_y = y;
max_y = y;
@@ -1150,10 +1161,6 @@
max_y = y;
}
- if (1 != sscanf(tokens[zcol - 1], "%lf", &z))
- G_fatal_error(_("Bad z-coordinate line %lu column %d. <%s>"), line,
- zcol, tokens[zcol - 1]);
-
if (first) {
min_z = z;
max_z = z;
@@ -1165,9 +1172,6 @@
if (z > max_z)
max_z = z;
}
-
-
- G_free_tokens(tokens);
}
if (!shell_style) {
@@ -1180,7 +1184,7 @@
fprintf(stdout, "n=%f s=%f e=%f w=%f b=%f t=%f\n",
max_y, min_y, max_x, min_x, min_z * zscale, max_z * zscale);
- G_debug(1, "Processed %lu lines.", line);
+ G_debug(1, "Processed %lu points.", line);
G_debug(1, "region template: g.region n=%f s=%f e=%f w=%f",
max_y, min_y, max_x, min_x);
Copied: grass/trunk/raster/r.in.lidar/r.in.lidar.html (from rev 46482, grass/trunk/raster/r.in.xyz/r.in.xyz.html)
===================================================================
--- grass/trunk/raster/r.in.lidar/r.in.lidar.html (rev 0)
+++ grass/trunk/raster/r.in.lidar/r.in.lidar.html 2011-06-05 09:07:11 UTC (rev 46570)
@@ -0,0 +1,275 @@
+<h2>DESCRIPTION</h2>
+
+The <em>r.in.lidar</em> module will load and bin LAS LiDAR point clouds
+into a new raster map. The user may choose from a variety of statistical
+methods in creating the new raster. Gridded data provided as a stream of
+x,y,z points may also be imported.
+<p>
+
+<em>r.in.lidar</em> is designed for processing massive point cloud datasets,
+for example raw LIDAR or sidescan sonar swath data. It has been tested with
+datasets as large as tens of billion of points (705GB in a single file).
+ <!-- Doug Newcomb, US Fish & Wildlife Service -->
+<p>
+
+Available statistics for populating the raster are:<br>
+<ul>
+<li>
+<table>
+<tr><td><em>n</em></td> <td>number of points in cell</td></tr>
+<tr><td><em>min</em></td> <td>minimum value of points in cell</td></tr>
+<tr><td><em>max</em></td> <td>maximum value of points in cell</td></tr>
+<tr><td><em>range</em></td> <td>range of points in cell</td></tr>
+<tr><td><em>sum</em></td> <td>sum of points in cell</td></tr>
+<tr><td><em>mean</em></td> <td>average value of points in cell</td></tr>
+<tr><td><em>stddev</em></td> <td>standard deviation of points in cell</td></tr>
+<tr><td><em>variance</em></td> <td>variance of points in cell</td></tr>
+<tr><td><em>coeff_var</em></td><td>coefficient of variance of points in cell</td></tr>
+<tr><td><em>median</em></td> <td>median value of points in cell</td></tr>
+<tr valign="baseline"><td><em>percentile</em> </td>
+ <td>p<sup><i>th</i></sup> percentile of points in cell</td></tr>
+<tr><td><em>skewness</em></td> <td>skewness of points in cell</td></tr>
+<tr><td><em>trimmean</em></td> <td>trimmed mean of points in cell</td></tr>
+</table><br>
+
+<li><em>Variance</em> and derivatives use the biased estimator (n). [subject to change]
+<li><em>Coefficient of variance</em> is given in percentage and defined as
+<tt>(stddev/mean)*100</tt>.
+</ul>
+<br>
+
+<h2>NOTES</h2>
+
+<h3>Gridded data</h3>
+
+If data is known to be on a regular grid <em>r.in.lidar</em> can reconstruct
+the map perfectly as long as some care is taken to set up the region
+correctly and that the data's native map projection is used. A typical
+method would involve determining the grid resolution either by examining
+the data's associated documentation or by studying the text file. Next scan
+the data with <em>r.in.lidar</em>'s <b>-s</b> (or <b>-g</b>) flag to find the
+input data's bounds. GRASS uses the cell-center raster convention where
+data points fall within the center of a cell, as opposed to the grid-node
+convention. Therefore you will need to grow the region out by half a cell
+in all directions beyond what the scan found in the file. After the region
+bounds and resolution are set correctly with <em>g.region</em>, run
+<em>r.in.lidar</em> using the <i>n</i> method and verify that n=1 at all places.
+<em>r.univar</em> can help. Once you are confident that the region exactly
+matches the data proceed to run <em>r.in.lidar</em> using one of the <i>mean,
+min, max</i>, or <i>median</i> methods. With n=1 throughout, the result
+should be identical regardless of which of those methods are used.
+
+
+<h3>Memory use</h3>
+
+While the <b>input</b> file can be arbitrarily large, <em>r.in.lidar</em>
+will use a large amount of system memory for large raster regions (10000x10000).
+If the module refuses to start complaining that there isn't enough memory,
+use the <b>percent</b> parameter to run the module in several passes.
+In addition using a less precise map format (<tt>CELL</tt> [integer] or
+<tt>FCELL</tt> [floating point]) will use less memory than a <tt>DCELL</tt>
+[double precision floating point] <b>output</b> map. Methods such as <em>n,
+min, max, sum</em> will also use less memory, while <em>stddev, variance,
+and coeff_var</em> will use more.
+
+The aggregate functions <em>median, percentile, skewness</em> and
+<em>trimmed mean</em> will use even more memory and may not be appropriate
+for use with arbitrarily large input files<!-- without a small value for percent= -->.
+<!-- explained: memory use for regular stats will be based solely on region size,
+ but for the aggregate fns it will also depend on the number of data points. (?) -->
+
+<p>
+
+The default map <b>type</b>=<tt>FCELL</tt> is intended as compromise between
+preserving data precision and limiting system resource consumption.
+If reading data from a <tt>stdin</tt> stream, the program can only run using
+a single pass.
+
+<h3>Setting region bounds and resolution</h3>
+
+You can use the <b>-s</b> scan flag to find the extent of the input data
+(and thus point density) before performing the full import. Use
+<em>g.region</em> to adjust the region bounds to match. The <b>-g</b> shell
+style flag prints the extent suitable as parameters for <em>g.region</em>.
+A suitable resolution can be found by dividing the number of input points
+by the area covered. e.g.
+
+<div class="code"><pre>
+ wc -l inputfile.txt
+ g.region -p
+ # points_per_cell = n_points / (rows * cols)
+
+ g.region -e
+ # UTM location:
+ # points_per_sq_m = n_points / (ns_extent * ew_extent)
+
+ # Lat/Lon location:
+ # points_per_sq_m = n_points / (ns_extent * ew_extent*cos(lat) * (1852*60)^2)
+</pre></div>
+
+<p>
+
+If you only intend to interpolate the data with <em>r.to.vect</em> and
+<em>v.surf.rst</em>, then there is little point to setting the region
+resolution so fine that you only catch one data point per cell -- you might
+as well use "<tt>v.in.ascii -zbt</tt>" directly.
+
+
+<h3>Filtering</h3>
+
+Points falling outside the current region will be skipped. This includes
+points falling <em>exactly</em> on the southern region bound.
+(to capture those adjust the region with "<tt>g.region s=s-0.000001</tt>";
+see <em>g.region</em>)
+<p>
+Blank lines and comment lines starting with the hash symbol (<tt>#</tt>)
+will be skipped.
+
+<p>
+
+The <b>zrange</b> parameter may be used for filtering the input data by
+vertical extent. Example uses might include preparing multiple raster
+sections to be combined into a 3D raster array with <em>r.to.rast3</em>, or
+for filtering outliers on relatively flat terrain.
+
+<p>
+
+In varied terrain the user may find that <em>min</em> maps make for a good
+noise filter as most LIDAR noise is from premature hits. The <em>min</em> map
+may also be useful to find the underlying topography in a forested or urban
+environment if the cells are over sampled.
+
+<p>
+
+The user can use a combination of <em>r.in.lidar</em> <b>output</b> maps to create
+custom filters. e.g. use <em>r.mapcalc</em> to create a <tt>mean-(2*stddev)</tt>
+map. [In this example the user may want to include a lower bound filter in
+<em>r.mapcalc</em> to remove highly variable points (small <em>n</em>) or run
+<em>r.neighbors</em> to smooth the stddev map before further use.]
+
+
+<h3>Interpolation into a DEM</h3>
+
+The vector engine's topological abilities introduce a finite memory overhead
+per vector point which will limit the size a vector map relative to
+available RAM. If you want more, use the <em>r.to.vect</em>
+<b>-b</b> flag to skip building topology. Without topology, however, all
+you'll be able to do with the vector map is display with <em>d.vect</em> and
+interpolate with <em>v.surf.rst</em>.
+Run <em>r.univar</em> on your raster map to check the number of non-NULL cells
+and adjust bounds and/or resolution as needed before proceeding.
+
+<p>
+
+Typical commands to create a DEM using a regularized spline fit:
+<div class="code"><pre>
+ r.univar lidar_min
+ r.to.vect -z feature=point in=lidar_min out=lidar_min_pt
+ v.surf.rst layer=0 in=lidar_min_pt elev=lidar_min.rst
+</pre></div>
+
+Typical commands to create a DEM using bsplines:
+<div class="code"><pre>
+ r.resamp.bspline in=lidar_min out=lidar_min.bspline
+</pre></div>
+<br>
+
+<h2>EXAMPLE</h2>
+
+Import the <a href="http://www.grassbook.org/data_menu2nd.phtml">Jockey's
+Ridge, NC, LIDAR dataset</a>, and process into a clean DEM:
+
+<div class="code"><pre>
+ # scan and set region bounds
+ r.in.xyz -s fs=, in=lidaratm2.txt out=test
+ g.region n=35.969493 s=35.949693 e=-75.620999 w=-75.639999
+ g.region res=0:00:00.075 -a
+ # create "n" map containing count of points per cell for checking density
+ r.in.xyz in=lidaratm2.txt out=lidar_n fs=, method=n zrange=-2,50
+ # check point density [rho = n_sum / (rows*cols)]
+ r.univar lidar_n | grep sum
+ # create "min" map (elevation filtered for premature hits)
+ r.in.xyz in=lidaratm2.txt out=lidar_min fs=, method=min zrange=-2,50
+ # zoom to area of interest
+ g.region n=35:57:56.25N s=35:57:13.575N w=75:38:23.7W e=75:37:15.675W
+ # check number of non-null cells (try and keep under a few million)
+ r.univar lidar_min | grep '^n:'
+ # convert to points
+ r.to.vect -z feature=point in=lidar_min out=lidar_min_pt
+ # interpolate using a regularized spline fit
+ v.surf.rst layer=0 in=lidar_min_pt elev=lidar_min.rst
+ # set color scale to something interesting
+ r.colors lidar_min.rst rule=bcyr -n -e
+ # prepare a 1:1:1 scaled version for NVIZ visualization (for lat/lon input)
+ r.mapcalc "lidar_min.rst_scaled = lidar_min.rst / (1852*60)"
+ r.colors lidar_min.rst_scaled rule=bcyr -n -e
+</pre></div>
+<br>
+
+
+<h2>TODO</h2>
+
+<ul>
+<li> Support for multiple map output from a single run.<br>
+ <tt>method=string[,string,...] output=name[,name,...]</tt>
+<li> Merge with r.in.xyz.
+ <!-- Bob Covill has supplied patches for MBIO interface already -->
+</ul>
+
+<h2>BUGS</h2>
+
+<ul>
+<li> <em>n</em> map sum can be ever-so-slightly more than `<tt>wc -l</tt>`
+ with e.g. <tt>percent=10</tt> or less.
+ <br>Cause unknown.
+
+<li> <em>n</em> map <tt>percent=100</tt> and <tt>percent=xx</tt> maps
+ differ slightly (point will fall above/below the segmentation line)
+ <br>Investigate with "<tt>r.mapcalc diff = bin_n.100 - bin_n.33</tt>" etc.
+ <br>Cause unknown.
+
+<li> "<tt>nan</tt>" can leak into <em>coeff_var</em> maps.
+ <br>Cause unknown. Possible work-around: "<tt>r.null setnull=nan</tt>"
+<!-- Another method: r.mapcalc 'No_nan = if(map == map, map, null() )' -->
+</ul>
+
+If you encounter any problems (or solutions!) please contact the GRASS
+Development Team.
+
+<h2>SEE ALSO</h2>
+
+<i>
+<a href="g.region.html">g.region</a><br>
+<a href="m.proj.html">m.proj</a><br>
+<a href="r.in.xyz.html">r.in.xyz</a><br>
+<a href="r.fillnulls.html">r.fillnulls</a><br>
+<a href="r.in.ascii.html">r.in.ascii</a><br>
+<a href="r.mapcalc.html">r.mapcalc</a><br>
+<a href="r.neighbors.html">r.neighbors</a><br>
+<a href="r.out.xyz.html">r.out.xyz</a><br>
+<a href="r.to.rast3.html">r.to.rast3</a><br>
+<a href="r.to.vect.html">r.to.vect</a><br>
+<a href="r.univar.html">r.univar</a><br>
+<a href="v.in.ascii.html">v.in.ascii</a><br>
+<a href="v.surf.rst.html">v.surf.rst</a><br>
+<br>
+<a href="v.lidar.correction.html">v.lidar.correction</a>,
+<a href="v.lidar.edgedetection.html">v.lidar.edgedetection</a>,
+<a href="v.lidar.growing.html">v.lidar.growing</a>,
+<a href="v.outlier.html">v.outlier</a>,
+<a href="v.surf.bspline.html">v.surf.bspline</a>
+</i>
+<p>
+<i><a href="http://www.ivarch.com/programs/pv.shtml">pv</a></i>
+ - The UNIX pipe viewer utility
+<BR><BR>
+
+
+<h2>AUTHORS</h2>
+
+Markus Metz<br>
+based on r.in.xyz by Hamish Bowman and Volker Wichmann<br>
+
+<br>
+<p>
+<i>Last changed: $Date$</i></p>
Deleted: grass/trunk/raster/r.in.lidar/r.in.xyz.html
===================================================================
--- grass/trunk/raster/r.in.xyz/r.in.xyz.html 2011-06-01 06:31:05 UTC (rev 46482)
+++ grass/trunk/raster/r.in.lidar/r.in.xyz.html 2011-06-05 09:07:11 UTC (rev 46570)
@@ -1,281 +0,0 @@
-<h2>DESCRIPTION</h2>
-
-The <em>r.in.xyz</em> module will load and bin ungridded x,y,z ASCII data
-into a new raster map. The user may choose from a variety of statistical
-methods in creating the new raster. Gridded data provided as a stream of
-x,y,z points may also be imported.
-<p>
-
-<em>r.in.xyz</em> is designed for processing massive point cloud datasets,
-for example raw LIDAR or sidescan sonar swath data. It has been tested with
-datasets as large as tens of billion of points (705GB in a single file).
- <!-- Doug Newcomb, US Fish & Wildlife Service -->
-<p>
-
-Available statistics for populating the raster are:<br>
-<ul>
-<li>
-<table>
-<tr><td><em>n</em></td> <td>number of points in cell</td></tr>
-<tr><td><em>min</em></td> <td>minimum value of points in cell</td></tr>
-<tr><td><em>max</em></td> <td>maximum value of points in cell</td></tr>
-<tr><td><em>range</em></td> <td>range of points in cell</td></tr>
-<tr><td><em>sum</em></td> <td>sum of points in cell</td></tr>
-<tr><td><em>mean</em></td> <td>average value of points in cell</td></tr>
-<tr><td><em>stddev</em></td> <td>standard deviation of points in cell</td></tr>
-<tr><td><em>variance</em></td> <td>variance of points in cell</td></tr>
-<tr><td><em>coeff_var</em></td><td>coefficient of variance of points in cell</td></tr>
-<tr><td><em>median</em></td> <td>median value of points in cell</td></tr>
-<tr valign="baseline"><td><em>percentile</em> </td>
- <td>p<sup><i>th</i></sup> percentile of points in cell</td></tr>
-<tr><td><em>skewness</em></td> <td>skewness of points in cell</td></tr>
-<tr><td><em>trimmean</em></td> <td>trimmed mean of points in cell</td></tr>
-</table><br>
-
-<li><em>Variance</em> and derivatives use the biased estimator (n). [subject to change]
-<li><em>Coefficient of variance</em> is given in percentage and defined as
-<tt>(stddev/mean)*100</tt>.
-</ul>
-<br>
-
-<h2>NOTES</h2>
-
-<h3>Gridded data</h3>
-
-If data is known to be on a regular grid <em>r.in.xyz</em> can reconstruct
-the map perfectly as long as some care is taken to set up the region
-correctly and that the data's native map projection is used. A typical
-method would involve determining the grid resolution either by examining
-the data's associated documentation or by studying the text file. Next scan
-the data with <em>r.in.xyz</em>'s <b>-s</b> (or <b>-g</b>) flag to find the
-input data's bounds. GRASS uses the cell-center raster convention where
-data points fall within the center of a cell, as opposed to the grid-node
-convention. Therefore you will need to grow the region out by half a cell
-in all directions beyond what the scan found in the file. After the region
-bounds and resolution are set correctly with <em>g.region</em>, run
-<em>r.in.xyz</em> using the <i>n</i> method and verify that n=1 at all places.
-<em>r.univar</em> can help. Once you are confident that the region exactly
-matches the data proceed to run <em>r.in.xyz</em> using one of the <i>mean,
-min, max</i>, or <i>median</i> methods. With n=1 throughout, the result
-should be identical regardless of which of those methods are used.
-
-
-<h3>Memory use</h3>
-
-While the <b>input</b> file can be arbitrarily large, <em>r.in.xyz</em>
-will use a large amount of system memory for large raster regions (10000x10000).
-If the module refuses to start complaining that there isn't enough memory,
-use the <b>percent</b> parameter to run the module in several passes.
-In addition using a less precise map format (<tt>CELL</tt> [integer] or
-<tt>FCELL</tt> [floating point]) will use less memory than a <tt>DCELL</tt>
-[double precision floating point] <b>output</b> map. Methods such as <em>n,
-min, max, sum</em> will also use less memory, while <em>stddev, variance,
-and coeff_var</em> will use more.
-
-The aggregate functions <em>median, percentile, skewness</em> and
-<em>trimmed mean</em> will use even more memory and may not be appropriate
-for use with arbitrarily large input files<!-- without a small value for percent= -->.
-<!-- explained: memory use for regular stats will be based solely on region size,
- but for the aggregate fns it will also depend on the number of data points. (?) -->
-
-<p>
-
-The default map <b>type</b>=<tt>FCELL</tt> is intended as compromise between
-preserving data precision and limiting system resource consumption.
-If reading data from a <tt>stdin</tt> stream, the program can only run using
-a single pass.
-
-<h3>Setting region bounds and resolution</h3>
-
-You can use the <b>-s</b> scan flag to find the extent of the input data
-(and thus point density) before performing the full import. Use
-<em>g.region</em> to adjust the region bounds to match. The <b>-g</b> shell
-style flag prints the extent suitable as parameters for <em>g.region</em>.
-A suitable resolution can be found by dividing the number of input points
-by the area covered. e.g.
-
-<div class="code"><pre>
- wc -l inputfile.txt
- g.region -p
- # points_per_cell = n_points / (rows * cols)
-
- g.region -e
- # UTM location:
- # points_per_sq_m = n_points / (ns_extent * ew_extent)
-
- # Lat/Lon location:
- # points_per_sq_m = n_points / (ns_extent * ew_extent*cos(lat) * (1852*60)^2)
-</pre></div>
-
-<p>
-
-If you only intend to interpolate the data with <em>r.to.vect</em> and
-<em>v.surf.rst</em>, then there is little point to setting the region
-resolution so fine that you only catch one data point per cell -- you might
-as well use "<tt>v.in.ascii -zbt</tt>" directly.
-
-
-<h3>Filtering</h3>
-
-Points falling outside the current region will be skipped. This includes
-points falling <em>exactly</em> on the southern region bound.
-(to capture those adjust the region with "<tt>g.region s=s-0.000001</tt>";
-see <em>g.region</em>)
-<p>
-Blank lines and comment lines starting with the hash symbol (<tt>#</tt>)
-will be skipped.
-
-<p>
-
-The <b>zrange</b> parameter may be used for filtering the input data by
-vertical extent. Example uses might include preparing multiple raster
-sections to be combined into a 3D raster array with <em>r.to.rast3</em>, or
-for filtering outliers on relatively flat terrain.
-
-<p>
-
-In varied terrain the user may find that <em>min</em> maps make for a good
-noise filter as most LIDAR noise is from premature hits. The <em>min</em> map
-may also be useful to find the underlying topography in a forested or urban
-environment if the cells are over sampled.
-
-<p>
-
-The user can use a combination of <em>r.in.xyz</em> <b>output</b> maps to create
-custom filters. e.g. use <em>r.mapcalc</em> to create a <tt>mean-(2*stddev)</tt>
-map. [In this example the user may want to include a lower bound filter in
-<em>r.mapcalc</em> to remove highly variable points (small <em>n</em>) or run
-<em>r.neighbors</em> to smooth the stddev map before further use.]
-
-
-<h3>Reprojection</h3>
-
-If the raster map is to be reprojected, it may be more appropriate to reproject
-the input points with <em>m.proj</em> or <em>cs2cs</em> before running
-<em>r.in.xyz</em>.
-
-<h3>Interpolation into a DEM</h3>
-
-The vector engine's topographic abilities introduce a finite memory overhead
-per vector point which will typically limit a vector map to approximately
-3 million points (~ 1750^2 cells). If you want more, use the <em>r.to.vect</em>
-<b>-b</b> flag to skip building topology. Without topology, however, all
-you'll be able to do with the vector map is display with <em>d.vect</em> and
-interpolate with <em>v.surf.rst</em>.
-Run <em>r.univar</em> on your raster map to check the number of non-NULL cells
-and adjust bounds and/or resolution as needed before proceeding.
-
-<p>
-
-Typical commands to create a DEM using a regularized spline fit:
-<div class="code"><pre>
- r.univar lidar_min
- r.to.vect -z feature=point in=lidar_min out=lidar_min_pt
- v.surf.rst layer=0 in=lidar_min_pt elev=lidar_min.rst
-</pre></div>
-<br>
-
-<h2>EXAMPLE</h2>
-
-Import the <a href="http://www.grassbook.org/data_menu2nd.phtml">Jockey's
-Ridge, NC, LIDAR dataset</a>, and process into a clean DEM:
-
-<div class="code"><pre>
- # scan and set region bounds
- r.in.xyz -s fs=, in=lidaratm2.txt out=test
- g.region n=35.969493 s=35.949693 e=-75.620999 w=-75.639999
- g.region res=0:00:00.075 -a
- # create "n" map containing count of points per cell for checking density
- r.in.xyz in=lidaratm2.txt out=lidar_n fs=, method=n zrange=-2,50
- # check point density [rho = n_sum / (rows*cols)]
- r.univar lidar_n | grep sum
- # create "min" map (elevation filtered for premature hits)
- r.in.xyz in=lidaratm2.txt out=lidar_min fs=, method=min zrange=-2,50
- # zoom to area of interest
- g.region n=35:57:56.25N s=35:57:13.575N w=75:38:23.7W e=75:37:15.675W
- # check number of non-null cells (try and keep under a few million)
- r.univar lidar_min | grep '^n:'
- # convert to points
- r.to.vect -z feature=point in=lidar_min out=lidar_min_pt
- # interpolate using a regularized spline fit
- v.surf.rst layer=0 in=lidar_min_pt elev=lidar_min.rst
- # set color scale to something interesting
- r.colors lidar_min.rst rule=bcyr -n -e
- # prepare a 1:1:1 scaled version for NVIZ visualization (for lat/lon input)
- r.mapcalc "lidar_min.rst_scaled = lidar_min.rst / (1852*60)"
- r.colors lidar_min.rst_scaled rule=bcyr -n -e
-</pre></div>
-<br>
-
-
-<h2>TODO</h2>
-
-<ul>
-<li> Support for multiple map output from a single run.<br>
- <tt>method=string[,string,...] output=name[,name,...]</tt>
-<li> Add two new flags for support for direct binary input from libLAS
- for LIDAR data and MB-System's mbio for multi-beam bathymetry data.
- <!-- Bob Covill has supplied patches for MBIO interface already -->
-</ul>
-
-<h2>BUGS</h2>
-
-<ul>
-<li> <em>n</em> map sum can be ever-so-slightly more than `<tt>wc -l</tt>`
- with e.g. <tt>percent=10</tt> or less.
- <br>Cause unknown.
-
-<li> <em>n</em> map <tt>percent=100</tt> and <tt>percent=xx</tt> maps
- differ slightly (point will fall above/below the segmentation line)
- <br>Investigate with "<tt>r.mapcalc diff = bin_n.100 - bin_n.33</tt>" etc.
- <br>Cause unknown.
-
-<li> "<tt>nan</tt>" can leak into <em>coeff_var</em> maps.
- <br>Cause unknown. Possible work-around: "<tt>r.null setnull=nan</tt>"
-<!-- Another method: r.mapcalc 'No_nan = if(map == map, map, null() )' -->
-</ul>
-
-If you encounter any problems (or solutions!) please contact the GRASS
-Development Team.
-
-<h2>SEE ALSO</h2>
-
-<i>
-<a href="g.region.html">g.region</a><br>
-<a href="m.proj.html">m.proj</a><br>
-<a href="r.fillnulls.html">r.fillnulls</a><br>
-<a href="r.in.ascii.html">r.in.ascii</a><br>
-<a href="r.mapcalc.html">r.mapcalc</a><br>
-<a href="r.neighbors.html">r.neighbors</a><br>
-<a href="r.out.xyz.html">r.out.xyz</a><br>
-<a href="r.to.rast3.html">r.to.rast3</a><br>
-<a href="r.to.vect.html">r.to.vect</a><br>
-<a href="r.univar.html">r.univar</a><br>
-<a href="v.in.ascii.html">v.in.ascii</a><br>
-<a href="v.surf.rst.html">v.surf.rst</a><br>
-<br>
-<a href="v.lidar.correction.html">v.lidar.correction</a>,
-<a href="v.lidar.edgedetection.html">v.lidar.edgedetection</a>,
-<a href="v.lidar.growing.html">v.lidar.growing</a>,
-<a href="v.outlier.html">v.outlier</a>,
-<a href="v.surf.bspline.html">v.surf.bspline</a>
-</i>
-<p>
-<i><a href="http://www.ivarch.com/programs/pv.shtml">pv</a></i>
- - The UNIX pipe viewer utility
-<BR><BR>
-
-
-<h2>AUTHORS</h2>
-
-Hamish Bowman<br> <i>
-Department of Marine Science<br>
-University of Otago<br>
-New Zealand</i><br>
-<br>
-Extended by Volker Wichmann to support the aggregate functions
-<i>median, percentile, skewness</i> and <i>trimmed mean</i>.
-
-<br>
-<p>
-<i>Last changed: $Date$</i></p>
Modified: grass/trunk/raster/r.in.lidar/support.c
===================================================================
--- grass/trunk/raster/r.in.xyz/support.c 2011-06-01 06:31:05 UTC (rev 46482)
+++ grass/trunk/raster/r.in.lidar/support.c 2011-06-05 09:07:11 UTC (rev 46570)
@@ -1,5 +1,5 @@
/*
- * r.in.xyz support fns.
+ * r.in.lidar support fns, from r.in.xyz.
* Copyright 2006 by M. Hamish Bowman, and The GRASS Development Team
* Author: M. Hamish Bowman, University of Otago, Dunedin, New Zealand
*
More information about the grass-commit
mailing list