[GRASS-SVN] r46418 - in grass/trunk/vector: . v.in.lidar
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed May 25 05:16:36 EDT 2011
Author: mmetz
Date: 2011-05-25 02:16:36 -0700 (Wed, 25 May 2011)
New Revision: 46418
Added:
grass/trunk/vector/v.in.lidar/
Removed:
grass/trunk/vector/v.in.lidar/geom.c
grass/trunk/vector/v.in.lidar/global.h
Modified:
grass/trunk/vector/v.in.lidar/Makefile
grass/trunk/vector/v.in.lidar/main.c
Log:
new module to import LiDAR LAS format
Modified: grass/trunk/vector/v.in.lidar/Makefile
===================================================================
--- grass/trunk/vector/v.in.ogr/Makefile 2011-05-22 11:51:32 UTC (rev 46356)
+++ grass/trunk/vector/v.in.lidar/Makefile 2011-05-25 09:16:36 UTC (rev 46418)
@@ -1,16 +1,16 @@
MODULE_TOPDIR = ../..
-PGM=v.in.ogr
+PGM=v.in.lidar
-LIBES = $(GPROJLIB) $(VECTORLIB) $(DBMILIB) $(GISLIB) $(GDALLIBS) $(MATHLIB)
+LIBES = $(GPROJLIB) $(VECTORLIB) $(DBMILIB) $(GISLIB) $(LASLIBS) $(MATHLIB)
DEPENDENCIES = $(GPROJDEP) $(VECTORDEP) $(DBMIDEP) $(GISDEP)
-EXTRA_INC = $(VECT_INC) $(PROJINC)
+EXTRA_INC = $(VECT_INC) $(PROJINC) $(LASINC)
EXTRA_CFLAGS = $(VECT_CFLAGS)
include $(MODULE_TOPDIR)/include/Make/Module.make
-ifneq ($(USE_OGR),)
+ifneq ($(USE_LIBLAS),)
default: cmd
endif
Deleted: grass/trunk/vector/v.in.lidar/geom.c
===================================================================
--- grass/trunk/vector/v.in.ogr/geom.c 2011-05-22 11:51:32 UTC (rev 46356)
+++ grass/trunk/vector/v.in.lidar/geom.c 2011-05-25 09:16:36 UTC (rev 46418)
@@ -1,477 +0,0 @@
-
-/****************************************************************
- *
- * MODULE: v.in.ogr
- *
- * AUTHOR(S): Radim Blazek
- *
- * PURPOSE: Import OGR vectors
- *
- * COPYRIGHT: (C) 2001 by the GRASS Development Team
- *
- * This program is free software under the
- * GNU General Public License (>=v2).
- * Read the file COPYING that comes with GRASS
- * for details.
- *
- **************************************************************/
-
-#include <stdlib.h>
-#include <string.h>
-#include <grass/gis.h>
-#include <grass/dbmi.h>
-#include <grass/vector.h>
-#include <grass/glocale.h>
-#include "ogr_api.h"
-#include "global.h"
-
-int split_line(struct Map_info *Map, int otype, struct line_pnts *Points,
- struct line_cats *Cats);
-
-/* Add categories to centroids inside polygon */
-int
-centroid(OGRGeometryH hGeom, CENTR * Centr, struct spatial_index *Sindex,
- int field, int cat, double min_area, int type)
-{
- int i, valid_isles, j, np, nr, ret;
- static int first = 1;
- static struct line_pnts *Points;
- struct line_pnts **IPoints;
- static struct line_cats *BCats, *Cats;
- OGRwkbGeometryType eType;
- OGRGeometryH hRing;
- double size;
- static struct ilist *List;
- struct bound_box box;
-
- G_debug(3, "centroid() cat = %d", cat);
-
- if (first) {
- Points = Vect_new_line_struct();
- BCats = Vect_new_cats_struct();
- Cats = Vect_new_cats_struct();
- List = Vect_new_list();
- first = 0;
- }
- else {
- Vect_reset_line(Points);
- Vect_reset_cats(Cats);
- Vect_reset_cats(BCats);
- Vect_cat_set(Cats, field, cat);
- }
-
- eType = wkbFlatten(OGR_G_GetGeometryType(hGeom));
-
- if (eType == wkbPolygon) {
- nr = OGR_G_GetGeometryCount(hGeom);
-
- G_debug(3, "polygon: %d rings", nr);
-
- /* SFS: 1 exterior boundary and 0 or more interior boundaries.
- * So I hope that exterior is the first one, even if it is not explicitly told */
-
- /* Area */
- hRing = OGR_G_GetGeometryRef(hGeom, 0);
- np = OGR_G_GetPointCount(hRing);
- Vect_reset_line(Points);
- for (j = 0; j < np; j++) {
- Vect_append_point(Points, OGR_G_GetX(hRing, j),
- OGR_G_GetY(hRing, j), OGR_G_GetZ(hRing, j));
- }
-
- /* Degenerate is ignored */
- if (Points->n_points < 4)
- return 0;
-
- /* Small areas ignored because boundaries are not imported */
- size = G_area_of_polygon(Points->x, Points->y, Points->n_points);
- if (size < min_area)
- return 0;
-
- /* Isles */
- IPoints =
- (struct line_pnts **)G_malloc((nr - 1) *
- sizeof(struct line_pnts *));
- valid_isles = 0;
- for (i = 1; i < nr; i++) {
-
- hRing = OGR_G_GetGeometryRef(hGeom, i);
- if ((np = OGR_G_GetPointCount(hRing)) > 0) {
- IPoints[valid_isles] = Vect_new_line_struct();
- for (j = 0; j < np; j++) {
- Vect_append_point(IPoints[valid_isles],
- OGR_G_GetX(hRing, j),
- OGR_G_GetY(hRing, j),
- OGR_G_GetZ(hRing, j));
- }
- size =
- G_area_of_polygon(IPoints[valid_isles]->x,
- IPoints[valid_isles]->y,
- IPoints[valid_isles]->n_points);
- if (size < min_area)
- Vect_destroy_line_struct(IPoints[valid_isles]);
- else
- valid_isles++;
- }
- }
-
- /* Find centroids */
- if (Points->n_points >= 4) {
- int centr, in;
- double x, y;
-
- Vect_line_box(Points, &box);
- Vect_spatial_index_select(Sindex, &box, List);
-
- for (i = 0; i < List->n_values; i++) {
- centr = List->value[i];
- x = Centr[centr].x;
- y = Centr[centr].y;
- ret = Vect_point_in_poly(x, y, Points);
- if (ret == 0)
- continue; /* outside */
-
- in = 1;
- for (j = 0; j < valid_isles; j++) {
- ret = Vect_point_in_poly(x, y, IPoints[j]);
- if (ret == 1) { /* centroid in inner ring */
- in = 0;
- break; /* inside isle */
- }
- }
- if (!in)
- continue;
-
- G_debug(3, "Centroid %d : layer %d cat %d", centr, field,
- cat);
- Vect_cat_set(Centr[centr].cats, field, cat);
- }
- }
-
- for (i = 0; i < valid_isles; i++) {
- Vect_destroy_line_struct(IPoints[i]);
- }
- G_free(IPoints);
- }
-
- /* I did not test this because I did not have files of these types */
- else if (eType == wkbGeometryCollection || eType == wkbMultiPolygon) {
- G_debug(3, "GeometryCollection or MultiPolygon/LineString/Point");
- nr = OGR_G_GetGeometryCount(hGeom);
- for (i = 0; i < nr; i++) {
- hRing = OGR_G_GetGeometryRef(hGeom, i);
-
- ret = centroid(hRing, Centr, Sindex, field, cat, min_area, type);
- }
- }
-
- return 0;
-}
-
-/* count polygons and isles */
-int poly_count(OGRGeometryH hGeom)
-{
- int i, nr, ret;
- OGRwkbGeometryType eType;
- OGRGeometryH hRing;
- eType = wkbFlatten(OGR_G_GetGeometryType(hGeom));
-
- if (eType == wkbPolygon) {
- G_debug(3, "Polygon");
- nr = OGR_G_GetGeometryCount(hGeom);
- n_polygon_boundaries += nr;
-
- }
- else if (eType == wkbGeometryCollection || eType == wkbMultiPolygon) {
- G_debug(3, "GeometryCollection or MultiPolygon");
- nr = OGR_G_GetGeometryCount(hGeom);
- for (i = 0; i < nr; i++) {
- hRing = OGR_G_GetGeometryRef(hGeom, i);
-
- ret = poly_count(hRing);
- if (ret == -1) {
- G_warning(_("Cannot read part of geometry"));
- }
- }
- }
- return 0;
-}
-
-/* Write geometry to output map */
-int
-geom(OGRGeometryH hGeom, struct Map_info *Map, int field, int cat,
- double min_area, int type, int mk_centr)
-{
- int i, valid_isles, j, np, nr, ret, otype;
- static int first = 1;
- static struct line_pnts *Points;
- struct line_pnts **IPoints;
- static struct line_cats *BCats, *Cats;
- OGRwkbGeometryType eType;
- OGRGeometryH hRing;
- double x, y;
- double size;
-
- G_debug(3, "geom() cat = %d", cat);
-
- if (first) {
- Points = Vect_new_line_struct();
- BCats = Vect_new_cats_struct();
- Cats = Vect_new_cats_struct();
- first = 0;
- }
- Vect_reset_line(Points);
- Vect_reset_cats(Cats);
- Vect_reset_cats(BCats);
- Vect_cat_set(Cats, field, cat);
-
- eType = wkbFlatten(OGR_G_GetGeometryType(hGeom));
-
- if (eType == wkbPoint) {
- if ((np = OGR_G_GetPointCount(hGeom)) == 0) {
- G_warning(_("Skipping empty geometry feature"));
- return 0;
- }
- Vect_append_point(Points, OGR_G_GetX(hGeom, 0), OGR_G_GetY(hGeom, 0),
- OGR_G_GetZ(hGeom, 0));
- if (type & GV_CENTROID)
- otype = GV_CENTROID;
- else
- otype = GV_POINT;
- Vect_write_line(Map, otype, Points, Cats);
- }
- else if (eType == wkbLineString) {
- if ((np = OGR_G_GetPointCount(hGeom)) == 0) {
- G_warning(_("Skipping empty geometry feature"));
- return 0;
- }
-
- for (i = 0; i < np; i++) {
- Vect_append_point(Points, OGR_G_GetX(hGeom, i),
- OGR_G_GetY(hGeom, i), OGR_G_GetZ(hGeom, i));
- }
- if (type & GV_BOUNDARY)
- otype = GV_BOUNDARY;
- else
- otype = GV_LINE;
-
- if (split_distance > 0 && otype == GV_BOUNDARY)
- split_line(Map, otype, Points, Cats);
- else
- Vect_write_line(Map, otype, Points, Cats);
- }
-
- else if (eType == wkbPolygon) {
- G_debug(3, "Polygon");
-
- /* SFS: 1 exterior boundary and 0 or more interior boundaries.
- * So I hope that exterior is the first one, even if it is not explicitly told */
-
- /* Area */
- hRing = OGR_G_GetGeometryRef(hGeom, 0);
- if ((np = OGR_G_GetPointCount(hRing)) == 0) {
- G_warning(_("Skipping empty geometry feature"));
- return 0;
- }
-
- n_polygons++;
- nr = OGR_G_GetGeometryCount(hGeom);
-
- Vect_reset_line(Points);
- for (j = 0; j < np; j++) {
- Vect_append_point(Points, OGR_G_GetX(hRing, j),
- OGR_G_GetY(hRing, j), OGR_G_GetZ(hRing, j));
- }
-
- /* Degenerate is not ignored because it may be useful to see where it is,
- * but may be eliminated by min_area option */
- if (Points->n_points < 4)
- G_warning(_("Degenerate polygon ([%d] vertices)"),
- Points->n_points);
-
- size = G_area_of_polygon(Points->x, Points->y, Points->n_points);
- if (size < min_area) {
- G_debug(2, "Area size [%.1e], area not imported", size);
- return 0;
- }
-
- if (type & GV_LINE)
- otype = GV_LINE;
- else
- otype = GV_BOUNDARY;
-
- if (split_distance > 0 && otype == GV_BOUNDARY)
- split_line(Map, otype, Points, BCats);
- else
- Vect_write_line(Map, otype, Points, BCats);
-
- /* Isles */
- IPoints =
- (struct line_pnts **)G_malloc((nr - 1) *
- sizeof(struct line_pnts *));
- valid_isles = 0;
- for (i = 1; i < nr; i++) {
- G_debug(3, "Inner ring %d", i);
-
- hRing = OGR_G_GetGeometryRef(hGeom, i);
-
- if ((np = OGR_G_GetPointCount(hRing)) == 0) {
- G_warning(_("Skipping empty geometry feature"));
- }
- else {
- IPoints[valid_isles] = Vect_new_line_struct();
-
- for (j = 0; j < np; j++) {
- Vect_append_point(IPoints[valid_isles],
- OGR_G_GetX(hRing, j),
- OGR_G_GetY(hRing, j),
- OGR_G_GetZ(hRing, j));
- }
-
- if (IPoints[valid_isles]->n_points < 4)
- G_warning(_("Degenerate island ([%d] vertices)"),
- IPoints[i - 1]->n_points);
-
- size =
- G_area_of_polygon(IPoints[valid_isles]->x,
- IPoints[valid_isles]->y,
- IPoints[valid_isles]->n_points);
- if (size < min_area) {
- G_debug(2, "Island size [%.1e], island not imported",
- size);
- }
- else {
- if (type & GV_LINE)
- otype = GV_LINE;
- else
- otype = GV_BOUNDARY;
- if (split_distance > 0 && otype == GV_BOUNDARY)
- split_line(Map, otype, IPoints[valid_isles], BCats);
- else
- Vect_write_line(Map, otype, IPoints[valid_isles],
- BCats);
- }
- valid_isles++;
- }
- } /* inner rings done */
-
- /* Centroid */
- /* Vect_get_point_in_poly_isl() would fail for degenerate polygon */
- if (mk_centr) {
- if (Points->n_points >= 4) {
- ret =
- Vect_get_point_in_poly_isl(Points, IPoints, valid_isles,
- &x, &y);
- if (ret == -1) {
- G_warning(_("Cannot calculate centroid"));
- }
- else {
- Vect_reset_line(Points);
- Vect_append_point(Points, x, y, 0.0);
- if (type & GV_POINT)
- otype = GV_POINT;
- else
- otype = GV_CENTROID;
- Vect_write_line(Map, otype, Points, Cats);
- }
- }
- else if (Points->n_points > 0) {
- if (Points->n_points >= 2) {
- /* center of 1. segment ( 2. point is not best for 3 vertices as 3. may be the same as 1.) */
- x = (Points->x[0] + Points->x[1]) / 2;
- y = (Points->y[0] + Points->y[1]) / 2;
- }
- else { /* one point */
- x = Points->x[0];
- y = Points->y[0];
- }
- Vect_reset_line(Points);
- Vect_append_point(Points, x, y, 0.0);
- if (type & GV_POINT)
- otype = GV_POINT;
- else
- otype = GV_CENTROID;
- Vect_write_line(Map, otype, Points, Cats);
- }
- else { /* 0 points */
- G_warning(_("No centroid written for polygon with 0 vertices"));
- }
- }
-
- for (i = 0; i < valid_isles; i++) {
- Vect_destroy_line_struct(IPoints[i]);
- }
- G_free(IPoints);
- }
-
- /* I did not test this because I did not have files of these types */
- else if (eType == wkbGeometryCollection
- || eType == wkbMultiPolygon
- || eType == wkbMultiLineString || eType == wkbMultiPoint) {
- G_debug(3, "GeometryCollection or MultiPolygon/LineString/Point");
- nr = OGR_G_GetGeometryCount(hGeom);
- for (i = 0; i < nr; i++) {
- hRing = OGR_G_GetGeometryRef(hGeom, i);
-
- ret = geom(hRing, Map, field, cat, min_area, type, mk_centr);
- if (ret == -1) {
- G_warning(_("Cannot write part of geometry"));
- }
- }
- }
-
- else {
- G_fatal_error(_("Unknown geometry type"));
- }
-
- return 0;
-}
-
-int split_line(struct Map_info *Map, int otype, struct line_pnts *Points,
- struct line_cats *Cats)
-{
- int i;
- double dist = 0., seg_dist, dx, dy;
- struct line_pnts *OutPoints;
-
- /* don't write zero length boundaries */
- Vect_line_prune(Points);
- if (Points->n_points < 2)
- return 0;
-
- /* can't split boundaries with only 2 vertices */
- if (Points->n_points == 2) {
- Vect_write_line(Map, otype, Points, Cats);
- return 0;
- }
-
- OutPoints = Vect_new_line_struct();
- Vect_append_point(OutPoints, Points->x[0], Points->y[0], Points->z[0]);
- Vect_append_point(OutPoints, Points->x[1], Points->y[1], Points->z[1]);
- dx = Points->x[1] - Points->x[0];
- dy = Points->y[1] - Points->y[0];
- dist = sqrt(dx * dx + dy * dy);
-
- /* trying to keep line length smaller than split_distance
- * alternative, less code: write line as soon as split_distance is exceeded */
- for (i = 2; i < Points->n_points; i++) {
- dx = Points->x[i] - Points->x[i - 1];
- dy = Points->y[i] - Points->y[i - 1];
- seg_dist = sqrt(dx * dx + dy * dy);
- dist += seg_dist;
- if (dist > split_distance) {
- Vect_write_line(Map, otype, OutPoints, Cats);
- Vect_reset_line(OutPoints);
- dist = seg_dist;
- Vect_append_point(OutPoints, Points->x[i - 1], Points->y[i - 1],
- Points->z[i - 1]);
- }
- Vect_append_point(OutPoints, Points->x[i], Points->y[i],
- Points->z[i]);
- }
- Vect_write_line(Map, otype, OutPoints, Cats);
-
- Vect_destroy_line_struct(OutPoints);
-
- return 0;
-}
Deleted: grass/trunk/vector/v.in.lidar/global.h
===================================================================
--- grass/trunk/vector/v.in.ogr/global.h 2011-05-22 11:51:32 UTC (rev 46356)
+++ grass/trunk/vector/v.in.lidar/global.h 2011-05-25 09:16:36 UTC (rev 46418)
@@ -1,40 +0,0 @@
-
-/****************************************************************
- *
- * MODULE: v.in.ogr
- *
- * AUTHOR(S): Radim Blazek
- * Markus Neteler (spatial parm, projection support)
- * Paul Kelly (projection support)
- *
- * PURPOSE: Import OGR vectors
- *
- * COPYRIGHT: (C) 2003 by the GRASS Development Team
- *
- * This program is free software under the
- * GNU General Public License (>=v2).
- * Read the file COPYING that comes with GRASS
- * for details.
- *
- * TODO: - make fixed field length of OFTIntegerList dynamic
- * - several other TODOs below
-**************************************************************/
-
-#ifndef __GLOBAL_H__
-#define __GLOBAL_H__
-
-
-extern int n_polygons;
-extern int n_polygon_boundaries;
-extern double split_distance;
-
-/* centroid structure */
-typedef struct
-{
- double x, y;
- struct line_cats *cats;
- int valid;
-} CENTR;
-
-
-#endif /* __GLOBAL_H__ */
Modified: grass/trunk/vector/v.in.lidar/main.c
===================================================================
--- grass/trunk/vector/v.in.ogr/main.c 2011-05-22 11:51:32 UTC (rev 46356)
+++ grass/trunk/vector/v.in.lidar/main.c 2011-05-25 09:16:36 UTC (rev 46418)
@@ -1,13 +1,12 @@
/****************************************************************
*
- * MODULE: v.in.ogr
+ * MODULE: v.in.lidar
*
- * AUTHOR(S): Radim Blazek
- * Markus Neteler (spatial parm, projection support)
- * Paul Kelly (projection support)
+ * AUTHOR(S): Markus Metz
+ * based on v.in.ogr
*
- * PURPOSE: Import OGR vectors
+ * PURPOSE: Import LiDAR LAS points
*
* COPYRIGHT: (C) 2003 by the GRASS Development Team
*
@@ -27,108 +26,134 @@
#include <grass/vector.h>
#include <grass/gprojects.h>
#include <grass/glocale.h>
-#include <gdal_version.h> /* needed for OFTDate */
-#include "ogr_api.h"
-#include "global.h"
+#include <liblas/capi/liblas.h>
#ifndef MAX
# define MIN(a,b) ((a<b) ? a : b)
# define MAX(a,b) ((a>b) ? a : b)
#endif
-int n_polygons;
-int n_polygon_boundaries;
-double split_distance;
+/*
+ * ASPRS Standard LIDAR Point Classes
+ * Classification Value (bits 0:4) : Meaning
+ * 0 : Created, never classified
+ * 1 : Unclassified
+ * 2 : Ground
+ * 3 : Low Vegetation
+ * 4 : Medium Vegetation
+ * 5 : High Vegetation
+ * 6 : Building
+ * 7 : Low Point (noise)
+ * 8 : Model Key-point (mass point)
+ * 9 : Water
+ * 10 : Reserved for ASPRS Definition
+ * 11 : Reserved for ASPRS Definition
+ * 12 : Overlap Points
+ * 13-31 : Reserved for ASPRS Definition
+ */
-int geom(OGRGeometryH hGeom, struct Map_info *Map, int field, int cat,
- double min_area, int type, int mk_centr);
-int centroid(OGRGeometryH hGeom, CENTR * Centr, struct spatial_index * Sindex,
- int field, int cat, double min_area, int type);
-int poly_count(OGRGeometryH hGeom);
+/* Classification Bit Field Encoding
+ * Bits | Field Name | Description
+ * 0-4 | Classification | Standard ASPRS classification as defined in the
+ * above classification table.
+ * 5 | Synthetic | If set then this point was created by a technique
+ * other than LIDAR collection such as digitized from
+ * a photogrammetric stereo model or by traversing
+ * a waveform.
+ * 6 | Key-point | If set, this point is considered to be a model
+ * key-point and thus generally should not be withheld
+ * in a thinning algorithm.
+ * 7 | Withheld | If set, this point should not be included in
+ * processing (synonymous with Deleted).
+*/
+struct class_table
+{
+ int code;
+ char *name;
+};
+
+static struct class_table class_val[] = {
+ {0, "Created, never classified"},
+ {1, "Unclassified"},
+ {2, "Ground"},
+ {3, "Low Vegetation"},
+ {4, "Medium Vegetation"},
+ {5, "High Vegetation"},
+ {6, "Building"},
+ {7, "Low Point (noise)"},
+ {8, "Model Key-point (mass point)"},
+ {9, "Water"},
+ {10, "Reserved for ASPRS Definition"},
+ {11, "Reserved for ASPRS Definition"},
+ {12, "Overlap Points"},
+ {13 /* 13 - 31 */, "Reserved for ASPRS Definition"},
+ {0, 0}
+};
+
+static struct class_table class_type[] = {
+ {5, "Synthetic"},
+ {6, "Key-point"},
+ {7, "withheld"},
+ {0, 0}
+};
+
+void print_lasinfo(LASHeaderH LAS_header, LASSRSH LAS_srs);
+
int main(int argc, char *argv[])
{
- int i, j, layer, arg_s_num, nogeom, ncnames;
+ int i;
float xmin = 0., ymin = 0., xmax = 0., ymax = 0.;
- int ncols = 0, type;
struct GModule *module;
- double min_area, snap;
- struct Option *dsn_opt, *out_opt, *layer_opt, *spat_opt, *where_opt,
- *min_area_opt;
- struct Option *snap_opt, *type_opt, *outloc_opt, *cnames_opt;
- struct Flag *list_flag, *no_clean_flag, *z_flag, *notab_flag,
- *region_flag;
- struct Flag *over_flag, *extend_flag, *formats_flag, *tolower_flag, *no_import_flag;
- char buf[2000], namebuf[2000], tempvect[2000];
- char *separator;
+ struct Option *in_opt, *out_opt, *spat_opt;
+ struct Option *outloc_opt;
+ struct Flag *print_flag, *notab_flag, *region_flag, *notopo_flag;
+ struct Flag *over_flag, *extend_flag, *no_import_flag;
+ char buf[2000];
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, cur_wind;
char error_msg[8192];
/* Vector */
- struct Map_info Map, Tmp;
+ struct Map_info Map;
int cat;
/* Attributes */
struct field_info *Fi;
dbDriver *driver;
dbString sql, strval;
- int dim, with_z;
+
+ /* LAS */
+ LASReaderH LAS_reader;
+ LASHeaderH LAS_header;
+ LASSRSH LAS_srs;
+ LASPointH LAS_point;
+ double scale_x, scale_y, scale_z, offset_x, offset_y, offset_z;
+ int las_point_format;
+ int have_time, have_color;
+ unsigned int not_valid;
- /* OGR */
- OGRDataSourceH Ogr_ds = NULL;
- OGRLayerH Ogr_layer;
- OGRFieldDefnH Ogr_field;
- char *Ogr_fieldname;
- OGRFieldType Ogr_ftype;
- OGRFeatureH Ogr_feature;
- OGRFeatureDefnH Ogr_featuredefn;
- OGRGeometryH Ogr_geometry, Ogr_oRing = NULL, poSpatialFilter = NULL;
- OGRSpatialReferenceH Ogr_projection;
- OGREnvelope oExt;
- int OFTIntegerListlength = 40; /* hack due to limitation in OGR */
+ struct line_pnts *Points;
+ struct line_cats *Cats;
- char **layer_names; /* names of layers to be imported */
- int *layers; /* layer indexes */
- int nlayers; /* number of layers to import */
- char **available_layer_names; /* names of layers to be imported */
- int navailable_layers;
- int layer_id;
- unsigned int n_features, feature_count;
+ unsigned int n_features, feature_count, n_outside;
int overwrite;
- double area_size = 0.;
G_gisinit(argv[0]);
module = G_define_module();
G_add_keyword(_("vector"));
G_add_keyword(_("import"));
- module->description = _("Converts vector daya into a GRASS vector map using OGR library.");
+ module->description = _("Converts LAS LiDAR point clouds to a GRASS vector map with libLAS.");
- dsn_opt = G_define_option();
- dsn_opt->key = "dsn";
- dsn_opt->type = TYPE_STRING;
- dsn_opt->required =YES;
- dsn_opt->label = _("OGR datasource name");
- dsn_opt->description = _("Examples:\n"
- "\t\tESRI Shapefile: directory containing shapefiles\n"
- "\t\tMapInfo File: directory containing mapinfo files");
+ in_opt = G_define_standard_option(G_OPT_F_INPUT);
+ in_opt->label = _("LAS input file");
+ in_opt->description = _("LiDAR input files in LAS format (*.las or *.laz)");
out_opt = G_define_standard_option(G_OPT_V_OUTPUT);
- layer_opt = G_define_option();
- layer_opt->key = "layer";
- layer_opt->type = TYPE_STRING;
- layer_opt->required = NO;
- layer_opt->multiple = YES;
- layer_opt->label =
- _("OGR layer name. If not given, all available layers are imported");
- layer_opt->description =
- _("Examples:\n" "\t\tESRI Shapefile: shapefile name\n"
- "\t\tMapInfo File: mapinfo file name");
- layer_opt->guisection = _("Selection");
-
spat_opt = G_define_option();
spat_opt->key = "spatial";
spat_opt->type = TYPE_DOUBLE;
@@ -140,40 +165,6 @@
spat_opt->description =
_("Format: xmin,ymin,xmax,ymax - usually W,S,E,N");
- where_opt = G_define_standard_option(G_OPT_DB_WHERE);
- where_opt->guisection = _("Selection");
-
- min_area_opt = G_define_option();
- min_area_opt->key = "min_area";
- min_area_opt->type = TYPE_DOUBLE;
- min_area_opt->required = NO;
- min_area_opt->answer = "0.0001";
- min_area_opt->label =
- _("Minimum size of area to be imported (square units)");
- min_area_opt->guisection = _("Min-area & snap");
- min_area_opt->description = _("Smaller areas and "
- "islands are ignored. Should be greater than snap^2");
-
- type_opt = G_define_standard_option(G_OPT_V_TYPE);
- type_opt->options = "point,line,boundary,centroid";
- type_opt->answer = "";
- type_opt->description = _("Optionally change default input type");
- type_opt->descriptions =
- _("point;import area centroids as points;"
- "line;import area boundaries as lines;"
- "boundary;import lines as area boundaries;"
- "centroid;import points as centroids");
- type_opt->guisection = _("Selection");
-
- snap_opt = G_define_option();
- snap_opt->key = "snap";
- snap_opt->type = TYPE_DOUBLE;
- snap_opt->required = NO;
- snap_opt->answer = "-1";
- snap_opt->label = _("Snapping threshold for boundaries");
- snap_opt->guisection = _("Min-area & snap");
- snap_opt->description = _("'-1' for no snap");
-
outloc_opt = G_define_option();
outloc_opt->key = "location";
outloc_opt->type = TYPE_STRING;
@@ -181,36 +172,12 @@
outloc_opt->description = _("Name for new location to create");
outloc_opt->key_desc = "name";
- cnames_opt = G_define_option();
- cnames_opt->key = "cnames";
- cnames_opt->type = TYPE_STRING;
- cnames_opt->required = NO;
- cnames_opt->multiple = YES;
- cnames_opt->description =
- _("List of column names to be used instead of original names, "
- "first is used for category column");
- cnames_opt->guisection = _("Attributes");
-
- list_flag = G_define_flag();
- list_flag->key = 'l';
- list_flag->description =
- _("List available layers in data source and exit");
- list_flag->suppress_required = YES;
+ print_flag = G_define_flag();
+ print_flag->key = 'p';
+ print_flag->description =
+ _("Print LAS file info and exit");
+ print_flag->suppress_required = YES;
- formats_flag = G_define_flag();
- formats_flag->key = 'f';
- formats_flag->description = _("List supported formats and exit");
- formats_flag->suppress_required = YES;
-
- /* if using -c, you lose topological information ! */
- no_clean_flag = G_define_flag();
- no_clean_flag->key = 'c';
- no_clean_flag->description = _("Do not clean polygons (not recommended)");
-
- z_flag = G_define_flag();
- z_flag->key = 'z';
- z_flag->description = _("Create 3D output");
-
notab_flag = G_define_flag();
notab_flag->key = 't';
notab_flag->description = _("Do not create attribute table");
@@ -231,18 +198,16 @@
extend_flag->description =
_("Extend location extents based on new dataset");
- tolower_flag = G_define_flag();
- tolower_flag->key = 'w';
- tolower_flag->description =
- _("Change column names to lowercase characters");
- tolower_flag->guisection = _("Attributes");
-
no_import_flag = G_define_flag();
no_import_flag->key = 'i';
no_import_flag->description =
_("Create the location specified by the \"location\" parameter and exit."
" Do not import the vector file.");
+ notopo_flag = G_define_flag();
+ notopo_flag->key = 'b';
+ notopo_flag->description = _("Do not build topology");
+
/* The parser checks if the map already exists in current mapset, this is
* wrong if location options is used, so we switch out the check and do it
* in the module after the parser */
@@ -251,77 +216,33 @@
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
- G_begin_polygon_area_calculations(); /* Used in geom() */
+ /* Open LAS file*/
+ LAS_reader = LASReader_Create(in_opt->answer);
+ LAS_header = LASReader_GetHeader(LAS_reader);
+ LAS_srs = LASHeader_GetSRS(LAS_header);
- OGRRegisterAll();
+ scale_x = LASHeader_GetScaleX(LAS_header);
+ scale_y = LASHeader_GetScaleY(LAS_header);
+ scale_z = LASHeader_GetScaleZ(LAS_header);
- /* list supported formats */
- if (formats_flag->answer) {
- int iDriver;
+ offset_x = LASHeader_GetOffsetX(LAS_header);
+ offset_y = LASHeader_GetOffsetY(LAS_header);
+ offset_z = LASHeader_GetOffsetZ(LAS_header);
- G_important_message(_("Available OGR Drivers:"));
+ xmin = LASHeader_GetMinX(LAS_header) * scale_x + offset_x;
+ xmax = LASHeader_GetMaxX(LAS_header) * scale_x + offset_x;
+ ymin = LASHeader_GetMinY(LAS_header) * scale_y + offset_y;
+ ymax = LASHeader_GetMaxY(LAS_header) * scale_y + offset_y;
- for (iDriver = 0; iDriver < OGRGetDriverCount(); iDriver++) {
- OGRSFDriverH poDriver = OGRGetDriver(iDriver);
- const char *pszRWFlag;
-
- if (OGR_Dr_TestCapability(poDriver, ODrCCreateDataSource))
- pszRWFlag = "rw";
- else
- pszRWFlag = "ro";
+ /* Print LAS header */
+ if (print_flag->answer) {
+ /* print... */
+ print_lasinfo(LAS_header, LAS_srs);
- fprintf(stdout, " %s (%s): %s\n",
- OGR_Dr_GetName(poDriver),
- pszRWFlag, OGR_Dr_GetName(poDriver));
- }
- exit(EXIT_SUCCESS);
- }
+ LASSRS_Destroy(LAS_srs);
+ LASHeader_Destroy(LAS_header);
+ LASReader_Destroy(LAS_reader);
- if (dsn_opt->answer == NULL) {
- G_fatal_error(_("Required parameter <%s> not set"), dsn_opt->key);
- }
-
- min_area = atof(min_area_opt->answer);
- snap = atof(snap_opt->answer);
- type = Vect_option_to_types(type_opt);
-
- ncnames = 0;
- if (cnames_opt->answers) {
- i = 0;
- while (cnames_opt->answers[i++]) {
- ncnames++;
- }
- }
-
- /* Open OGR DSN */
- Ogr_ds = NULL;
- if (strlen(dsn_opt->answer) > 0)
- Ogr_ds = OGROpen(dsn_opt->answer, FALSE, NULL);
-
- if (Ogr_ds == NULL)
- G_fatal_error(_("Unable to open data source <%s>"), dsn_opt->answer);
-
- /* Make a list of available layers */
- navailable_layers = OGR_DS_GetLayerCount(Ogr_ds);
- available_layer_names =
- (char **)G_malloc(navailable_layers * sizeof(char *));
-
- if (list_flag->answer)
- G_important_message(_("Data source contains %d layers:"),
- navailable_layers);
-
- for (i = 0; i < navailable_layers; i++) {
- Ogr_layer = OGR_DS_GetLayer(Ogr_ds, i);
- Ogr_featuredefn = OGR_L_GetLayerDefn(Ogr_layer);
- available_layer_names[i] =
- G_store((char *)OGR_FD_GetName(Ogr_featuredefn));
-
- if (list_flag->answer) {
- fprintf(stdout, "%s\n", available_layer_names[i]);
- }
- }
- if (list_flag->answer) {
- fflush(stdout);
exit(EXIT_SUCCESS);
}
@@ -336,40 +257,7 @@
}
}
- /* Make a list of layers to be imported */
- if (layer_opt->answer) { /* From option */
- nlayers = 0;
- while (layer_opt->answers[nlayers])
- nlayers++;
- layer_names = (char **)G_malloc(nlayers * sizeof(char *));
- layers = (int *)G_malloc(nlayers * sizeof(int));
-
- for (i = 0; i < nlayers; i++) {
- layer_names[i] = G_store(layer_opt->answers[i]);
- /* Find it in the source */
- layers[i] = -1;
- for (j = 0; j < navailable_layers; j++) {
- if (strcmp(available_layer_names[j], layer_names[i]) == 0) {
- layers[i] = j;
- break;
- }
- }
- if (layers[i] == -1)
- G_fatal_error(_("Layer <%s> not available"), layer_names[i]);
- }
- }
- else { /* use list of all layers */
- nlayers = navailable_layers;
- layer_names = available_layer_names;
- layers = (int *)G_malloc(nlayers * sizeof(int));
- for (i = 0; i < nlayers; i++)
- layers[i] = i;
- }
-
- /* Get first imported layer to use for extents and projection check */
- Ogr_layer = OGR_DS_GetLayer(Ogr_ds, layers[0]);
-
if (region_flag->answer) {
if (spat_opt->answer)
G_fatal_error(_("Select either the current region flag or the spatial option, not both"));
@@ -385,7 +273,7 @@
/* cut out a piece of the map */
/* order: xmin,ymin,xmax,ymax */
- arg_s_num = 0;
+ int arg_s_num = 0;
i = 0;
while (spat_opt->answers[i]) {
if (i == 0)
@@ -405,78 +293,30 @@
if (spat_opt->answer || region_flag->answer) {
G_debug(2, "cut out with boundaries: xmin:%f ymin:%f xmax:%f ymax:%f",
xmin, ymin, xmax, ymax);
-
- /* in theory this could be an irregular polygon */
- poSpatialFilter = OGR_G_CreateGeometry(wkbPolygon);
- Ogr_oRing = OGR_G_CreateGeometry(wkbLinearRing);
- OGR_G_AddPoint(Ogr_oRing, xmin, ymin, 0.0);
- OGR_G_AddPoint(Ogr_oRing, xmin, ymax, 0.0);
- OGR_G_AddPoint(Ogr_oRing, xmax, ymax, 0.0);
- OGR_G_AddPoint(Ogr_oRing, xmax, ymin, 0.0);
- OGR_G_AddPoint(Ogr_oRing, xmin, ymin, 0.0);
- OGR_G_AddGeometryDirectly(poSpatialFilter, Ogr_oRing);
-
- OGR_L_SetSpatialFilter(Ogr_layer, poSpatialFilter);
}
- if (where_opt->answer) {
-
- /* select by attribute */
- OGR_L_SetAttributeFilter(Ogr_layer, where_opt->answer);
- }
-
/* fetch boundaries */
- if ((OGR_L_GetExtent(Ogr_layer, &oExt, 1)) == OGRERR_NONE) {
- G_get_window(&cellhd);
- cellhd.north = oExt.MaxY;
- cellhd.south = oExt.MinY;
- cellhd.west = oExt.MinX;
- cellhd.east = oExt.MaxX;
- cellhd.rows = 20; /* TODO - calculate useful values */
- cellhd.cols = 20;
- cellhd.ns_res = (cellhd.north - cellhd.south) / cellhd.rows;
- cellhd.ew_res = (cellhd.east - cellhd.west) / cellhd.cols;
- }
- else {
- cellhd.north = 1.;
- cellhd.south = 0.;
- cellhd.west = 0.;
- cellhd.east = 1.;
- cellhd.top = 1.;
- cellhd.bottom = 1.;
- cellhd.rows = 1;
- cellhd.rows3 = 1;
- cellhd.cols = 1;
- cellhd.cols3 = 1;
- cellhd.depths = 1;
- cellhd.ns_res = 1.;
- cellhd.ns_res3 = 1.;
- cellhd.ew_res = 1.;
- cellhd.ew_res3 = 1.;
- cellhd.tb_res = 1.;
- }
+ G_get_window(&cellhd);
+ cellhd.north = ymax;
+ cellhd.south = ymin;
+ cellhd.west = xmin;
+ cellhd.east = xmax;
+ cellhd.rows = 20; /* TODO - calculate useful values */
+ cellhd.cols = 20;
+ cellhd.ns_res = (cellhd.north - cellhd.south) / cellhd.rows;
+ cellhd.ew_res = (cellhd.east - cellhd.west) / cellhd.cols;
- /* suppress boundary splitting ? */
- if (no_clean_flag->answer) {
- split_distance = -1.;
- }
- else {
- split_distance = 0.;
- area_size =
- sqrt((cellhd.east - cellhd.west) * (cellhd.north - cellhd.south));
- }
-
/* Fetch input map projection in GRASS form. */
proj_info = NULL;
proj_units = NULL;
- Ogr_projection = OGR_L_GetSpatialRef(Ogr_layer); /* should not be freed later */
+ projstr = LASSRS_GetWKT_CompoundOK(LAS_srs);
/* Do we need to create a new location? */
if (outloc_opt->answer != NULL) {
/* Convert projection information non-interactively as we can't
* assume the user has a terminal open */
- if (GPJ_osr_to_grass(&cellhd, &proj_info,
- &proj_units, Ogr_projection, 0) < 0) {
+ if (GPJ_wkt_to_grass(&cellhd, &proj_info,
+ &proj_units, projstr, 0) < 0) {
G_fatal_error(_("Unable to convert input map projection to GRASS "
"format; cannot create new location."));
}
@@ -496,14 +336,17 @@
int err = 0;
/* Projection only required for checking so convert non-interactively */
- if (GPJ_osr_to_grass(&cellhd, &proj_info,
- &proj_units, Ogr_projection, 0) < 0)
+ 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_window(&loc_wind, "", "DEFAULT_WIND", "PERMANENT");
+ if (outloc_opt->answer)
+ G__get_window(&loc_wind, "", "DEFAULT_WIND", "PERMANENT");
+ else
+ G_get_set_window(&loc_wind);
/* fetch LOCATION PROJ info */
if (loc_wind.proj != PROJECTION_XY) {
loc_proj_info = G_get_projinfo();
@@ -605,587 +448,286 @@
db_init_string(&strval);
/* open output vector */
- /* open temporary vector, do the work in the temporary vector
- * at the end copy alive lines to output vector
- * in case of polygons this reduces the coor file size by a factor of 2 to 5
- * only needed for polygons, but the presence of polygons can be detected
- * only during OGR feature import, not before */
sprintf(buf, "%s", out_opt->answer);
/* strip any @mapset from vector output name */
G_find_vector(buf, G_mapset());
- sprintf(tempvect, "%s_tmp", buf);
- G_verbose_message(_("Using temporary vector <%s>"), tempvect);
- Vect_open_new(&Map, out_opt->answer, z_flag->answer != 0);
- Vect_open_new(&Tmp, tempvect, z_flag->answer != 0);
+ Vect_open_new(&Map, out_opt->answer, 1);
Vect_hist_command(&Map);
+
+ n_features = LASHeader_GetPointRecordsCount(LAS_header);
+ las_point_format = LASHeader_GetDataFormatId(LAS_header);
- /* Points and lines are written immediately with categories. Boundaries of polygons are
- * written to the vector then cleaned and centroids are calculated for all areas in cleaan vector.
- * Then second pass through finds all centroids in each polygon feature and adds its category
- * to the centroid. The result is that one centroids may have 0, 1 ore more categories
- * of one ore more (more input layers) fields. */
- with_z = 0;
- for (layer = 0; layer < nlayers; layer++) {
- G_message(_("Layer: %s"), layer_names[layer]);
- layer_id = layers[layer];
+ have_time = (las_point_format == 1 || las_point_format == 3 ||
+ las_point_format == 4 || las_point_format == 5);
- Ogr_layer = OGR_DS_GetLayer(Ogr_ds, layer_id);
- Ogr_featuredefn = OGR_L_GetLayerDefn(Ogr_layer);
+ have_color = (las_point_format == 2 || las_point_format == 3 ||
+ las_point_format == 5);
- /* Add DB link */
- if (!notab_flag->answer) {
- char *cat_col_name = "cat";
+ /* Add DB link */
+ if (!notab_flag->answer) {
+ char *cat_col_name = "cat";
- if (nlayers == 1) { /* one layer only */
- Fi = Vect_default_field_info(&Map, layer + 1, NULL,
- GV_1TABLE);
- }
- else {
- Fi = Vect_default_field_info(&Map, layer + 1, NULL,
- GV_MTABLE);
- }
+ Fi = Vect_default_field_info(&Map, 1, NULL, GV_1TABLE);
- if (ncnames > 0) {
- cat_col_name = cnames_opt->answers[0];
- }
- Vect_map_add_dblink(&Map, layer + 1, layer_names[layer], Fi->table,
- cat_col_name, Fi->database, Fi->driver);
+ Vect_map_add_dblink(&Map, 1, out_opt->answer, Fi->table,
+ cat_col_name, Fi->database, Fi->driver);
- ncols = OGR_FD_GetFieldCount(Ogr_featuredefn);
- G_debug(2, "%d columns", ncols);
-
- /* Create table */
- sprintf(buf, "create table %s (%s integer", Fi->table,
- cat_col_name);
- db_set_string(&sql, buf);
- for (i = 0; i < ncols; i++) {
-
- Ogr_field = OGR_FD_GetFieldDefn(Ogr_featuredefn, i);
- Ogr_ftype = OGR_Fld_GetType(Ogr_field);
-
- G_debug(3, "Ogr_ftype: %i", Ogr_ftype); /* look up below */
-
- if (i < ncnames - 1) {
- Ogr_fieldname = G_store(cnames_opt->answers[i + 1]);
- }
- else {
- /* Change column names to [A-Za-z][A-Za-z0-9_]* */
- Ogr_fieldname = G_store(OGR_Fld_GetNameRef(Ogr_field));
- G_debug(3, "Ogr_fieldname: '%s'", Ogr_fieldname);
-
- G_str_to_sql(Ogr_fieldname);
-
- G_debug(3, "Ogr_fieldname: '%s'", Ogr_fieldname);
-
- }
-
- /* avoid that we get the 'cat' column twice */
- if (strcmp(Ogr_fieldname, "cat") == 0) {
- sprintf(namebuf, "%s_", Ogr_fieldname);
- Ogr_fieldname = G_store(namebuf);
- }
-
- /* captial column names are a pain in SQL */
- if (tolower_flag->answer)
- G_str_to_lower(Ogr_fieldname);
-
- if (strcmp(OGR_Fld_GetNameRef(Ogr_field), Ogr_fieldname) != 0) {
- G_warning(_("Column name changed: '%s' -> '%s'"),
- OGR_Fld_GetNameRef(Ogr_field), Ogr_fieldname);
- }
-
- /** Simple 32bit integer OFTInteger = 0 **/
-
- /** List of 32bit integers OFTIntegerList = 1 **/
-
- /** Double Precision floating point OFTReal = 2 **/
-
- /** List of doubles OFTRealList = 3 **/
-
- /** String of ASCII chars OFTString = 4 **/
-
- /** Array of strings OFTStringList = 5 **/
-
- /** Double byte string (unsupported) OFTWideString = 6 **/
-
- /** List of wide strings (unsupported) OFTWideStringList = 7 **/
-
- /** Raw Binary data (unsupported) OFTBinary = 8 **/
-
- /** OFTDate = 9 **/
-
- /** OFTTime = 10 **/
-
- /** OFTDateTime = 11 **/
-
-
- if (Ogr_ftype == OFTInteger) {
- sprintf(buf, ", %s integer", Ogr_fieldname);
- }
- else if (Ogr_ftype == OFTIntegerList) {
- /* hack: treat as string */
- sprintf(buf, ", %s varchar ( %d )", Ogr_fieldname,
- OFTIntegerListlength);
- G_warning(_("Writing column <%s> with fixed length %d chars (may be truncated)"),
- Ogr_fieldname, OFTIntegerListlength);
- }
- else if (Ogr_ftype == OFTReal) {
- sprintf(buf, ", %s double precision", Ogr_fieldname);
-#if GDAL_VERSION_NUM >= 1320
- }
- else if (Ogr_ftype == OFTDate) {
- sprintf(buf, ", %s date", Ogr_fieldname);
- }
- else if (Ogr_ftype == OFTTime) {
- sprintf(buf, ", %s time", Ogr_fieldname);
- }
- else if (Ogr_ftype == OFTDateTime) {
- sprintf(buf, ", %s datetime", Ogr_fieldname);
-#endif
- }
- else if (Ogr_ftype == OFTString) {
- int fwidth;
-
- fwidth = OGR_Fld_GetWidth(Ogr_field);
- /* TODO: read all records first and find the longest string length */
- if (fwidth == 0) {
- G_warning(_("Width for column %s set to 255 (was not specified by OGR), "
- "some strings may be truncated!"),
- Ogr_fieldname);
- fwidth = 255;
- }
- sprintf(buf, ", %s varchar ( %d )", Ogr_fieldname,
- fwidth);
- }
- else if (Ogr_ftype == OFTStringList) {
- /* hack: treat as string */
- sprintf(buf, ", %s varchar ( %d )", Ogr_fieldname,
- OFTIntegerListlength);
- G_warning(_("Writing column %s with fixed length %d chars (may be truncated)"),
- Ogr_fieldname, OFTIntegerListlength);
- }
- else {
- G_warning(_("Column type not supported (%s)"),
- Ogr_fieldname);
- buf[0] = 0;
- }
- db_append_string(&sql, buf);
- G_free(Ogr_fieldname);
- }
- db_append_string(&sql, ")");
- G_debug(3, db_get_string(&sql));
-
- driver =
- db_start_driver_open_database(Fi->driver,
- Vect_subst_var(Fi->database,
- &Map));
- if (driver == NULL) {
- G_fatal_error(_("Unable open database <%s> by driver <%s>"),
- Vect_subst_var(Fi->database, &Map), Fi->driver);
- }
-
- if (db_execute_immediate(driver, &sql) != DB_OK) {
- db_close_database(driver);
- db_shutdown_driver(driver);
- G_fatal_error(_("Unable to create table: '%s'"),
- db_get_string(&sql));
- }
-
- if (db_create_index2(driver, Fi->table, cat_col_name) != DB_OK)
- G_warning(_("Unable to create index for table <%s>, key <%s>"),
- Fi->table, cat_col_name);
-
- if (db_grant_on_table
- (driver, Fi->table, DB_PRIV_SELECT,
- DB_GROUP | DB_PUBLIC) != DB_OK)
- G_fatal_error(_("Unable to grant privileges on table <%s>"),
- Fi->table);
-
- db_begin_transaction(driver);
+ /* check available LAS info, depends on POINT DATA RECORD FORMAT [0-5] */
+ /* X (double),
+ * Y (double),
+ * Z (double),
+ * intensity (double),
+ * return number (int),
+ * number of returns (int),
+ * scan direction (int),
+ * flight line edge (int),
+ * classification type (char),
+ * class (char),
+ * time (double) (FORMAT 1, 3, 4, 5),
+ * scan angle rank (int),
+ * source ID (int),
+ * user data (char), ???
+ * red (int) (FORMAT 2, 3, 5),
+ * green (int) (FORMAT 2, 3, 5),
+ * blue (int) (FORMAT 2, 3, 5)*/
+
+ /* Create table */
+ sprintf(buf, "create table %s (%s integer", Fi->table,
+ cat_col_name);
+ db_set_string(&sql, buf);
+
+ /* x, y, z */
+ sprintf(buf, ", x_coord double precision");
+ db_append_string(&sql, buf);
+ sprintf(buf, ", y_coord double precision");
+ db_append_string(&sql, buf);
+ sprintf(buf, ", z_coord double precision");
+ db_append_string(&sql, buf);
+ /* intensity */
+ sprintf(buf, ", intensity integer");
+ db_append_string(&sql, buf);
+ /* return number */
+ sprintf(buf, ", return integer");
+ db_append_string(&sql, buf);
+ /* number of returns */
+ sprintf(buf, ", n_returns integer");
+ db_append_string(&sql, buf);
+ /* scan direction */
+ sprintf(buf, ", scan_dir integer");
+ db_append_string(&sql, buf);
+ /* flight line edge */
+ sprintf(buf, ", edge integer");
+ db_append_string(&sql, buf);
+ /* classification type */
+ sprintf(buf, ", cl_type varchar(20)");
+ db_append_string(&sql, buf);
+ /* classification class */
+ sprintf(buf, ", class varchar(40)");
+ db_append_string(&sql, buf);
+ /* GPS time */
+ if (have_time) {
+ sprintf(buf, ", gps_time double precision");
+ db_append_string(&sql, buf);
}
+ /* scan angle */
+ sprintf(buf, ", angle integer");
+ db_append_string(&sql, buf);
+ /* source id */
+ sprintf(buf, ", src_id integer");
+ db_append_string(&sql, buf);
+ /* user data */
+ sprintf(buf, ", usr_data integer");
+ db_append_string(&sql, buf);
+ /* colors */
+ if (have_color) {
+ sprintf(buf, ", red integer, green integer, blue integer");
+ db_append_string(&sql, buf);
+ sprintf(buf, ", GRASSRGB varchar(11)");
+ db_append_string(&sql, buf);
+ }
- /* Import feature */
- cat = 1;
- nogeom = 0;
- OGR_L_ResetReading(Ogr_layer);
- n_features = feature_count = 0;
+ db_append_string(&sql, ")");
+ G_debug(3, db_get_string(&sql));
- n_polygon_boundaries = 0;
- n_features = OGR_L_GetFeatureCount(Ogr_layer, 1);
-
- /* estimate distance for boundary splitting --> */
-
- if (split_distance > -0.5) {
- /* count polygons and isles */
- G_message(_("Counting polygons for %d features..."), n_features);
- while ((Ogr_feature = OGR_L_GetNextFeature(Ogr_layer)) != NULL) {
- G_percent(feature_count++, n_features, 1); /* show something happens */
- /* Geometry */
- Ogr_geometry = OGR_F_GetGeometryRef(Ogr_feature);
- if (Ogr_geometry != NULL) {
- poly_count(Ogr_geometry);
- }
- OGR_F_Destroy(Ogr_feature);
- }
- /* rewind layer */
- OGR_L_ResetReading(Ogr_layer);
- feature_count = 0;
+ driver =
+ db_start_driver_open_database(Fi->driver,
+ Vect_subst_var(Fi->database,
+ &Map));
+ if (driver == NULL) {
+ G_fatal_error(_("Unable open database <%s> by driver <%s>"),
+ Vect_subst_var(Fi->database, &Map), Fi->driver);
}
- G_debug(1, "n polygon boundaries: %d", n_polygon_boundaries);
- if (split_distance > -0.5 && n_polygon_boundaries > 50) {
- split_distance =
- area_size / log(n_polygon_boundaries);
- /* divisor is the handle: increase divisor to decrease split_distance */
- split_distance = split_distance / 4.;
- G_debug(1, "root of area size: %f", area_size);
- G_verbose_message(_("Boundary splitting distance in map units: %G"),
- split_distance);
+ if (db_execute_immediate(driver, &sql) != DB_OK) {
+ db_close_database(driver);
+ db_shutdown_driver(driver);
+ G_fatal_error(_("Unable to create table: '%s'"),
+ db_get_string(&sql));
}
- /* <-- estimate distance for boundary splitting */
-
- G_important_message(_("Importing map %d features..."), n_features);
- while ((Ogr_feature = OGR_L_GetNextFeature(Ogr_layer)) != NULL) {
- G_percent(feature_count++, n_features, 1); /* show something happens */
- /* Geometry */
- Ogr_geometry = OGR_F_GetGeometryRef(Ogr_feature);
- if (Ogr_geometry == NULL) {
- nogeom++;
- }
- else {
- dim = OGR_G_GetCoordinateDimension(Ogr_geometry);
- if (dim > 2)
- with_z = 1;
- geom(Ogr_geometry, &Tmp, layer + 1, cat, min_area, type,
- no_clean_flag->answer);
- }
+ if (db_create_index2(driver, Fi->table, cat_col_name) != DB_OK)
+ G_warning(_("Unable to create index for table <%s>, key <%s>"),
+ Fi->table, cat_col_name);
- /* Attributes */
- if (!notab_flag->answer) {
- sprintf(buf, "insert into %s values ( %d", Fi->table, cat);
- db_set_string(&sql, buf);
- for (i = 0; i < ncols; i++) {
- Ogr_field = OGR_FD_GetFieldDefn(Ogr_featuredefn, i);
- Ogr_ftype = OGR_Fld_GetType(Ogr_field);
- if (OGR_F_IsFieldSet(Ogr_feature, i)) {
- if (Ogr_ftype == OFTInteger || Ogr_ftype == OFTReal) {
- sprintf(buf, ", %s",
- OGR_F_GetFieldAsString(Ogr_feature, i));
-#if GDAL_VERSION_NUM >= 1320
- /* should we use OGR_F_GetFieldAsDateTime() here ? */
- }
- else if (Ogr_ftype == OFTDate || Ogr_ftype == OFTTime
- || Ogr_ftype == OFTDateTime) {
- char *newbuf;
+ if (db_grant_on_table
+ (driver, Fi->table, DB_PRIV_SELECT,
+ DB_GROUP | DB_PUBLIC) != DB_OK)
+ G_fatal_error(_("Unable to grant privileges on table <%s>"),
+ Fi->table);
- db_set_string(&strval, (char *)
- OGR_F_GetFieldAsString(Ogr_feature,
- i));
- db_double_quote_string(&strval);
- sprintf(buf, ", '%s'", db_get_string(&strval));
- newbuf = G_str_replace(buf, "/", "-"); /* fix 2001/10/21 to 2001-10-21 */
- sprintf(buf, "%s", newbuf);
-#endif
- }
- else if (Ogr_ftype == OFTString ||
- Ogr_ftype == OFTIntegerList) {
- db_set_string(&strval, (char *)
- OGR_F_GetFieldAsString(Ogr_feature,
- i));
- db_double_quote_string(&strval);
- sprintf(buf, ", '%s'", db_get_string(&strval));
- }
-
- }
- else {
- /* G_warning (_("Column value not set" )); */
- if (Ogr_ftype == OFTInteger || Ogr_ftype == OFTReal) {
- sprintf(buf, ", NULL");
-#if GDAL_VERSION_NUM >= 1320
- }
- else if (Ogr_ftype == OFTString ||
- Ogr_ftype == OFTIntegerList ||
- Ogr_ftype == OFTDate) {
-#else
- }
- else if (Ogr_ftype == OFTString ||
- Ogr_ftype == OFTIntegerList) {
-#endif
- sprintf(buf, ", ''");
- }
- }
- db_append_string(&sql, buf);
- }
- db_append_string(&sql, " )");
- G_debug(3, db_get_string(&sql));
-
- if (db_execute_immediate(driver, &sql) != DB_OK) {
- db_close_database(driver);
- db_shutdown_driver(driver);
- G_fatal_error(_("Cannot insert new row: %s"),
- db_get_string(&sql));
- }
- }
-
- OGR_F_Destroy(Ogr_feature);
- cat++;
- }
- G_percent(n_features, n_features, 1); /* finish it */
-
- if (!notab_flag->answer) {
- db_commit_transaction(driver);
- db_close_database_shutdown_driver(driver);
- }
-
- if (nogeom > 0)
- G_warning(_("%d %s without geometry"), nogeom,
- nogeom == 1 ? "feature" : "features");
+ db_begin_transaction(driver);
}
+ /* Import feature */
+ cat = 1;
+ not_valid = 0;
+ feature_count = 0;
+ n_outside = 0;
- separator = "-----------------------------------------------------";
- G_message("%s", separator);
+ Points = Vect_new_line_struct();
+ Cats = Vect_new_cats_struct();
+
+ G_important_message(_("Importing %d points..."), n_features);
+ while ((LAS_point = LASReader_GetNextPoint(LAS_reader)) != NULL) {
+ double x, y, z;
- /* TODO: is it necessary to build here? probably not, consumes time */
- /* GV_BUILD_BASE is sufficient to toggle boundary cleaning */
- Vect_build_partial(&Tmp, GV_BUILD_BASE);
-
- if (!no_clean_flag->answer &&
- Vect_get_num_primitives(&Tmp, GV_BOUNDARY) > 0) {
- int ret, centr, ncentr, otype, n_overlaps, n_nocat;
- CENTR *Centr;
- struct spatial_index si;
- double x, y, total_area, overlap_area, nocat_area;
- struct bound_box box;
- struct line_pnts *Points;
- int nmodif;
-
- Points = Vect_new_line_struct();
-
- G_message("%s", separator);
-
- G_warning(_("Cleaning polygons, result is not guaranteed!"));
-
- if (snap >= 0) {
- G_message("%s", separator);
- G_message(_("Snap boundaries (threshold = %.3e):"), snap);
- Vect_snap_lines(&Tmp, GV_BOUNDARY, snap, NULL);
+ G_percent(feature_count++, n_features, 1); /* show something happens */
+
+ if (!LASPoint_IsValid(LAS_point)) {
+ not_valid++;
+ continue;
}
- /* It is not to clean to snap centroids, but I have seen data with 2 duplicate polygons
- * (as far as decimal places were printed) and centroids were not identical */
- /* Disabled, because overlapping polygons result in many duplicate centroids anyway */
- /*
- fprintf ( stderr, separator );
- fprintf ( stderr, "Snap centroids (threshold 0.000001):\n" );
- Vect_snap_lines ( &Map, GV_CENTROID, 0.000001, NULL, stderr );
- */
-
- G_message("%s", separator);
- G_message(_("Break polygons:"));
- Vect_break_polygons(&Tmp, GV_BOUNDARY, NULL);
-
- /* It is important to remove also duplicate centroids in case of duplicate input polygons */
- G_message("%s", separator);
- G_message(_("Remove duplicates:"));
- Vect_remove_duplicates(&Tmp, GV_BOUNDARY | GV_CENTROID, NULL);
-
- /* in non-pathological cases, the bulk of the cleaning is now done */
-
- /* Vect_clean_small_angles_at_nodes() can change the geometry so that new intersections
- * are created. We must call Vect_break_lines(), Vect_remove_duplicates()
- * and Vect_clean_small_angles_at_nodes() until no more small angles are found */
- do {
- G_message("%s", separator);
- G_message(_("Break boundaries:"));
- Vect_break_lines(&Tmp, GV_BOUNDARY, NULL);
-
- G_message("%s", separator);
- G_message(_("Remove duplicates:"));
- Vect_remove_duplicates(&Tmp, GV_BOUNDARY, NULL);
-
- G_message("%s", separator);
- G_message(_("Clean boundaries at nodes:"));
- nmodif =
- Vect_clean_small_angles_at_nodes(&Tmp, GV_BOUNDARY, NULL);
- } while (nmodif > 0);
-
- G_message("%s", separator);
- if (type & GV_BOUNDARY) { /* that means lines were converted boundaries */
- G_message(_("Change boundary dangles to lines:"));
- Vect_chtype_dangles(&Tmp, -1.0, NULL);
- }
- else {
- G_message(_("Remove dangles:"));
- Vect_remove_dangles(&Tmp, GV_BOUNDARY, -1.0, NULL);
- }
-
- G_message("%s", separator);
- if (type & GV_BOUNDARY) {
- G_message(_("Change boundary bridges to lines:"));
- Vect_chtype_bridges(&Tmp, NULL);
- }
- else {
- G_message(_("Remove bridges:"));
- Vect_remove_bridges(&Tmp, NULL);
- }
-
- /* merge boundaries */
- G_message("%s", separator);
- G_message(_("Merge boundaries:"));
- Vect_merge_lines(&Tmp, GV_BOUNDARY, NULL, NULL);
-
- /* Boundaries are hopefully clean, build areas */
- G_message("%s", separator);
- Vect_build_partial(&Tmp, GV_BUILD_ATTACH_ISLES);
-
- /* Calculate new centroids for all areas, centroids have the same id as area */
- ncentr = Vect_get_num_areas(&Tmp);
- G_debug(3, "%d centroids/areas", ncentr);
-
- Centr = (CENTR *) G_calloc(ncentr + 1, sizeof(CENTR));
- Vect_spatial_index_init(&si, 0);
- for (centr = 1; centr <= ncentr; centr++) {
- Centr[centr].valid = 0;
- Centr[centr].cats = Vect_new_cats_struct();
- ret = Vect_get_point_in_area(&Tmp, centr, &x, &y);
- if (ret < 0) {
- G_warning(_("Cannot calculate area centroid"));
+ Vect_reset_line(Points);
+ Vect_reset_cats(Cats);
+
+ x = LASPoint_GetX(LAS_point) * scale_x + offset_x;
+ y = LASPoint_GetY(LAS_point) * scale_y + offset_y;
+ z = LASPoint_GetZ(LAS_point) * scale_z + offset_z;
+
+ if (spat_opt->answer || region_flag->answer) {
+ if (x < xmin || x > xmax || y < ymin || y > ymax) {
+ n_outside++;
continue;
}
-
- Centr[centr].x = x;
- Centr[centr].y = y;
- Centr[centr].valid = 1;
- box.N = box.S = y;
- box.E = box.W = x;
- box.T = box.B = 0;
- Vect_spatial_index_add_item(&si, centr, &box);
}
- /* Go through all layers and find centroids for each polygon */
- for (layer = 0; layer < nlayers; layer++) {
- G_message("%s", separator);
- G_message(_("Find centroids for layer: %s"), layer_names[layer]);
- layer_id = layers[layer];
- Ogr_layer = OGR_DS_GetLayer(Ogr_ds, layer_id);
- n_features = OGR_L_GetFeatureCount(Ogr_layer, 1);
- OGR_L_ResetReading(Ogr_layer);
+ Vect_append_point(Points, x, y, z);
+ Vect_cat_set(Cats, 1, cat);
+ Vect_write_line(&Map, GV_POINT, Points, Cats);
- cat = 0; /* field = layer + 1 */
- G_percent(cat, n_features, 2);
- while ((Ogr_feature = OGR_L_GetNextFeature(Ogr_layer)) != NULL) {
- cat++;
- G_percent(cat, n_features, 2);
- /* Geometry */
- Ogr_geometry = OGR_F_GetGeometryRef(Ogr_feature);
- if (Ogr_geometry != NULL) {
- centroid(Ogr_geometry, Centr, &si, layer + 1, cat,
- min_area, type);
- }
+ /* Attributes */
+ if (!notab_flag->answer) {
+ char class_flag;
+ int las_class_type, las_class;
- OGR_F_Destroy(Ogr_feature);
- }
- }
+ /* use LASPoint_Validate (LASPointH hPoint) to check for
+ * return number, number of returns, scan direction, flight line edge,
+ * classification, scan angle rank */
+ sprintf(buf, "insert into %s values ( %d", Fi->table, cat);
+ db_set_string(&sql, buf);
- /* Write centroids */
- G_message("%s", separator);
- G_message(_("Write centroids:"));
-
- n_overlaps = n_nocat = 0;
- total_area = overlap_area = nocat_area = 0.0;
- for (centr = 1; centr <= ncentr; centr++) {
- double area;
-
- G_percent(centr, ncentr, 2);
-
- area = Vect_get_area_area(&Tmp, centr);
- total_area += area;
-
- if (!(Centr[centr].valid)) {
- continue;
+ /* x, y, z */
+ sprintf(buf, ", %f", x);
+ db_append_string(&sql, buf);
+ sprintf(buf, ", %f", y);
+ db_append_string(&sql, buf);
+ sprintf(buf, ", %f", z);
+ db_append_string(&sql, buf);
+ /* intensity */
+ sprintf(buf, ", %d", LASPoint_GetIntensity(LAS_point));
+ db_append_string(&sql, buf);
+ /* return number */
+ sprintf(buf, ", %d", LASPoint_GetReturnNumber(LAS_point));
+ db_append_string(&sql, buf);
+ /* number of returns */
+ sprintf(buf, ", %d", LASPoint_GetNumberOfReturns(LAS_point));
+ db_append_string(&sql, buf);
+ /* scan direction */
+ sprintf(buf, ", %d", LASPoint_GetScanDirection(LAS_point));
+ db_append_string(&sql, buf);
+ /* flight line edge */
+ sprintf(buf, ", %d", LASPoint_GetFlightLineEdge(LAS_point));
+ db_append_string(&sql, buf);
+ class_flag = LASPoint_GetClassification(LAS_point);
+ /* classification type int or char ? */
+ las_class_type = class_flag / 32;
+ sprintf(buf, ", \"%s\"", class_type[las_class_type].name);
+ db_append_string(&sql, buf);
+ /* classification class int or char ? */
+ las_class = class_flag % 32;
+ if (las_class > 13)
+ las_class = 13;
+ sprintf(buf, ", \"%s\"", class_val[las_class].name);
+ db_append_string(&sql, buf);
+ /* GPS time */
+ if (have_time) {
+ sprintf(buf, ", %f", LASPoint_GetTime(LAS_point));
+ db_append_string(&sql, buf);
}
+ /* scan angle */
+ sprintf(buf, ", %d", LASPoint_GetScanAngleRank(LAS_point));
+ db_append_string(&sql, buf);
+ /* source id */
+ sprintf(buf, ", %d", LASPoint_GetPointSourceId(LAS_point));
+ db_append_string(&sql, buf);
+ /* user data */
+ sprintf(buf, ", %d", LASPoint_GetUserData(LAS_point));
+ db_append_string(&sql, buf);
+ /* colors */
+ if (have_color) {
+ LASColorH LAS_color = LASPoint_GetColor(LAS_point);
+ int red = LASColor_GetRed(LAS_color);
+ int green = LASColor_GetGreen(LAS_color);
+ int blue = LASColor_GetBlue(LAS_color);
- if (Centr[centr].cats->n_cats == 0) {
- nocat_area += area;
- n_nocat++;
- continue;
+ sprintf(buf, ", %d, %d, %d", red, green, blue);
+ db_append_string(&sql, buf);
+ sprintf(buf, ", \"%03d:%03d:%03d\"", red, green, blue);
+ db_append_string(&sql, buf);
}
+ db_append_string(&sql, " )");
+ G_debug(3, db_get_string(&sql));
- if (Centr[centr].cats->n_cats > 1) {
- Vect_cat_set(Centr[centr].cats, nlayers + 1,
- Centr[centr].cats->n_cats);
- overlap_area += area;
- n_overlaps++;
+ if (db_execute_immediate(driver, &sql) != DB_OK) {
+ db_close_database(driver);
+ db_shutdown_driver(driver);
+ G_fatal_error(_("Cannot insert new row: %s"),
+ db_get_string(&sql));
}
-
- Vect_reset_line(Points);
- Vect_append_point(Points, Centr[centr].x, Centr[centr].y, 0.0);
- if (type & GV_POINT)
- otype = GV_POINT;
- else
- otype = GV_CENTROID;
- Vect_write_line(&Tmp, otype, Points, Centr[centr].cats);
}
- if (Centr)
- G_free(Centr);
-
- Vect_spatial_index_destroy(&si);
- if (n_overlaps > 0) {
- G_warning(_("%d areas represent more (overlapping) features, because polygons overlap "
- "in input layer(s). Such areas are linked to more than 1 row in attribute table. "
- "The number of features for those areas is stored as category in layer %d"),
- n_overlaps, nlayers + 1);
- }
+ cat++;
+ LASPoint_Destroy(LAS_point);
+ }
+ G_percent(n_features, n_features, 1); /* finish it */
- G_message("%s", separator);
-
- Vect_hist_write(&Map, separator);
- Vect_hist_write(&Map, "\n");
- sprintf(buf, _("%d input polygons\n"), n_polygons);
- G_message(_("%d input polygons"), n_polygons);
- Vect_hist_write(&Map, buf);
-
- sprintf(buf, _("Total area: %G (%d areas)\n"), total_area, ncentr);
- G_message(_("Total area: %G (%d areas)"), total_area, ncentr);
- Vect_hist_write(&Map, buf);
-
- sprintf(buf, _("Overlapping area: %G (%d areas)\n"), overlap_area,
- n_overlaps);
- G_message(_("Overlapping area: %G (%d areas)"), overlap_area,
- n_overlaps);
- Vect_hist_write(&Map, buf);
-
- sprintf(buf, _("Area without category: %G (%d areas)\n"), nocat_area,
- n_nocat);
- G_message(_("Area without category: %G (%d areas)"), nocat_area,
- n_nocat);
- Vect_hist_write(&Map, buf);
- G_message("%s", separator);
+ if (!notab_flag->answer) {
+ db_commit_transaction(driver);
+ db_close_database_shutdown_driver(driver);
}
+
+ LASSRS_Destroy(LAS_srs);
+ LASHeader_Destroy(LAS_header);
+ LASReader_Destroy(LAS_reader);
- /* needed?
- * OGR_DS_Destroy( Ogr_ds );
- */
-
- /* Copy temporary vector to output vector */
- Vect_copy_map_lines(&Tmp, &Map);
- /* release memory occupied by topo, we may need that memory for main output */
- Vect_set_release_support(&Tmp);
- Vect_close(&Tmp);
- Vect_delete(tempvect);
-
- Vect_build(&Map);
+ /* close map */
+ if (!notopo_flag->answer)
+ Vect_build(&Map);
Vect_close(&Map);
+
+ G_message(_("%d points imported"), n_features - not_valid - n_outside);
+ if (not_valid)
+ G_message(_("%d input points were not valid"), not_valid);
+ if (n_outside)
+ G_message(_("%d input points were outside of the selected area"), n_outside);
/* -------------------------------------------------------------------- */
/* Extend current window based on dataset. */
/* -------------------------------------------------------------------- */
if (extend_flag->answer) {
- G_get_default_window(&loc_wind);
+ G_get_set_window(&loc_wind);
loc_wind.north = MAX(loc_wind.north, cellhd.north);
loc_wind.south = MIN(loc_wind.south, cellhd.south);
@@ -1200,12 +742,79 @@
/ loc_wind.ew_res);
loc_wind.east = loc_wind.west + loc_wind.cols * loc_wind.ew_res;
- G__put_window(&loc_wind, "../PERMANENT", "DEFAULT_WIND");
+ if (outloc_opt->answer)
+ G__put_window(&loc_wind, "../PERMANENT", "DEFAULT_WIND");
+ else
+ G_put_window(&loc_wind);
}
- if (with_z && !z_flag->answer)
- G_warning(_("Input data contains 3D features. Created vector is 2D only, "
- "use -z flag to import 3D vector"));
-
exit(EXIT_SUCCESS);
}
+
+void print_lasinfo(LASHeaderH LAS_header, LASSRSH LAS_srs)
+{
+ char *las_srs_proj4 = LASSRS_GetProj4(LAS_srs);
+ int las_point_format = LASHeader_GetDataFormatId(LAS_header);
+
+ fprintf(stdout, "\nUsing LAS Library Version '%s'\n\n",
+ LAS_GetFullVersion());
+ fprintf(stdout, "LAS File Version: %d.%d\n",
+ LASHeader_GetVersionMajor(LAS_header),
+ LASHeader_GetVersionMinor(LAS_header));
+ fprintf(stdout, "System ID: '%s'\n",
+ LASHeader_GetSystemId(LAS_header));
+ fprintf(stdout, "Generating Software: '%s'\n",
+ LASHeader_GetSoftwareId(LAS_header));
+ fprintf(stdout, "File Creation Day/Year: %d/%d\n",
+ LASHeader_GetCreationDOY(LAS_header),
+ LASHeader_GetCreationYear(LAS_header));
+ fprintf(stdout, "Point Data Format: %d\n",
+ las_point_format);
+ fprintf(stdout, "Number of Point Records: %d\n",
+ LASHeader_GetPointRecordsCount(LAS_header));
+ fprintf(stdout, "Number of Points by Return: %d %d %d %d %d\n",
+ LASHeader_GetPointRecordsByReturnCount(LAS_header, 0),
+ LASHeader_GetPointRecordsByReturnCount(LAS_header, 1),
+ LASHeader_GetPointRecordsByReturnCount(LAS_header, 2),
+ LASHeader_GetPointRecordsByReturnCount(LAS_header, 3),
+ LASHeader_GetPointRecordsByReturnCount(LAS_header, 4));
+ fprintf(stdout, "Scale Factor X Y Z: %g %g %g\n",
+ LASHeader_GetScaleX(LAS_header),
+ LASHeader_GetScaleY(LAS_header),
+ LASHeader_GetScaleZ(LAS_header));
+ fprintf(stdout, "Offset X Y Z: %g %g %g\n",
+ LASHeader_GetOffsetX(LAS_header),
+ LASHeader_GetOffsetY(LAS_header),
+ LASHeader_GetOffsetZ(LAS_header));
+ fprintf(stdout, "Min X Y Z: %g %g %g\n",
+ LASHeader_GetMinX(LAS_header),
+ LASHeader_GetMinY(LAS_header),
+ LASHeader_GetMinZ(LAS_header));
+ fprintf(stdout, "Max X Y Z: %g %g %g\n",
+ LASHeader_GetMaxX(LAS_header),
+ LASHeader_GetMaxY(LAS_header),
+ LASHeader_GetMaxZ(LAS_header));
+ if (las_srs_proj4 && strlen(las_srs_proj4) > 0) {
+ fprintf(stdout, "Spatial Reference:\n");
+ fprintf(stdout, "%s\n", las_srs_proj4);
+ }
+ else {
+ fprintf(stdout, "Spatial Reference: None\n");
+ }
+
+ fprintf(stdout, "\nData Fields:\n");
+ fprintf(stdout, " 'X'\n 'Y'\n 'Z'\n 'Intensity'\n 'Return Number'\n");
+ fprintf(stdout, " 'Number of Returns'\n 'Scan Direction'\n");
+ fprintf(stdout, " 'Flighline Edge'\n 'Classification'\n 'Scan Angle Rank'\n");
+ fprintf(stdout, " 'User Data'\n 'Point Source ID'\n");
+ if (las_point_format == 1 || las_point_format == 3 || las_point_format == 4 || las_point_format == 5) {
+ fprintf(stdout, " 'GPS Time'\n");
+ }
+ if (las_point_format == 2 || las_point_format == 3 || las_point_format == 5) {
+ fprintf(stdout, " 'Red'\n 'Green'\n 'Blue'\n");
+ }
+ fprintf(stdout, "\n");
+ fflush(stdout);
+
+ return;
+}
More information about the grass-commit
mailing list