00001
00021 #include <string.h>
00022 #include <unistd.h>
00023 #include <math.h>
00024 #include <grass/gis.h>
00025 #include <grass/glocale.h>
00026
00027
00028 static double scancatlabel(const char *);
00029 static void raster_row_error(const struct Cell_head *window, double north,
00030 double east);
00031
00032
00050 DCELL G_get_raster_sample(
00051 int fd,
00052 const struct Cell_head *window,
00053 struct Categories *cats,
00054 double north, double east,
00055 int usedesc, INTERP_TYPE itype)
00056 {
00057 double retval;
00058
00059 switch (itype) {
00060 case NEAREST:
00061 retval = G_get_raster_sample_nearest(fd, window, cats, north, east, usedesc);
00062 break;
00063 case BILINEAR:
00064 retval = G_get_raster_sample_bilinear(fd, window, cats, north, east, usedesc);
00065 break;
00066 case CUBIC:
00067 retval = G_get_raster_sample_cubic(fd, window, cats, north, east, usedesc);
00068 break;
00069 default:
00070 G_fatal_error("G_get_raster_sample: %s",
00071 _("Unknown interpolation type"));
00072 }
00073
00074 return retval;
00075 }
00076
00077
00078 DCELL G_get_raster_sample_nearest(
00079 int fd,
00080 const struct Cell_head *window,
00081 struct Categories *cats,
00082 double north, double east, int usedesc)
00083 {
00084 int row, col;
00085 DCELL result;
00086 DCELL *maprow = G_allocate_d_raster_buf();
00087
00088
00089 row = (int)floor(G_northing_to_row(north, window));
00090 col = (int)floor(G_easting_to_col(east, window));
00091
00092 if (row < 0 || row >= G_window_rows() ||
00093 col < 0 || col >= G_window_cols()) {
00094 G_set_d_null_value(&result, 1);
00095 goto done;
00096 }
00097
00098 if (G_get_d_raster_row(fd, maprow, row) < 0)
00099 raster_row_error(window, north, east);
00100
00101 if (G_is_d_null_value(&maprow[col])) {
00102 G_set_d_null_value(&result, 1);
00103 goto done;
00104 }
00105
00106 if (usedesc) {
00107 char *buf = G_get_cat(maprow[col], cats);
00108
00109 G_squeeze(buf);
00110 result = scancatlabel(buf);
00111 }
00112 else
00113 result = maprow[col];
00114
00115 done:
00116 G_free(maprow);
00117
00118 return result;
00119 }
00120
00121
00122 DCELL G_get_raster_sample_bilinear(
00123 int fd,
00124 const struct Cell_head *window,
00125 struct Categories *cats,
00126 double north, double east, int usedesc)
00127 {
00128 int row, col;
00129 double grid[2][2];
00130 DCELL *arow = G_allocate_d_raster_buf();
00131 DCELL *brow = G_allocate_d_raster_buf();
00132 double frow, fcol, trow, tcol;
00133 DCELL result;
00134
00135 frow = G_northing_to_row(north, window);
00136 fcol = G_easting_to_col(east, window);
00137
00138
00139 row = (int)floor(frow - 0.5);
00140 col = (int)floor(fcol - 0.5);
00141
00142 trow = frow - row - 0.5;
00143 tcol = fcol - col - 0.5;
00144
00145 if (row < 0 || row + 1 >= G_window_rows() ||
00146 col < 0 || col + 1 >= G_window_cols()) {
00147 G_set_d_null_value(&result, 1);
00148 goto done;
00149 }
00150
00151 if (G_get_d_raster_row(fd, arow, row) < 0)
00152 raster_row_error(window, north, east);
00153 if (G_get_d_raster_row(fd, brow, row + 1) < 0)
00154 raster_row_error(window, north, east);
00155
00156 if (G_is_d_null_value(&arow[col]) || G_is_d_null_value(&arow[col + 1]) ||
00157 G_is_d_null_value(&brow[col]) || G_is_d_null_value(&brow[col + 1])) {
00158 G_set_d_null_value(&result, 1);
00159 goto done;
00160 }
00161
00162
00163
00164
00165
00166
00167
00168 if (usedesc) {
00169 char *buf;
00170
00171 G_squeeze(buf = G_get_cat((int)arow[col], cats));
00172 grid[0][0] = scancatlabel(buf);
00173 G_squeeze(buf = G_get_cat((int)arow[col + 1], cats));
00174 grid[0][1] = scancatlabel(buf);
00175 G_squeeze(buf = G_get_cat((int)brow[col], cats));
00176 grid[1][0] = scancatlabel(buf);
00177 G_squeeze(buf = G_get_cat((int)brow[col + 1], cats));
00178 grid[1][1] = scancatlabel(buf);
00179 }
00180 else {
00181 grid[0][0] = arow[col];
00182 grid[0][1] = arow[col + 1];
00183 grid[1][0] = brow[col];
00184 grid[1][1] = brow[col + 1];
00185 }
00186
00187 result = G_interp_bilinear(tcol, trow,
00188 grid[0][0], grid[0][1], grid[1][0], grid[1][1]);
00189
00190 done:
00191 G_free(arow);
00192 G_free(brow);
00193
00194 return result;
00195 }
00196
00197 DCELL G_get_raster_sample_cubic(
00198 int fd,
00199 const struct Cell_head *window,
00200 struct Categories *cats,
00201 double north, double east, int usedesc)
00202 {
00203 int i, j, row, col;
00204 double grid[4][4];
00205 DCELL *rows[4];
00206 double frow, fcol, trow, tcol;
00207 DCELL result;
00208
00209 for (i = 0; i < 4; i++)
00210 rows[i] = G_allocate_d_raster_buf();
00211
00212 frow = G_northing_to_row(north, window);
00213 fcol = G_easting_to_col(east, window);
00214
00215
00216 row = (int)floor(frow - 1.5);
00217 col = (int)floor(fcol - 1.5);
00218
00219 trow = frow - row - 1.5;
00220 tcol = fcol - col - 1.5;
00221
00222 if (row < 0 || row + 3 >= G_window_rows() ||
00223 col < 0 || col + 3 >= G_window_cols()) {
00224 G_set_d_null_value(&result, 1);
00225 goto done;
00226 }
00227
00228 for (i = 0; i < 4; i++)
00229 if (G_get_d_raster_row(fd, rows[i], row + i) < 0)
00230 raster_row_error(window, north, east);
00231
00232 for (i = 0; i < 4; i++)
00233 for (j = 0; j < 4; j++)
00234 if (G_is_d_null_value(&rows[i][col + j])) {
00235 G_set_d_null_value(&result, 1);
00236 goto done;
00237 }
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247 if (usedesc) {
00248 char *buf;
00249
00250 for (i = 0; i < 4; i++) {
00251 for (j = 0; j < 4; j++) {
00252 G_squeeze(buf = G_get_cat(rows[i][col + j], cats));
00253 grid[i][j] = scancatlabel(buf);
00254 }
00255 }
00256 }
00257 else {
00258 for (i = 0; i < 4; i++)
00259 for (j = 0; j < 4; j++)
00260 grid[i][j] = rows[i][col + j];
00261 }
00262
00263 result = G_interp_bicubic(tcol, trow,
00264 grid[0][0], grid[0][1], grid[0][2], grid[0][3],
00265 grid[1][0], grid[1][1], grid[1][2], grid[1][3],
00266 grid[2][0], grid[2][1], grid[2][2], grid[2][3],
00267 grid[3][0], grid[3][1], grid[3][2], grid[3][3]);
00268
00269 done:
00270 for (i = 0; i < 4; i++)
00271 G_free(rows[i]);
00272
00273 return result;
00274 }
00275
00276
00277 static double scancatlabel(const char *str)
00278 {
00279 double val;
00280
00281 if (strcmp(str, "no data") != 0)
00282 sscanf(str, "%lf", &val);
00283 else {
00284 G_warning(_("\"no data\" label found; setting to zero"));
00285 val = 0.0;
00286 }
00287
00288 return val;
00289 }
00290
00291
00292 static void raster_row_error(const struct Cell_head *window, double north,
00293 double east)
00294 {
00295 G_debug(3, "DIAG: \tRegion is: n=%g s=%g e=%g w=%g",
00296 window->north, window->south, window->east, window->west);
00297 G_debug(3, " \tData point is north=%g east=%g", north, east);
00298
00299 G_fatal_error(_("Problem reading raster map"));
00300 }