main.c 4.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187
  1. /****************************************************************************
  2. *
  3. * MODULE: r.circle
  4. *
  5. * AUTHOR(S): Bill Brown - CERL (Jan, 1993)
  6. * Markus Neteler
  7. *
  8. * PURPOSE: Creates a raster map containing concentric rings
  9. * around a given point.
  10. *
  11. * COPYRIGHT: (C) 2006-2008 by the GRASS Development Team
  12. *
  13. * This program is free software under the GNU General Public
  14. * License (>=v2). Read the file COPYING that comes with GRASS
  15. * for details.
  16. *
  17. ***************************************************************************/
  18. #include <stdlib.h>
  19. #include <strings.h>
  20. #include <math.h>
  21. #include <grass/gis.h>
  22. #include <grass/raster.h>
  23. #include <grass/glocale.h>
  24. static double distance(double *, double *, double, double, int);
  25. #ifndef HUGE_VAL
  26. #define HUGE_VAL 1.7976931348623157e+308
  27. #endif
  28. int main(int argc, char *argv[])
  29. {
  30. struct GModule *module;
  31. struct Option *coord, *out_file, *min, *max, *mult;
  32. struct Flag *flag;
  33. int *int_buf;
  34. struct Cell_head w;
  35. struct History history;
  36. int cellfile;
  37. double east, north, pt[2], cur[2], row, col, fmult;
  38. double fmin, fmax;
  39. int binary;
  40. G_gisinit(argv[0]);
  41. module = G_define_module();
  42. G_add_keyword(_("raster"));
  43. G_add_keyword(_("geometry"));
  44. module->description =
  45. _("Creates a raster map containing concentric "
  46. "rings around a given point.");
  47. out_file = G_define_standard_option(G_OPT_R_OUTPUT);
  48. coord = G_define_option();
  49. coord->key = "coordinate";
  50. coord->type = TYPE_STRING;
  51. coord->required = YES;
  52. coord->key_desc = "x,y";
  53. coord->description = _("The coordinate of the center (east,north)");
  54. min = G_define_option();
  55. min->key = "min";
  56. min->type = TYPE_DOUBLE;
  57. min->required = NO;
  58. min->description = _("Minimum radius for ring/circle map (in meters)");
  59. max = G_define_option();
  60. max->key = "max";
  61. max->type = TYPE_DOUBLE;
  62. max->required = NO;
  63. max->description = _("Maximum radius for ring/circle map (in meters)");
  64. mult = G_define_option();
  65. mult->key = "mult";
  66. mult->type = TYPE_DOUBLE;
  67. mult->required = NO;
  68. mult->description = _("Data value multiplier");
  69. flag = G_define_flag();
  70. flag->key = 'b';
  71. flag->description = _("Generate binary raster map");
  72. if (G_parser(argc, argv))
  73. exit(EXIT_FAILURE);
  74. G_scan_easting(coord->answers[0], &east, G_projection());
  75. G_scan_northing(coord->answers[1], &north, G_projection());
  76. pt[0] = east;
  77. pt[1] = north;
  78. fmult = 1.0;
  79. if (min->answer)
  80. sscanf(min->answer, "%lf", &fmin);
  81. else
  82. fmin = 0;
  83. if (max->answer)
  84. sscanf(max->answer, "%lf", &fmax);
  85. else
  86. fmax = HUGE_VAL;
  87. if (fmin > fmax)
  88. G_fatal_error(_("Please specify a radius in which min < max"));
  89. if (mult->answer)
  90. if (1 != sscanf(mult->answer, "%lf", &fmult))
  91. fmult = 1.0;
  92. /* nonsense test */
  93. if (flag->answer && (!min->answer && !max->answer))
  94. G_fatal_error(_("Please specify min and/or max radius when "
  95. "using the binary flag"));
  96. if (flag->answer)
  97. binary = 1; /* generate binary pattern only, useful for MASK */
  98. else
  99. binary = 0;
  100. G_get_set_window(&w);
  101. cellfile = Rast_open_c_new(out_file->answer);
  102. int_buf = (int *)G_malloc(w.cols * sizeof(int));
  103. {
  104. int c;
  105. for (row = 0; row < w.rows; row++) {
  106. G_percent(row, w.rows, 2);
  107. cur[1] = Rast_row_to_northing(row + 0.5, &w);
  108. for (col = 0; col < w.cols; col++) {
  109. c = col;
  110. cur[0] = Rast_col_to_easting(col + 0.5, &w);
  111. int_buf[c] =
  112. (int)(distance(pt, cur, fmin, fmax, binary) * fmult);
  113. if (int_buf[c] == 0)
  114. Rast_set_null_value(&int_buf[c], 1, CELL_TYPE);
  115. }
  116. Rast_put_row(cellfile, int_buf, CELL_TYPE);
  117. }
  118. }
  119. G_free(int_buf);
  120. Rast_close(cellfile);
  121. Rast_short_history(out_file->answer, "raster", &history);
  122. Rast_command_history(&history);
  123. Rast_write_history(out_file->answer, &history);
  124. G_done_msg(_("Raster map <%s> created."),
  125. out_file->answer);
  126. return (EXIT_SUCCESS);
  127. }
  128. /*******************************************************************/
  129. static double distance(double from[2], double to[2], double min, double max,
  130. int binary)
  131. {
  132. static int first = 1;
  133. double dist;
  134. if (first) {
  135. first = 0;
  136. G_begin_distance_calculations();
  137. }
  138. dist = G_distance(from[0], from[1], to[0], to[1]);
  139. if ((dist >= min) && (dist <= max))
  140. if (!binary)
  141. return dist;
  142. else
  143. return 1;
  144. else
  145. return 0; /* should be NULL ? */
  146. }
  147. /**********************************************************************/