1 #include <stdio.h>
2 #include <math.h>
3
4 #include <grass/raster.h>
5 #include "raster3d_intern.h"
6
7
8 /*---------------------------------------------------------------------------*/
9
10
11 /*!
12 * \brief
13 *
14 * Returns in <em>region2d</em> the <em>2d</em> portion of <em>region3d</em>.
15 *
16 * \param region3d
17 * \param region2d
18 * \return void
19 */
20
Rast3d_extract2d_region(RASTER3D_Region * region3d,struct Cell_head * region2d)21 void Rast3d_extract2d_region(RASTER3D_Region * region3d, struct Cell_head *region2d)
22 {
23 region2d->proj = region3d->proj;
24 region2d->zone = region3d->zone;
25
26 region2d->north = region3d->north;
27 region2d->south = region3d->south;
28 region2d->east = region3d->east;
29 region2d->west = region3d->west;
30
31 region2d->rows = region3d->rows;
32 region2d->cols = region3d->cols;
33
34 region2d->ns_res = region3d->ns_res;
35 region2d->ew_res = region3d->ew_res;
36 }
37
38 /*!
39 * \brief
40 *
41 * Returns in <em>region2d</em> the <em>2d</em> portion of <em>region3d</em>.
42 *
43 * \param region3d
44 * \param region2d
45 * \return void
46 */
47
Rast3d_region_to_cell_head(RASTER3D_Region * region3d,struct Cell_head * region2d)48 void Rast3d_region_to_cell_head(RASTER3D_Region * region3d, struct Cell_head *region2d)
49 {
50 region2d->proj = region3d->proj;
51 region2d->zone = region3d->zone;
52
53 region2d->north = region3d->north;
54 region2d->south = region3d->south;
55 region2d->east = region3d->east;
56 region2d->west = region3d->west;
57 region2d->top = region3d->top;
58 region2d->bottom = region3d->bottom;
59
60 region2d->rows = region3d->rows;
61 region2d->rows3 = region3d->rows;
62 region2d->cols = region3d->cols;
63 region2d->cols3 = region3d->cols;
64 region2d->depths = region3d->depths;
65
66 region2d->ns_res = region3d->ns_res;
67 region2d->ns_res3 = region3d->ns_res;
68 region2d->ew_res = region3d->ew_res;
69 region2d->ew_res3 = region3d->ew_res;
70 region2d->tb_res = region3d->tb_res;
71 }
72
73 /*---------------------------------------------------------------------------*/
74
75
76 /*!
77 * \brief
78 *
79 * Replaces the <em>2d</em> portion of <em>region3d</em> with the
80 * values stored in <em>region2d</em>.
81 *
82 * \param region2d
83 * \param region3d
84 * \return void
85 */
86
87 void
Rast3d_incorporate2d_region(struct Cell_head * region2d,RASTER3D_Region * region3d)88 Rast3d_incorporate2d_region(struct Cell_head *region2d, RASTER3D_Region * region3d)
89 {
90 region3d->proj = region2d->proj;
91 region3d->zone = region2d->zone;
92
93 region3d->north = region2d->north;
94 region3d->south = region2d->south;
95 region3d->east = region2d->east;
96 region3d->west = region2d->west;
97
98 region3d->rows = region2d->rows;
99 region3d->cols = region2d->cols;
100
101 region3d->ns_res = region2d->ns_res;
102 region3d->ew_res = region2d->ew_res;
103 }
104
105 /*!
106 * \brief
107 *
108 * Replaces the <em>2d</em> portion of <em>region3d</em> with the
109 * values stored in <em>region2d</em>.
110 *
111 * \param region2d
112 * \param region3d
113 * \return void
114 */
115
116 void
Rast3d_region_from_to_cell_head(struct Cell_head * region2d,RASTER3D_Region * region3d)117 Rast3d_region_from_to_cell_head(struct Cell_head *region2d, RASTER3D_Region * region3d)
118 {
119 region3d->proj = region2d->proj;
120 region3d->zone = region2d->zone;
121
122 region3d->north = region2d->north;
123 region3d->south = region2d->south;
124 region3d->east = region2d->east;
125 region3d->west = region2d->west;
126 region3d->top = region2d->top;
127 region3d->bottom = region2d->bottom;
128
129 region3d->rows = region2d->rows3;
130 region3d->cols = region2d->cols3;
131 region3d->depths = region2d->depths;
132
133 region3d->ns_res = region2d->ns_res3;
134 region3d->ew_res = region2d->ew_res3;
135 region3d->tb_res = region2d->tb_res;
136 }
137
138 /*---------------------------------------------------------------------------*/
139
140
141 /*!
142 * \brief
143 *
144 * Computes an adjusts the resolutions in the region structure from the region
145 * boundaries and number of cells per dimension.
146 *
147 * \param region
148 * \return void
149 */
150
Rast3d_adjust_region(RASTER3D_Region * region)151 void Rast3d_adjust_region(RASTER3D_Region * region)
152 {
153 struct Cell_head region2d;
154
155 Rast3d_region_to_cell_head(region, ®ion2d);
156 G_adjust_Cell_head3(®ion2d, 1, 1, 1);
157 Rast3d_region_from_to_cell_head(®ion2d, region);
158
159 if (region->depths <= 0)
160 Rast3d_fatal_error("Rast3d_adjust_region: depths <= 0");
161 region->tb_res = (region->top - region->bottom) / region->depths;
162 }
163
164 /*---------------------------------------------------------------------------*/
165
166
167 /*!
168 * \brief
169 *
170 * Computes an adjusts the number of cells per dimension in the region
171 * structure from the region boundaries and resolutions.
172 *
173 * \param region
174 * \return void
175 */
176
Rast3d_adjust_region_res(RASTER3D_Region * region)177 void Rast3d_adjust_region_res(RASTER3D_Region * region)
178 {
179 struct Cell_head region2d;
180
181 Rast3d_region_to_cell_head(region, ®ion2d);
182 G_adjust_Cell_head3(®ion2d, 1, 1, 1);
183 Rast3d_region_from_to_cell_head(®ion2d, region);
184
185 if (region->tb_res <= 0)
186 Rast3d_fatal_error("Rast3d_adjust_region_res: tb_res <= 0");
187
188 region->depths = (region->top - region->bottom + region->tb_res / 2.0) /
189 region->tb_res;
190 if (region->depths == 0)
191 region->depths = 1;
192 }
193
194 /*---------------------------------------------------------------------------*/
195
196
197 /*!
198 * \brief
199 *
200 * Copies the values of <em>regionSrc</em> into <em>regionDst</em>.
201 *
202 * \param regionDest
203 * \param regionSrc
204 * \return void
205 */
206
Rast3d_region_copy(RASTER3D_Region * regionDest,RASTER3D_Region * regionSrc)207 void Rast3d_region_copy(RASTER3D_Region * regionDest, RASTER3D_Region * regionSrc)
208 {
209 regionDest->proj = regionSrc->proj;
210 regionDest->zone = regionSrc->zone;
211
212 regionDest->north = regionSrc->north;
213 regionDest->south = regionSrc->south;
214 regionDest->east = regionSrc->east;
215 regionDest->west = regionSrc->west;
216 regionDest->top = regionSrc->top;
217 regionDest->bottom = regionSrc->bottom;
218
219 regionDest->rows = regionSrc->rows;
220 regionDest->cols = regionSrc->cols;
221 regionDest->depths = regionSrc->depths;
222
223 regionDest->ns_res = regionSrc->ns_res;
224 regionDest->ew_res = regionSrc->ew_res;
225 regionDest->tb_res = regionSrc->tb_res;
226 }
227
228
229 /*---------------------------------------------------------------------------*/
230
231 int
Rast3d_read_region_map(const char * name,const char * mapset,RASTER3D_Region * region)232 Rast3d_read_region_map(const char *name, const char *mapset, RASTER3D_Region * region)
233 {
234 char fullName[GPATH_MAX];
235 char xname[GNAME_MAX], xmapset[GMAPSET_MAX];
236
237 if (G_name_is_fully_qualified(name, xname, xmapset))
238 Rast3d_filename(fullName, RASTER3D_HEADER_ELEMENT, xname, xmapset);
239 else {
240 if (!mapset || !*mapset)
241 mapset = G_find_raster3d(name, "");
242 Rast3d_filename(fullName, RASTER3D_HEADER_ELEMENT, name, mapset);
243 }
244 return Rast3d_read_window(region, fullName);
245 }
246
247 /*---------------------------------------------------------------------------*/
248
249
250 /*!
251 * \brief
252 *
253 * Returns 1 if region-coordinates <em>(north, east, top)</em> are
254 * inside the region of <em>map</em>. Returns 0 otherwise.
255 *
256 * \param region
257 * \param north
258 * \param east
259 * \param top
260 * \return int
261 */
262
Rast3d_is_valid_location(RASTER3D_Region * region,double north,double east,double top)263 int Rast3d_is_valid_location(RASTER3D_Region *region, double north, double east, double top)
264 {
265 return ((north >= region->south) && (north <= region->north) &&
266 (east >= region->west) && (east <= region->east) &&
267 (((top >= region->bottom) && (top <= region->top)) ||
268 ((top <= region->bottom) && (top >= region->top))));
269 }
270
271 /*---------------------------------------------------------------------------*/
272
273 /*!
274 * \brief
275 *
276 * Converts region-coordinates <em>(north, east,
277 * top)</em> into cell-coordinates <em>(x, y, z)</em>.
278 *
279 * \param region
280 * \param north
281 * \param east
282 * \param top
283 * \param x
284 * \param y
285 * \param z
286 * \return void
287 */
288
289 void
Rast3d_location2coord(RASTER3D_Region * region,double north,double east,double top,int * x,int * y,int * z)290 Rast3d_location2coord(RASTER3D_Region *region, double north, double east, double top,
291 int *x, int *y, int *z)
292 {
293 double col, row, depth;
294 LOCATION_TO_COORD(region, north, east, top, &col, &row, &depth);
295
296 *x = (int)floor(col);
297 *y = (int)floor(row);
298 *z = (int)floor(depth);
299 }
300
301 /*!
302 * \brief
303 *
304 * Converts region-coordinates <em>(north, east,
305 * top)</em> into cell-coordinates <em>(x, y, z)</em>.
306 *
307 * <b>Note:</b> The results are <i>double</i> numbers. Casting them to
308 * <i>int</i> will give the column, row and depth number.
309 *
310 * \param region
311 * \param north
312 * \param east
313 * \param top
314 * \param x
315 * \param y
316 * \param z
317 * \return void
318 */
319
320 void
Rast3d_location2coord_double(RASTER3D_Region * region,double north,double east,double top,double * x,double * y,double * z)321 Rast3d_location2coord_double(RASTER3D_Region *region, double north, double east, double top,
322 double *x, double *y, double *z)
323 {
324 LOCATION_TO_COORD(region, north, east, top, x, y, z);
325
326 G_debug(4, "Rast3d_location2coord_double x %f y %f z %f\n", *x, *y, *z);
327 }
328
329 /*!
330 * \brief
331 *
332 * Converts region-coordinates <em>(north, east,
333 * top)</em> into cell-coordinates <em>(x, y, z)</em>.
334 * This function calls Rast3d_fatal_error in case location is not in window.
335 *
336 * \param region
337 * \param north
338 * \param east
339 * \param top
340 * \param x
341 * \param y
342 * \param z
343 * \return void
344 */
345
346 void
Rast3d_location2coord2(RASTER3D_Region * region,double north,double east,double top,int * x,int * y,int * z)347 Rast3d_location2coord2(RASTER3D_Region *region, double north, double east, double top,
348 int *x, int *y, int *z)
349 {
350 if (!Rast3d_is_valid_location(region, north, east, top))
351 Rast3d_fatal_error("Rast3d_location2coord2: location not in region");
352
353 double col, row, depth;
354 LOCATION_TO_COORD(region, north, east, top, &col, &row, &depth);
355
356 *x = (int)floor(col);
357 *y = (int)floor(row);
358 *z = (int)floor(depth);
359 }
360
361 /*!
362 * \brief
363 *
364 * Converts cell-coordinates <em>(x, y, z)</em> into region-coordinates
365 * <em>(north, east, top)</em>.
366 *
367 * * <b>Note:</b> x, y and z is a double:
368 * - x+0.0 will return the easting for the western edge of the column.
369 * - x+0.5 will return the easting for the center of the column.
370 * - x+1.0 will return the easting for the eastern edge of the column.
371 *
372 * - y+0.0 will return the northing for the northern edge of the row.
373 * - y+0.5 will return the northing for the center of the row.
374 * - y+1.0 will return the northing for the southern edge of the row.
375 *
376 * - z+0.0 will return the top for the lower edge of the depth.
377 * - z+0.5 will return the top for the center of the depth.
378 * - z+1.0 will return the top for the upper edge of the column.
379 *
380 *
381 * \param region
382 * \param x
383 * \param y
384 * \param z
385 * \param north
386 * \param east
387 * \param top
388 * \return void
389 */
390
391 void
Rast3d_coord2location(RASTER3D_Region * region,double x,double y,double z,double * north,double * east,double * top)392 Rast3d_coord2location(RASTER3D_Region * region, double x, double y, double z, double *north, double *east, double *top)
393 {
394 COORD_TO_LOCATION(region, x, y, z, north, east, top);
395
396 G_debug(4, "Rast3d_coord2location north %g east %g top %g\n", *north, *east, *top);
397 }
398