GRASS Programmer's Manual
6.4.1(2011)
|
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 }