• Main Page
  • Related Pages
  • Data Structures
  • Files
  • File List
  • Globals

sample.c

Go to the documentation of this file.
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 /* prototypes */
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     /* convert northing and easting to row and col, resp */
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     /* convert northing and easting to row and col, resp */
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      * now were ready to do bilinear interpolation over
00164      * arow[col], arow[col+1],
00165      * brow[col], brow[col+1]
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     /* convert northing and easting to row and col, resp */
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      * now were ready to do cubic interpolation over
00241      * arow[col], arow[col+1], arow[col+2], arow[col+3],
00242      * brow[col], brow[col+1], brow[col+2], brow[col+3],
00243      * crow[col], crow[col+1], crow[col+2], crow[col+3],
00244      * drow[col], drow[col+1], drow[col+2], drow[col+3],
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 }

Generated on Wed Oct 13 2010 12:09:30 for GRASS Programmer's Manual by  doxygen 1.7.1