main.c 8.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292
  1. /*-
  2. * s.sample - GRASS program to sample a raster map at site locations.
  3. * Copyright (C) 1994. James Darrell McCauley.
  4. *
  5. * Author: James Darrell McCauley darrell@mccauley-usa.com
  6. * http://mccauley-usa.com/
  7. *
  8. * This program is free software; you can redistribute it and/or
  9. * modify it under the terms of the GNU General Public License
  10. * as published by the Free Software Foundation; either version 2
  11. * of the License, or (at your option) any later version.
  12. *
  13. * This program is distributed in the hope that it will be useful,
  14. * but WITHOUT ANY WARRANTY; without even the implied warranty of
  15. * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  16. * GNU General Public License for more details.
  17. *
  18. * You should have received a copy of the GNU General Public License
  19. * along with this program; if not, write to the Free Software
  20. * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
  21. *
  22. * Modification History:
  23. * <04 Jan 1994> - began coding (jdm)
  24. * <06 Jan 1994> - announced version 0.1B on pasture.ecn.purdue.edu (jdm)
  25. * <21 Jan 1994> - changed wording on help screen. Revised to 0.2B (jdm)
  26. * <24 Jan 1994> - got rid of diagnostic messages. Revised to 0.3B (jdm)
  27. * ?? Revised to 0.4B (jdm)
  28. * <19 Dec 1994> - fixed bug in readsites, added html. Revised to 0.5B (jdm)
  29. * <02 Jan 1995> - cleaned Gmakefile, man page, html.
  30. * fixed memory error in bilinear and cubic 0.6B (jdm)
  31. * <25 Feb 1995> - cleaned 'gcc -Wall' warnings 0.7B (jdm)
  32. * <15 Jun 1995> - fixed pointer error for G_{col,row}_to_{easting,northing}.
  33. * 0.8B (jdm)
  34. * <13 Sep 2000> - released under GPL
  35. *
  36. */
  37. /* s.sample v 0.8B <15 Jun 1995>; Copyright (c) 1994-1995. James Darrell McCauley" */
  38. #include <stdio.h>
  39. #include <stdlib.h>
  40. #include <string.h>
  41. #include <math.h>
  42. #include <grass/gis.h>
  43. #include <grass/raster.h>
  44. #include <grass/glocale.h>
  45. #include <grass/dbmi.h>
  46. #include <grass/vector.h>
  47. int main(int argc, char **argv)
  48. {
  49. double scale, predicted, actual;
  50. INTERP_TYPE method = UNKNOWN;
  51. int fdrast; /* file descriptor for raster map is int */
  52. struct Cell_head window;
  53. struct GModule *module;
  54. struct Map_info In, Out;
  55. struct
  56. {
  57. struct Option *input, *output, *rast, *z, *column;
  58. } parm;
  59. struct
  60. {
  61. struct Flag *B, *C;
  62. } flag;
  63. int line, nlines;
  64. struct line_pnts *Points;
  65. struct line_cats *Cats;
  66. /* Attributes */
  67. int field, nrecords;
  68. int ctype;
  69. struct field_info *Fi;
  70. dbDriver *Driver;
  71. dbCatValArray cvarr;
  72. char buf[2000];
  73. dbString sql;
  74. G_gisinit(argv[0]);
  75. module = G_define_module();
  76. G_add_keyword(_("vector"));
  77. module->description =
  78. _("Samples a raster map at vector point locations.");
  79. parm.input = G_define_standard_option(G_OPT_V_INPUT);
  80. parm.input->description = _("Vector map defining sample points");
  81. parm.column = G_define_option();
  82. parm.column->key = "column";
  83. parm.column->type = TYPE_STRING;
  84. parm.column->required = YES;
  85. parm.column->description =
  86. _("Vector map attribute column to use for comparison");
  87. parm.output = G_define_standard_option(G_OPT_V_OUTPUT);
  88. parm.output->description = _("Vector map to store differences");
  89. parm.rast = G_define_standard_option(G_OPT_R_INPUT);
  90. parm.rast->key = "raster";
  91. parm.rast->description = _("Raster map to be sampled");
  92. parm.z = G_define_option();
  93. parm.z->key = "z";
  94. parm.z->type = TYPE_DOUBLE;
  95. parm.z->required = NO;
  96. parm.z->answer = "1.0";
  97. parm.z->description =
  98. _("Option scaling factor for values read from raster map. "
  99. "Sampled values will be multiplied by this factor");
  100. flag.B = G_define_flag();
  101. flag.B->key = 'b';
  102. flag.B->description =
  103. _("Bilinear interpolation (default is nearest neighbor)");
  104. flag.C = G_define_flag();
  105. flag.C->key = 'c';
  106. flag.C->description =
  107. _("Cubic convolution interpolation (default is nearest neighbor)");
  108. if (G_parser(argc, argv))
  109. exit(EXIT_FAILURE);
  110. sscanf(parm.z->answer, "%lf", &scale);
  111. if (flag.B->answer || flag.C->answer) {
  112. if (flag.C->answer)
  113. method = CUBIC;
  114. if (flag.B->answer)
  115. method = BILINEAR;
  116. if (flag.B->answer && flag.C->answer)
  117. G_fatal_error(_("Flags -b & -c are mutually exclusive. Choose only one."));
  118. }
  119. else {
  120. method = NEAREST;
  121. }
  122. G_get_window(&window);
  123. /* Open input */
  124. Vect_set_open_level(2);
  125. Vect_open_old(&In, parm.input->answer, "");
  126. if ((fdrast = Rast_open_old(parm.rast->answer, "")) < 0)
  127. G_fatal_error(_("Unable to open raster map <%s>"), parm.rast->answer);
  128. /* Read attributes */
  129. field = 1;
  130. Fi = Vect_get_field(&In, field);
  131. if (Fi == NULL)
  132. G_fatal_error(_("Database connection not defined for layer %d"),
  133. field);
  134. Driver = db_start_driver_open_database(Fi->driver, Fi->database);
  135. if (Driver == NULL)
  136. G_fatal_error("Unable to open database <%s> by driver <%s>",
  137. Fi->database, Fi->driver);
  138. nrecords = db_select_CatValArray(Driver, Fi->table, Fi->key,
  139. parm.column->answer, NULL, &cvarr);
  140. G_debug(3, "nrecords = %d", nrecords);
  141. ctype = cvarr.ctype;
  142. if (ctype != DB_C_TYPE_INT && ctype != DB_C_TYPE_DOUBLE)
  143. G_fatal_error(_("Column type <%s> not supported (must be integer or double precision)"),
  144. db_sqltype_name(ctype));
  145. if (nrecords < 0)
  146. G_fatal_error(_("Unable to select data from table"));
  147. G_message(_("%d records selected from table"), nrecords);
  148. db_close_database_shutdown_driver(Driver);
  149. /* Open output */
  150. Vect_open_new(&Out, parm.output->answer, 0);
  151. Vect_hist_copy(&In, &Out);
  152. Vect_hist_command(&Out);
  153. /* Create table */
  154. db_init_string(&sql);
  155. Fi = Vect_default_field_info(&Out, 1, NULL, GV_1TABLE);
  156. Vect_map_add_dblink(&Out, Fi->number, Fi->name, Fi->table, Fi->key,
  157. Fi->database, Fi->driver);
  158. Driver =
  159. db_start_driver_open_database(Fi->driver,
  160. Vect_subst_var(Fi->database, &Out));
  161. if (Driver == NULL)
  162. G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
  163. Fi->database, Fi->driver);
  164. sprintf(buf,
  165. "create table %s ( cat integer, pnt_val double precision, rast_val double precision, "
  166. "diff double precision)", Fi->table);
  167. db_set_string(&sql, buf);
  168. if (db_execute_immediate(Driver, &sql) != DB_OK)
  169. G_fatal_error(_("Unable to create table <%s>"), db_get_string(&sql));
  170. if (db_create_index2(Driver, Fi->table, Fi->key) != DB_OK)
  171. G_warning(_("Cannot create index"));
  172. if (db_grant_on_table
  173. (Driver, Fi->table, DB_PRIV_SELECT, DB_GROUP | DB_PUBLIC) != DB_OK)
  174. G_fatal_error(_("Unable to grant privileges on table <%s>"),
  175. Fi->table);
  176. G_message(_("Checking vector points..."));
  177. Points = Vect_new_line_struct();
  178. Cats = Vect_new_cats_struct();
  179. nlines = Vect_get_num_lines(&In);
  180. for (line = 1; line <= nlines; line++) {
  181. int type, cat, ret, cval;
  182. double dval;
  183. G_debug(3, "line = %d", line);
  184. type = Vect_read_line(&In, Points, Cats, line);
  185. if (!(type & GV_POINT))
  186. continue;
  187. Vect_cat_get(Cats, 1, &cat);
  188. G_debug(4, "cat = %d", cat);
  189. /* find actual value */
  190. if (ctype == DB_C_TYPE_INT) {
  191. ret = db_CatValArray_get_value_int(&cvarr, cat, &cval);
  192. if (ret != DB_OK)
  193. G_warning(_("No record for category %d in table <%s>"), cat,
  194. Fi->table);
  195. actual = cval;
  196. }
  197. else if (ctype == DB_C_TYPE_DOUBLE) {
  198. ret = db_CatValArray_get_value_double(&cvarr, cat, &dval);
  199. if (ret != DB_OK)
  200. G_warning(_("No record for category %d in table <%s>"), cat,
  201. Fi->table);
  202. actual = dval;
  203. }
  204. else {
  205. G_fatal_error(_("Column type not supported"));
  206. }
  207. G_debug(4, "actual = %e", actual);
  208. /* find predicted value */
  209. predicted = Rast_get_sample(fdrast, &window, NULL, Points->y[0],
  210. Points->x[0], 0, method);
  211. if (Rast_is_d_null_value(&predicted))
  212. continue;
  213. predicted *= scale;
  214. G_debug(4, "predicted = %e", predicted);
  215. Vect_reset_cats(Cats);
  216. Vect_cat_set(Cats, 1, cat);
  217. sprintf(buf, "insert into %s values ( %d, %e, %e, %e )",
  218. Fi->table, cat, actual, predicted, predicted - actual);
  219. db_set_string(&sql, buf);
  220. if (db_execute_immediate(Driver, &sql) != DB_OK)
  221. G_fatal_error(_("Unable to insert row: %s"), db_get_string(&sql));
  222. Vect_write_line(&Out, GV_POINT, Points, Cats);
  223. }
  224. db_close_database_shutdown_driver(Driver);
  225. Rast_close(fdrast);
  226. Vect_close(&In);
  227. Vect_build(&Out);
  228. Vect_close(&Out);
  229. exit(EXIT_SUCCESS);
  230. }