main.c 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348
  1. /****************************************************************************
  2. *
  3. * MODULE: r3.cross.rast
  4. *
  5. * AUTHOR(S): Original author
  6. * Soeren Gebbert soerengebbert at gmx de
  7. * 23 Feb 2006 Berlin
  8. * PURPOSE: Creates a cross section 2D map from one RASTER3D raster map based on a 2D elevation map
  9. *
  10. * COPYRIGHT: (C) 2005 by the GRASS Development Team
  11. *
  12. * This program is free software under the GNU General Public
  13. * License (>=v2). Read the file COPYING that comes with GRASS
  14. * for details.
  15. *
  16. *****************************************************************************/
  17. #include <stdio.h>
  18. #include <stdlib.h>
  19. #include <string.h>
  20. #include <grass/gis.h>
  21. #include <grass/raster.h>
  22. #include <grass/raster3d.h>
  23. #include <grass/glocale.h>
  24. /*- param and global variables -----------------------------------------*/
  25. typedef struct {
  26. struct Option *input, *output, *elevation;
  27. struct Flag *mask;
  28. } paramType;
  29. paramType param; /*param */
  30. int globalElevMapType;
  31. /*- prototypes --------------------------------------------------------------*/
  32. void fatal_error(void *map, int elevfd, int outfd, char *errorMsg); /*Simple Error message */
  33. void set_params(); /*Fill the paramType structure */
  34. void rast3d_cross_section(void *map, RASTER3D_Region region, int elevfd, int outfd); /*Write the raster */
  35. void close_output_map(int fd); /*close the map */
  36. /* ************************************************************************* */
  37. /* Error handling ********************************************************** */
  38. /* ************************************************************************* */
  39. void fatal_error(void *map, int elevfd, int outfd, char *errorMsg)
  40. {
  41. /* Close files and exit */
  42. if (map != NULL) {
  43. if (!Rast3d_close(map))
  44. Rast3d_fatal_error(_("Unable to close 3D raster map"));
  45. }
  46. /*unopen the output map */
  47. if (outfd != -1)
  48. Rast_unopen(outfd);
  49. if (elevfd != -1)
  50. close_output_map(elevfd);
  51. Rast3d_fatal_error("%s", errorMsg);
  52. exit(EXIT_FAILURE);
  53. }
  54. /* ************************************************************************* */
  55. /* Close the raster map ********************************************* */
  56. /* ************************************************************************* */
  57. void close_output_map(int fd)
  58. {
  59. Rast_close(fd);
  60. }
  61. /* ************************************************************************* */
  62. /* Set up the arguments we are expecting *********************************** */
  63. /* ************************************************************************* */
  64. void set_params()
  65. {
  66. param.input = G_define_option();
  67. param.input->key = "input";
  68. param.input->type = TYPE_STRING;
  69. param.input->required = YES;
  70. param.input->gisprompt = "old,grid3,3d-raster";
  71. param.input->description = _("Input 3D raster map for cross section");
  72. param.elevation = G_define_option();
  73. param.elevation->key = "elevation";
  74. param.elevation->type = TYPE_STRING;
  75. param.elevation->required = YES;
  76. param.elevation->description =
  77. _("2D elevation map used to create the cross section map");
  78. param.elevation->gisprompt = "old,cell,raster";
  79. param.output = G_define_option();
  80. param.output->key = "output";
  81. param.output->type = TYPE_STRING;
  82. param.output->required = YES;
  83. param.output->description = _("Resulting cross section 2D raster map");
  84. param.output->gisprompt = "new,cell,raster";
  85. param.mask = G_define_flag();
  86. param.mask->key = 'm';
  87. param.mask->description = _("Use 3D raster mask (if exists) with input map");
  88. }
  89. /* ************************************************************************* */
  90. /* Compute the cross section raster map ************************************ */
  91. /* ************************************************************************* */
  92. void rast3d_cross_section(void *map,RASTER3D_Region region, int elevfd, int outfd)
  93. {
  94. int col, row;
  95. int rows, cols, depths, typeIntern;
  96. FCELL *fcell = NULL;
  97. DCELL *dcell = NULL;
  98. void *elevrast;
  99. void *ptr;
  100. int intvalue;
  101. float fvalue;
  102. double dvalue;
  103. double elevation = 0;
  104. double north, east;
  105. struct Cell_head window;
  106. Rast_get_window(&window);
  107. rows = region.rows;
  108. cols = region.cols;
  109. depths = region.depths;
  110. /*Typ of the RASTER3D Tile */
  111. typeIntern = Rast3d_tile_type_map(map);
  112. /*Allocate mem for the output maps row */
  113. if (typeIntern == FCELL_TYPE)
  114. fcell = Rast_allocate_f_buf();
  115. else if (typeIntern == DCELL_TYPE)
  116. dcell = Rast_allocate_d_buf();
  117. /*Mem for the input map row */
  118. elevrast = Rast_allocate_buf(globalElevMapType);
  119. for (row = 0; row < rows; row++) {
  120. G_percent(row, rows - 1, 10);
  121. /*Read the input map row */
  122. Rast_get_row(elevfd, elevrast, row, globalElevMapType);
  123. for (col = 0, ptr = elevrast; col < cols; col++, ptr =
  124. G_incr_void_ptr(ptr, Rast_cell_size(globalElevMapType))) {
  125. if (Rast_is_null_value(ptr, globalElevMapType)) {
  126. if (typeIntern == FCELL_TYPE)
  127. Rast_set_null_value(&fcell[col], 1, FCELL_TYPE);
  128. else if (typeIntern == DCELL_TYPE)
  129. Rast_set_null_value(&dcell[col], 1, DCELL_TYPE);
  130. continue;
  131. }
  132. /*Read the elevation value */
  133. if (globalElevMapType == CELL_TYPE) {
  134. intvalue = *(CELL *) ptr;
  135. elevation = intvalue;
  136. } else if (globalElevMapType == FCELL_TYPE) {
  137. fvalue = *(FCELL *) ptr;
  138. elevation = fvalue;
  139. } else if (globalElevMapType == DCELL_TYPE) {
  140. dvalue = *(DCELL *) ptr;
  141. elevation = dvalue;
  142. }
  143. /* Compute the coordinates */
  144. north = Rast_row_to_northing((double)row + 0.5, &window);
  145. east = Rast_col_to_easting((double)col + 0.5, &window);
  146. /* Get the voxel value */
  147. if (typeIntern == FCELL_TYPE)
  148. Rast3d_get_region_value(map, north, east, elevation, &fcell[col], FCELL_TYPE);
  149. if (typeIntern == DCELL_TYPE)
  150. Rast3d_get_region_value(map, north, east, elevation, &dcell[col], DCELL_TYPE);
  151. }
  152. /*Write the data to the output map */
  153. if (typeIntern == FCELL_TYPE)
  154. Rast_put_f_row(outfd, fcell);
  155. if (typeIntern == DCELL_TYPE)
  156. Rast_put_d_row(outfd, dcell);
  157. }
  158. G_debug(3, "\nDone\n");
  159. /*Free the mem */
  160. if (elevrast)
  161. G_free(elevrast);
  162. if (dcell)
  163. G_free(dcell);
  164. if (fcell)
  165. G_free(fcell);
  166. }
  167. /* ************************************************************************* */
  168. /* Main function, open the RASTER3D map and create the cross section map ******** */
  169. /* ************************************************************************* */
  170. int main(int argc, char *argv[])
  171. {
  172. RASTER3D_Region region;
  173. struct Cell_head window2d;
  174. struct GModule *module;
  175. void *map = NULL; /*The 3D Rastermap */
  176. int changemask = 0;
  177. int elevfd = -1, outfd = -1; /*file descriptors */
  178. int output_type, cols, rows;
  179. /* Initialize GRASS */
  180. G_gisinit(argv[0]);
  181. module = G_define_module();
  182. G_add_keyword(_("raster3d"));
  183. G_add_keyword(_("profile"));
  184. G_add_keyword(_("raster"));
  185. G_add_keyword(_("voxel"));
  186. module->description =
  187. _("Creates cross section 2D raster map from 3D raster map based on 2D elevation map");
  188. /* Get parameters from user */
  189. set_params();
  190. /* Have GRASS get inputs */
  191. if (G_parser(argc, argv))
  192. exit(EXIT_FAILURE);
  193. G_debug(3, "Open 3D raster map %s", param.input->answer);
  194. if (NULL == G_find_raster3d(param.input->answer, ""))
  195. Rast3d_fatal_error(_("3D raster map <%s> not found"),
  196. param.input->answer);
  197. /* Figure out the region from the map */
  198. Rast3d_init_defaults();
  199. Rast3d_get_window(&region);
  200. /*Check if the g3d-region is equal to the 2d rows and cols */
  201. rows = Rast_window_rows();
  202. cols = Rast_window_cols();
  203. /*If not equal, set the 2D windows correct */
  204. if (rows != region.rows || cols != region.cols) {
  205. G_message
  206. (_("The 2D and 3D region settings are different. Using the 3D raster map settings to adjust the 2D region."));
  207. G_get_set_window(&window2d);
  208. window2d.ns_res = region.ns_res;
  209. window2d.ew_res = region.ew_res;
  210. window2d.rows = region.rows;
  211. window2d.cols = region.cols;
  212. Rast_set_window(&window2d);
  213. }
  214. /*******************/
  215. /*Open the 3d raster map */
  216. /*******************/
  217. map = Rast3d_open_cell_old(param.input->answer,
  218. G_find_raster3d(param.input->answer, ""),
  219. &region, RASTER3D_TILE_SAME_AS_FILE,
  220. RASTER3D_USE_CACHE_DEFAULT);
  221. if (map == NULL)
  222. Rast3d_fatal_error(_("Unable to open 3D raster map <%s>"),
  223. param.input->answer);
  224. /*Get the output type */
  225. output_type = Rast3d_file_type_map(map);
  226. if (output_type == FCELL_TYPE || output_type == DCELL_TYPE) {
  227. /********************************/
  228. /*Open the elevation raster map */
  229. /********************************/
  230. elevfd = Rast_open_old(param.elevation->answer, "");
  231. globalElevMapType = Rast_get_map_type(elevfd);
  232. /**********************/
  233. /*Open the Outputmap */
  234. /**********************/
  235. if (G_find_raster2(param.output->answer, ""))
  236. G_message(_("Output map already exists. Will be overwritten!"));
  237. if (output_type == FCELL_TYPE)
  238. outfd = Rast_open_new(param.output->answer, FCELL_TYPE);
  239. else if (output_type == DCELL_TYPE)
  240. outfd = Rast_open_new(param.output->answer, DCELL_TYPE);
  241. /*if requested set the Mask on */
  242. if (param.mask->answer) {
  243. if (Rast3d_mask_file_exists()) {
  244. changemask = 0;
  245. if (Rast3d_mask_is_off(map)) {
  246. Rast3d_mask_on(map);
  247. changemask = 1;
  248. }
  249. }
  250. }
  251. /************************/
  252. /*Create the Rastermaps */
  253. /************************/
  254. rast3d_cross_section(map, region, elevfd, outfd);
  255. /*We set the Mask off, if it was off before */
  256. if (param.mask->answer) {
  257. if (Rast3d_mask_file_exists())
  258. if (Rast3d_mask_is_on(map) && changemask)
  259. Rast3d_mask_off(map);
  260. }
  261. Rast_close(outfd);
  262. Rast_close(elevfd);
  263. } else {
  264. fatal_error(map, -1, -1,
  265. _("Wrong 3D raster datatype! Cannot create raster map"));
  266. }
  267. /* Close files and exit */
  268. if (!Rast3d_close(map))
  269. Rast3d_fatal_error(_("Unable to close 3D raster map <%s>"),
  270. param.input->answer);
  271. return (EXIT_SUCCESS);
  272. }