00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #include <string.h>
00013 #include <unistd.h>
00014 #include <math.h>
00015 #include <grass/gis.h>
00016 #include <grass/glocale.h>
00017
00018
00019
00020 static double raster_sample_nearest (int fd,
00021 struct Cell_head *window,
00022 struct Categories *cats,
00023 double north, double east, int usedesc);
00024 static double raster_sample_bilinear (int fd,
00025 struct Cell_head *window,
00026 struct Categories *cats,
00027 double north, double east, int usedesc);
00028 static double raster_sample_cubic (int fd,
00029 struct Cell_head *window,
00030 struct Categories *cats,
00031 double north, double east, int usedesc);
00032 static double scancatlabel (char *);
00033 static void raster_row_error(struct Cell_head *window, double north, double east);
00034
00035
00054 double G_get_raster_sample (int fd,
00055 struct Cell_head *window,
00056 struct Categories *cats,
00057 double north, double east,
00058 int usedesc, INTERP_TYPE itype)
00059 {
00060 double retval;
00061
00062 switch (itype)
00063 {
00064 case NEAREST:
00065 retval = raster_sample_nearest(fd, window, cats, north, east, usedesc);
00066 break;
00067 case BILINEAR:
00068 retval = raster_sample_bilinear(fd, window, cats, north, east, usedesc);
00069 break;
00070 case CUBIC:
00071 retval = raster_sample_cubic(fd, window, cats, north, east, usedesc);
00072 break;
00073 default:
00074 G_fatal_error(_("unknown interpolation type"));
00075 }
00076
00077 return retval;
00078 }
00079
00080
00081 static double raster_sample_nearest (int fd,
00082 struct Cell_head *window,
00083 struct Categories *cats,
00084 double north, double east, int usedesc)
00085 {
00086 int row, col;
00087 double predicted;
00088 DCELL *maprow = G_allocate_d_raster_buf();
00089
00090
00091 row = (int) G_northing_to_row (north, window);
00092 col = (int) G_easting_to_col (east, window);
00093
00094 if (G_get_d_raster_row (fd, maprow, row) < 0)
00095 raster_row_error(window, north, east);
00096
00097 if (G_is_d_null_value(&(maprow[col])))
00098 {
00099 predicted = 0.0;
00100 }
00101 else if (usedesc)
00102 {
00103 char *buf = G_get_cat(maprow[col], cats);
00104
00105 G_squeeze(buf);
00106 predicted = scancatlabel(buf);
00107 }
00108 else
00109 {
00110 predicted = maprow[col];
00111 }
00112
00113 G_free(maprow);
00114
00115 return predicted;
00116 }
00117
00118
00119 static double raster_sample_bilinear (int fd,
00120 struct Cell_head *window,
00121 struct Categories *cats,
00122 double north, double east, int usedesc)
00123 {
00124 int i, row, col;
00125 double grid[2][2], tmp1, tmp2;
00126 DCELL *arow = G_allocate_d_raster_buf ();
00127 DCELL *brow = G_allocate_d_raster_buf ();
00128
00129
00130 row = (int) G_northing_to_row (north, window);
00131 col = (int) G_easting_to_col (east, window);
00132
00133 if (G_get_d_raster_row (fd, arow, row) < 0)
00134 raster_row_error(window, north, east);
00135
00136
00137
00138
00139
00140 if (row == 0)
00141 {
00142
00143 if (G_get_d_raster_row (fd, brow, row + 1) < 0)
00144 raster_row_error(window, north, east);
00145 }
00146 else if (row+1 == G_window_rows ())
00147 {
00148
00149
00150 for(i = 0; i < G_window_cols(); ++i)
00151 brow[i] = arow[i];
00152
00153 row--;
00154
00155 if (G_get_d_raster_row (fd, arow, row) < 0)
00156 raster_row_error(window, north, east);
00157 }
00158 else if (north - G_row_to_northing ((double) row + 0.5, window) > 0)
00159 {
00160
00161
00162 for(i = 0; i < G_window_cols(); ++i)
00163 brow[i] = arow[i];
00164
00165 row--;
00166
00167 if (G_get_d_raster_row (fd, arow, row) < 0)
00168 raster_row_error(window, north, east);
00169 }
00170 else
00171 {
00172
00173 if (G_get_d_raster_row (fd, brow, row+1) < 0)
00174 raster_row_error(window, north, east);
00175 }
00176
00177
00178
00179
00180
00181 if (col+1 == G_window_cols())
00182 col--;
00183 else if (east - G_col_to_easting((double) col + 0.5, window) < 0)
00184 col--;
00185
00186
00187
00188
00189
00190
00191
00192 if (usedesc)
00193 {
00194 char *buf;
00195
00196 G_squeeze(buf = G_get_cat ((int)arow[col], cats));
00197 grid[0][0] = scancatlabel (buf);
00198 G_squeeze(buf = G_get_cat ((int)arow[col+1], cats));
00199 grid[0][1] = scancatlabel (buf);
00200 G_squeeze(buf = G_get_cat ((int)brow[col], cats));
00201 grid[1][0] = scancatlabel (buf);
00202 G_squeeze(buf = G_get_cat ((int)brow[col+1], cats));
00203 grid[1][1] = scancatlabel (buf);
00204 }
00205 else
00206 {
00207 grid[0][0] = (double) arow[col];
00208 grid[0][1] = (double) arow[col + 1];
00209 grid[1][0] = (double) brow[col];
00210 grid[1][1] = (double) brow[col + 1];
00211 }
00212
00213
00214 if (G_is_d_null_value(&(arow[col])))
00215 grid[0][0] = 0.0;
00216 if (G_is_d_null_value(&(arow[col+1])))
00217 grid[0][1] = 0.0;
00218 if (G_is_d_null_value(&(brow[col])))
00219 grid[1][0] = 0.0;
00220 if (G_is_d_null_value(&(brow[col+1])))
00221 grid[1][1] = 0.0;
00222
00223 east = fabs(G_col_to_easting((double)col, window) - east);
00224 while (east > window->ew_res)
00225 east -= window->ew_res;
00226
00227 north = fabs(G_row_to_northing((double)row, window) - north);
00228 while (north > window->ns_res)
00229 north -= window->ns_res;
00230
00231
00232 tmp1 = east * grid[0][1] + (window->ew_res - east) * grid[0][0];
00233 tmp1 /= window->ew_res;
00234 tmp2 = east * grid[1][1] + (window->ew_res - east) * grid[1][0];
00235 tmp2 /= window->ew_res;
00236 G_debug(3, "DIAG: r=%d c=%d t1=%g t2=%g\te=%g n=%g",
00237 row, col, tmp1, tmp2, east, north);
00238 G_debug(3, "DIAG: %g %g\n %g %g",
00239 grid[0][0], grid[0][1], grid[1][0], grid[1][1]);
00240
00241
00242
00243
00244
00245 G_free(arow);
00246 G_free(brow);
00247
00248 return (north * tmp2 + (window->ns_res - north) * tmp1) / window->ns_res;
00249 }
00250
00251
00252 static double raster_sample_cubic (int fd,
00253 struct Cell_head *window,
00254 struct Categories *cats,
00255 double north, double east, int usedesc)
00256 {
00257 int i, row, col;
00258 double grid[4][4], tmp[4];
00259 DCELL *arow = G_allocate_d_raster_buf();
00260 DCELL *brow = G_allocate_d_raster_buf();
00261 DCELL *crow = G_allocate_d_raster_buf();
00262 DCELL *drow = G_allocate_d_raster_buf();
00263
00264
00265 row = G_northing_to_row(north, window);
00266 col = G_easting_to_col(east, window);
00267
00268 if (G_get_d_raster_row(fd, arow, row) < 0)
00269 raster_row_error(window, north, east);
00270
00271
00272 if (row == 0)
00273 {
00274
00275 if (G_get_d_raster_row(fd, brow, row + 1) < 0)
00276 raster_row_error(window, north, east);
00277 if (G_get_d_raster_row(fd, crow, row + 2) < 0)
00278 raster_row_error(window, north, east);
00279 if (G_get_d_raster_row(fd, drow, row + 3) < 0)
00280 raster_row_error(window, north, east);
00281 }
00282 else if (row == 1)
00283 {
00284
00285 for (i = 0; i < G_window_cols(); ++i)
00286 brow[i] = arow[i];
00287
00288 if (G_get_d_raster_row(fd, arow, row - 1) < 0)
00289 raster_row_error(window, north, east);
00290 if (G_get_d_raster_row(fd, crow, row + 1) < 0)
00291 raster_row_error(window, north, east);
00292 if (G_get_d_raster_row(fd, drow, row + 2) < 0)
00293 raster_row_error(window, north, east);
00294
00295 row--;
00296 }
00297 else if (row + 1 == G_window_rows())
00298 {
00299
00300 for (i = 0; i < G_window_cols(); ++i)
00301 drow[i] = arow[i];
00302
00303 if (G_get_d_raster_row(fd, arow, row - 3) < 0)
00304 raster_row_error(window, north, east);
00305 if (G_get_d_raster_row(fd, brow, row - 2) < 0)
00306 raster_row_error(window, north, east);
00307 if (G_get_d_raster_row(fd, crow, row - 1) < 0)
00308 raster_row_error(window, north, east);
00309
00310 row -= 3;
00311 }
00312 else if (row + 2 == G_window_rows())
00313 {
00314
00315 for (i = 0; i < G_window_cols(); ++i)
00316 crow[i] = arow[i];
00317
00318 if (G_get_d_raster_row(fd, arow, row - 2) < 0)
00319 raster_row_error(window, north, east);
00320 if (G_get_d_raster_row(fd, brow, row - 1) < 0)
00321 raster_row_error(window, north, east);
00322 if (G_get_d_raster_row(fd, drow, row + 1) < 0)
00323 raster_row_error(window, north, east);
00324
00325 row -= 2;
00326 }
00327 else if (north - G_row_to_northing((double)row + 0.5, window) > 0)
00328 {
00329
00330
00331
00332
00333 for (i = 0; i < G_window_cols(); ++i)
00334 crow[i] = arow[i];
00335
00336 if (G_get_d_raster_row(fd, arow, row - 2) < 0)
00337 raster_row_error(window, north, east);
00338 if (G_get_d_raster_row(fd, brow, row - 1) < 0)
00339 raster_row_error(window, north, east);
00340 if (G_get_d_raster_row(fd, drow, row + 1) < 0)
00341 raster_row_error(window, north, east);
00342
00343 row -= 2;
00344 }
00345 else
00346 {
00347
00348
00349
00350
00351 for (i = 0; i < G_window_cols(); ++i)
00352 brow[i] = arow[i];
00353
00354 if (G_get_d_raster_row(fd, arow, row - 1) < 0)
00355 raster_row_error(window, north, east);
00356 if (G_get_d_raster_row(fd, crow, row + 1) < 0)
00357 raster_row_error(window, north, east);
00358 if (G_get_d_raster_row(fd, drow, row + 2) < 0)
00359 raster_row_error(window, north, east);
00360
00361 row--;
00362 }
00363
00364
00365
00366
00367
00368 if (col == 0 || col == 1)
00369 col = 0;
00370 else if (col + 1 == G_window_cols())
00371 col -= 3;
00372 else if (col + 2 == G_window_cols())
00373 col -= 2;
00374 else if (east - G_col_to_easting((double)col + 0.5, window) < 0)
00375
00376 col -= 2;
00377 else
00378 col--;
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388 if (usedesc)
00389 {
00390 char *buf;
00391
00392 for (i = 0; i < 4; ++i) {
00393 G_squeeze(buf = G_get_cat(arow[col + i], cats));
00394 grid[0][i] = scancatlabel(buf);
00395 G_squeeze(buf = G_get_cat(brow[col + i], cats));
00396 grid[1][i] = scancatlabel(buf);
00397 G_squeeze(buf = G_get_cat(crow[col + i], cats));
00398 grid[2][i] = scancatlabel(buf);
00399 G_squeeze(buf = G_get_cat(drow[col + i], cats));
00400 grid[3][i] = scancatlabel(buf);
00401 }
00402 }
00403 else
00404 {
00405 for (i = 0; i < 4; ++i) {
00406 grid[0][i] = (double)arow[col + i];
00407 grid[1][i] = (double)brow[col + i];
00408 grid[2][i] = (double)crow[col + i];
00409 grid[3][i] = (double)drow[col + i];
00410 }
00411 }
00412
00413
00414 for (i = 0; i < 4; i++) {
00415 if (G_is_d_null_value(&(arow[col + i])))
00416 grid[0][i] = 0.0;
00417 if (G_is_d_null_value(&(brow[col + i])))
00418 grid[1][i] = 0.0;
00419 if (G_is_d_null_value(&(crow[col + i])))
00420 grid[2][i] = 0.0;
00421 if (G_is_d_null_value(&(drow[col + i])))
00422 grid[3][i] = 0.0;
00423 }
00424
00425
00426 east = fabs(G_col_to_easting((double)col + 1, window) - east);
00427 while (east > window->ew_res)
00428 east -= window->ew_res;
00429 east /= window->ew_res;
00430
00431 north = fabs(G_row_to_northing((double)row + 1, window) - north);
00432 while (north > window->ns_res)
00433 north -= window->ns_res;
00434 north /= window->ns_res;
00435
00436
00437 for (i = 0; i < 4; ++i)
00438 tmp[i] = east * (east * (east * (grid[i][3] - grid[i][2] + grid[i][1] - grid[i][0])
00439 + (grid[i][2] - grid[i][3] - 2 * grid[i][1] + 2 * grid[i][0]))
00440 + (grid[i][2] - grid[i][0])) + grid[i][1];
00441
00442 #ifdef DEBUG
00443 for (i = 0; i < 4; ++i) {
00444 for (j = 0; j < 4; ++j)
00445 fprintf(stderr, "%g ", grid[i][j]);
00446 fprintf(stderr, "\n");
00447 }
00448
00449 G_message("DIAG: (%d,%d) 1=%3.2g 2=%3.2g 3=%3.2g 4=%3.2g\te=%g n=%g",
00450 row, col, tmp[0], tmp[1], tmp[2], tmp[3], east, north);
00451 #endif
00452
00453 G_free(arow);
00454 G_free(brow);
00455 G_free(crow);
00456 G_free(drow);
00457
00458
00459 return (north * (north * (north * (tmp[3] - tmp[2] + tmp[1] - tmp[0])
00460 + (tmp[2] - tmp[3] - 2 * tmp[1] + 2 * tmp[0]))
00461 + (tmp[2] - tmp[0])) + tmp[1]);
00462 }
00463
00464
00465 static double scancatlabel (char *str)
00466 {
00467 double val;
00468
00469 if (strcmp(str,"no data") != 0)
00470 sscanf(str, "%lf", &val);
00471 else
00472 {
00473 G_warning(_("\"no data\" label found; setting to zero"));
00474 val = 0.0;
00475 }
00476
00477 return val;
00478 }
00479
00480
00481 static void raster_row_error(struct Cell_head *window, double north, double east)
00482 {
00483 G_debug(3, "DIAG: \tRegion is: n=%g s=%g e=%g w=%g",
00484 window->north, window->south, window->east, window->west);
00485 G_debug(3, " \tData point is north=%g east=%g", north, east);
00486
00487 G_fatal_error(_("problem reading raster cell file"));
00488 }