1 /*
2 * r.in.lidar projection-related functions
3 *
4 * Copyright 2011-2015 by Markus Metz, and The GRASS Development Team
5 * Authors:
6 * Markus Metz (r.in.lidar)
7 * Vaclav Petras (move functions to a file)
8 *
9 * This program is free software licensed under the GPL (>=v2).
10 * Read the COPYING file that comes with GRASS for details.
11 *
12 */
13
14 #include <string.h>
15
16 #include <grass/glocale.h>
17
18 #include <liblas/capi/liblas.h>
19
20 #include "local_proto.h"
21
22
print_lasinfo(LASHeaderH LAS_header,LASSRSH LAS_srs)23 void print_lasinfo(LASHeaderH LAS_header, LASSRSH LAS_srs)
24 {
25 char *las_srs_proj4 = LASSRS_GetProj4(LAS_srs);
26 int las_point_format = LASHeader_GetDataFormatId(LAS_header);
27
28 fprintf(stdout, "\nUsing LAS Library Version '%s'\n\n",
29 LAS_GetFullVersion());
30 fprintf(stdout, "LAS File Version: %d.%d\n",
31 LASHeader_GetVersionMajor(LAS_header),
32 LASHeader_GetVersionMinor(LAS_header));
33 fprintf(stdout, "System ID: '%s'\n",
34 LASHeader_GetSystemId(LAS_header));
35 fprintf(stdout, "Generating Software: '%s'\n",
36 LASHeader_GetSoftwareId(LAS_header));
37 fprintf(stdout, "File Creation Day/Year: %d/%d\n",
38 LASHeader_GetCreationDOY(LAS_header),
39 LASHeader_GetCreationYear(LAS_header));
40 fprintf(stdout, "Point Data Format: %d\n",
41 las_point_format);
42 fprintf(stdout, "Number of Point Records: %d\n",
43 LASHeader_GetPointRecordsCount(LAS_header));
44 fprintf(stdout, "Number of Points by Return: %d %d %d %d %d\n",
45 LASHeader_GetPointRecordsByReturnCount(LAS_header, 0),
46 LASHeader_GetPointRecordsByReturnCount(LAS_header, 1),
47 LASHeader_GetPointRecordsByReturnCount(LAS_header, 2),
48 LASHeader_GetPointRecordsByReturnCount(LAS_header, 3),
49 LASHeader_GetPointRecordsByReturnCount(LAS_header, 4));
50 fprintf(stdout, "Scale Factor X Y Z: %g %g %g\n",
51 LASHeader_GetScaleX(LAS_header),
52 LASHeader_GetScaleY(LAS_header), LASHeader_GetScaleZ(LAS_header));
53 fprintf(stdout, "Offset X Y Z: %g %g %g\n",
54 LASHeader_GetOffsetX(LAS_header),
55 LASHeader_GetOffsetY(LAS_header),
56 LASHeader_GetOffsetZ(LAS_header));
57 fprintf(stdout, "Min X Y Z: %g %g %g\n",
58 LASHeader_GetMinX(LAS_header),
59 LASHeader_GetMinY(LAS_header), LASHeader_GetMinZ(LAS_header));
60 fprintf(stdout, "Max X Y Z: %g %g %g\n",
61 LASHeader_GetMaxX(LAS_header),
62 LASHeader_GetMaxY(LAS_header), LASHeader_GetMaxZ(LAS_header));
63 if (las_srs_proj4 && strlen(las_srs_proj4) > 0) {
64 fprintf(stdout, "Spatial Reference:\n");
65 fprintf(stdout, "%s\n", las_srs_proj4);
66 }
67 else {
68 fprintf(stdout, "Spatial Reference: None\n");
69 }
70
71 fprintf(stdout, "\nData Fields:\n");
72 fprintf(stdout,
73 " 'X'\n 'Y'\n 'Z'\n 'Intensity'\n 'Return Number'\n");
74 fprintf(stdout, " 'Number of Returns'\n 'Scan Direction'\n");
75 fprintf(stdout,
76 " 'Flighline Edge'\n 'Classification'\n 'Scan Angle Rank'\n");
77 fprintf(stdout, " 'User Data'\n 'Point Source ID'\n");
78 if (las_point_format == 1 || las_point_format == 3 ||
79 las_point_format == 4 || las_point_format == 5) {
80 fprintf(stdout, " 'GPS Time'\n");
81 }
82 if (las_point_format == 2 || las_point_format == 3 ||
83 las_point_format == 5) {
84 fprintf(stdout, " 'Red'\n 'Green'\n 'Blue'\n");
85 }
86 fprintf(stdout, "\n");
87 fflush(stdout);
88
89 return;
90 }
91
92
scan_bounds(LASReaderH LAS_reader,int shell_style,int extents,int update,double zscale,struct Cell_head * region)93 int scan_bounds(LASReaderH LAS_reader, int shell_style, int extents, int update,
94 double zscale, struct Cell_head *region)
95 {
96 unsigned long line;
97 int first;
98 double min_x, max_x, min_y, max_y, min_z, max_z;
99 double x, y, z;
100 LASPointH LAS_point;
101
102 line = 0;
103 first = TRUE;
104
105 /* init to nan in case no points are found */
106 min_x = max_x = min_y = max_y = min_z = max_z = 0.0 / 0.0;
107
108 G_verbose_message(_("Scanning data ..."));
109
110 LASReader_Seek(LAS_reader, 0);
111
112 while ((LAS_point = LASReader_GetNextPoint(LAS_reader)) != NULL) {
113 line++;
114
115 /* we don't do any filtering here */
116
117 x = LASPoint_GetX(LAS_point);
118 y = LASPoint_GetY(LAS_point);
119 z = LASPoint_GetZ(LAS_point);
120
121 if (first) {
122 min_x = x;
123 max_x = x;
124 min_y = y;
125 max_y = y;
126 min_z = z;
127 max_z = z;
128 first = FALSE;
129 }
130 else {
131 if (x < min_x)
132 min_x = x;
133 if (x > max_x)
134 max_x = x;
135 if (y < min_y)
136 min_y = y;
137 if (y > max_y)
138 max_y = y;
139 if (z < min_z)
140 min_z = z;
141 if (z > max_z)
142 max_z = z;
143 }
144 }
145
146 if (!extents) {
147 if (!shell_style) {
148 fprintf(stderr, _("Range: min max\n"));
149 fprintf(stdout, "x: %11f %11f\n", min_x, max_x);
150 fprintf(stdout, "y: %11f %11f\n", min_y, max_y);
151 fprintf(stdout, "z: %11f %11f\n", min_z * zscale, max_z * zscale);
152 }
153 else
154 fprintf(stdout, "n=%f s=%f e=%f w=%f b=%f t=%f\n",
155 max_y, min_y, max_x, min_x, min_z * zscale,
156 max_z * zscale);
157
158 G_debug(1, "Processed %lu points.", line);
159 G_debug(1, "region template: g.region n=%f s=%f e=%f w=%f",
160 max_y, min_y, max_x, min_x);
161 }
162 else if (update) {
163 if (min_x < region->west)
164 region->west = min_x;
165 if (max_x > region->east)
166 region->east = max_x;
167 if (min_y < region->south)
168 region->south = min_y;
169 if (max_y > region->north)
170 region->north = max_y;
171 }
172 else {
173 region->east = max_x;
174 region->west = min_x;
175 region->north = max_y;
176 region->south = min_y;
177 }
178
179 return 0;
180 }
181