GRASS Programmer's Manual
6.4.2(2012)
|
00001 00002 /*- 00003 * Written by H. Mitasova, I. Kosinovsky, D. Gerdes Summer 1993 00004 * University of Illinois 00005 * US Army Construction Engineering Research Lab 00006 * Copyright 1993, H. Mitasova (University of Illinois), 00007 * I. Kosinovsky, (USA-CERL), and D.Gerdes (USA-CERL) 00008 * 00009 * modified by McCauley in August 1995 00010 * modified by Mitasova in August 1995 00011 * 00012 */ 00013 00014 #define MULT 100000 00015 00016 #include <stdio.h> 00017 #include <math.h> 00018 #include <grass/gis.h> 00019 #include <grass/bitmap.h> 00020 #include <grass/linkm.h> 00021 00022 #include <grass/interpf.h> 00023 00024 00025 /* output cell maps for elevation, aspect, slope and curvatures */ 00026 00027 int IL_resample_output_2d(struct interp_params *params, double zmin, double zmax, /* min,max input z-values */ 00028 double zminac, double zmaxac, /* min,max interpolated values */ 00029 double c1min, double c1max, double c2min, double c2max, double gmin, double gmax, double ertot, /* total interplating func. error */ 00030 char *input, /* input file name */ 00031 double *dnorm, struct Cell_head *outhd, /* Region with desired resolution */ 00032 struct Cell_head *winhd, /* Current region */ 00033 char *smooth, int n_points) 00034 00035 /* 00036 * Creates output files as well as history files and color tables for 00037 * them. 00038 */ 00039 { 00040 FCELL *cell1; /* cell buffer */ 00041 int cf1 = 0, cf2 = 0, cf3 = 0, cf4 = 0, cf5 = 0, cf6 = 0; /* cell file descriptors */ 00042 int nrows, ncols; /* current region rows and columns */ 00043 int i; /* loop counter */ 00044 char *mapset; 00045 float dat1, dat2; 00046 struct Colors colors, colors2; 00047 double value1, value2; 00048 struct History hist, hist1, hist2, hist3, hist4, hist5; 00049 struct _Color_Rule_ *rule; 00050 char *maps, *type; 00051 int cond1, cond2; 00052 00053 cond2 = ((params->pcurv != NULL) || 00054 (params->tcurv != NULL) || (params->mcurv != NULL)); 00055 cond1 = ((params->slope != NULL) || (params->aspect != NULL) || cond2); 00056 00057 /* change region to output cell file region */ 00058 fprintf(stderr, 00059 "Temporarily changing the region to desired resolution...\n"); 00060 if (G_set_window(outhd) < 0) { 00061 fprintf(stderr, "Cannot set region to output region!\n"); 00062 return -1; 00063 } 00064 mapset = G_mapset(); 00065 00066 cell1 = G_allocate_f_raster_buf(); 00067 00068 if (params->elev != NULL) { 00069 cf1 = G_open_fp_cell_new(params->elev); 00070 if (cf1 < 0) { 00071 fprintf(stderr, "unable to create raster map %s\n", params->elev); 00072 return -1; 00073 } 00074 } 00075 00076 if (params->slope != NULL) { 00077 cf2 = G_open_fp_cell_new(params->slope); 00078 if (cf2 < 0) { 00079 fprintf(stderr, "unable to create raster map %s\n", 00080 params->slope); 00081 return -1; 00082 } 00083 } 00084 00085 if (params->aspect != NULL) { 00086 cf3 = G_open_fp_cell_new(params->aspect); 00087 if (cf3 < 0) { 00088 fprintf(stderr, "unable to create raster map %s\n", 00089 params->aspect); 00090 return -1; 00091 } 00092 } 00093 00094 if (params->pcurv != NULL) { 00095 cf4 = G_open_fp_cell_new(params->pcurv); 00096 if (cf4 < 0) { 00097 fprintf(stderr, "unable to create raster map %s\n", 00098 params->pcurv); 00099 return -1; 00100 } 00101 } 00102 00103 if (params->tcurv != NULL) { 00104 cf5 = G_open_fp_cell_new(params->tcurv); 00105 if (cf5 < 0) { 00106 fprintf(stderr, "unable to create raster map %s\n", 00107 params->tcurv); 00108 return -1; 00109 } 00110 } 00111 00112 if (params->mcurv != NULL) { 00113 cf6 = G_open_fp_cell_new(params->mcurv); 00114 if (cf6 < 0) { 00115 fprintf(stderr, "unable to create raster map %s\n", 00116 params->mcurv); 00117 return -1; 00118 } 00119 } 00120 00121 nrows = outhd->rows; 00122 if (nrows != params->nsizr) { 00123 fprintf(stderr, "first change your rows number(%d) to %d!\n", 00124 nrows, params->nsizr); 00125 return -1; 00126 } 00127 00128 ncols = outhd->cols; 00129 if (ncols != params->nsizc) { 00130 fprintf(stderr, "first change your rows number(%d) to %d!\n", 00131 ncols, params->nsizc); 00132 return -1; 00133 } 00134 00135 if (params->elev != NULL) { 00136 fseek(params->Tmp_fd_z, 0L, 0); /* seek to the beginning */ 00137 for (i = 0; i < params->nsizr; i++) { 00138 /* seek to the right row */ 00139 if (fseek(params->Tmp_fd_z, (long) 00140 ((params->nsizr - 1 - 00141 i) * params->nsizc * sizeof(FCELL)), 0) == -1) { 00142 fprintf(stderr, "cannot fseek to the right spot\n"); 00143 return -1; 00144 } 00145 fread(cell1, sizeof(FCELL), params->nsizc, params->Tmp_fd_z); 00146 if (G_put_f_raster_row(cf1, cell1) < 0) { 00147 fprintf(stderr, "cannot write file\n"); 00148 return -1; 00149 } 00150 } 00151 } 00152 00153 if (params->slope != NULL) { 00154 fseek(params->Tmp_fd_dx, 0L, 0); /* seek to the beginning */ 00155 for (i = 0; i < params->nsizr; i++) { 00156 /* seek to the right row */ 00157 if (fseek(params->Tmp_fd_dx, (long) 00158 ((params->nsizr - 1 - 00159 i) * params->nsizc * sizeof(FCELL)), 0) == -1) { 00160 fprintf(stderr, "cannot fseek to the right spot\n"); 00161 return -1; 00162 } 00163 fread(cell1, sizeof(FCELL), params->nsizc, params->Tmp_fd_dx); 00164 /* 00165 * for (ii==0;ii<params->nsizc;ii++) { fprintf(stderr,"ii=%d ",ii); 00166 * fprintf(stderr,"%f ",cell1[ii]); } 00167 * fprintf(stderr,"params->nsizc=%d \n",params->nsizc); 00168 */ 00169 if (G_put_f_raster_row(cf2, cell1) < 0) { 00170 fprintf(stderr, "cannot write file\n"); 00171 return -1; 00172 } 00173 } 00174 } 00175 00176 if (params->aspect != NULL) { 00177 fseek(params->Tmp_fd_dy, 0L, 0); /* seek to the beginning */ 00178 for (i = 0; i < params->nsizr; i++) { 00179 /* seek to the right row */ 00180 if (fseek(params->Tmp_fd_dy, (long) 00181 ((params->nsizr - 1 - 00182 i) * params->nsizc * sizeof(FCELL)), 0) == -1) { 00183 fprintf(stderr, "cannot fseek to the right spot\n"); 00184 return -1; 00185 } 00186 fread(cell1, sizeof(FCELL), params->nsizc, params->Tmp_fd_dy); 00187 if (G_put_f_raster_row(cf3, cell1) < 0) { 00188 fprintf(stderr, "cannot write file\n"); 00189 return -1; 00190 } 00191 } 00192 } 00193 00194 if (params->pcurv != NULL) { 00195 fseek(params->Tmp_fd_xx, 0L, 0); /* seek to the beginning */ 00196 for (i = 0; i < params->nsizr; i++) { 00197 /* seek to the right row */ 00198 if (fseek(params->Tmp_fd_xx, (long) 00199 ((params->nsizr - 1 - 00200 i) * params->nsizc * sizeof(FCELL)), 0) == -1) { 00201 fprintf(stderr, "cannot fseek to the right spot\n"); 00202 return -1; 00203 } 00204 fread(cell1, sizeof(FCELL), params->nsizc, params->Tmp_fd_xx); 00205 if (G_put_f_raster_row(cf4, cell1) < 0) { 00206 fprintf(stderr, "cannot write file\n"); 00207 return -1; 00208 } 00209 } 00210 } 00211 00212 if (params->tcurv != NULL) { 00213 fseek(params->Tmp_fd_yy, 0L, 0); /* seek to the beginning */ 00214 for (i = 0; i < params->nsizr; i++) { 00215 /* seek to the right row */ 00216 if (fseek(params->Tmp_fd_yy, (long) 00217 ((params->nsizr - 1 - 00218 i) * params->nsizc * sizeof(FCELL)), 0) == -1) { 00219 fprintf(stderr, "cannot fseek to the right spot\n"); 00220 return -1; 00221 } 00222 fread(cell1, sizeof(FCELL), params->nsizc, params->Tmp_fd_yy); 00223 if (G_put_f_raster_row(cf5, cell1) < 0) { 00224 fprintf(stderr, "cannot write file\n"); 00225 return -1; 00226 } 00227 } 00228 } 00229 00230 if (params->mcurv != NULL) { 00231 fseek(params->Tmp_fd_xy, 0L, 0); /* seek to the beginning */ 00232 for (i = 0; i < params->nsizr; i++) { 00233 /* seek to the right row */ 00234 if (fseek(params->Tmp_fd_xy, (long) 00235 ((params->nsizr - 1 - 00236 i) * params->nsizc * sizeof(FCELL)), 0) == -1) { 00237 fprintf(stderr, "cannot fseek to the right spot\n"); 00238 return -1; 00239 } 00240 fread(cell1, sizeof(FCELL), params->nsizc, params->Tmp_fd_xy); 00241 if (G_put_f_raster_row(cf6, cell1) < 0) { 00242 fprintf(stderr, "cannot write file\n"); 00243 return -1; 00244 } 00245 } 00246 } 00247 00248 if (cf1) 00249 G_close_cell(cf1); 00250 if (cf2) 00251 G_close_cell(cf2); 00252 if (cf3) 00253 G_close_cell(cf3); 00254 if (cf4) 00255 G_close_cell(cf4); 00256 if (cf5) 00257 G_close_cell(cf5); 00258 if (cf6) 00259 G_close_cell(cf6); 00260 00261 /* write colormaps and history for output cell files */ 00262 /* colortable for elevations */ 00263 maps = G_find_file("cell", input, ""); 00264 00265 if (params->elev != NULL) { 00266 if (maps == NULL) { 00267 fprintf(stderr, "file [%s] not found\n", input); 00268 return -1; 00269 } 00270 G_init_colors(&colors2); 00271 /* 00272 * G_mark_colors_as_fp(&colors2); 00273 */ 00274 00275 if (G_read_colors(input, maps, &colors) >= 0) { 00276 if (colors.modular.rules) { 00277 rule = colors.modular.rules; 00278 00279 while (rule->next) 00280 rule = rule->next; 00281 00282 for (; rule; rule = rule->prev) { 00283 value1 = rule->low.value * params->zmult; 00284 value2 = rule->high.value * params->zmult; 00285 G_add_modular_d_raster_color_rule(&value1, rule->low.red, 00286 rule->low.grn, 00287 rule->low.blu, &value2, 00288 rule->high.red, 00289 rule->high.grn, 00290 rule->high.blu, 00291 &colors2); 00292 } 00293 } 00294 00295 if (colors.fixed.rules) { 00296 rule = colors.fixed.rules; 00297 00298 while (rule->next) 00299 rule = rule->next; 00300 00301 for (; rule; rule = rule->prev) { 00302 value1 = rule->low.value * params->zmult; 00303 value2 = rule->high.value * params->zmult; 00304 G_add_d_raster_color_rule(&value1, rule->low.red, 00305 rule->low.grn, rule->low.blu, 00306 &value2, rule->high.red, 00307 rule->high.grn, rule->high.blu, 00308 &colors2); 00309 } 00310 } 00311 00312 maps = NULL; 00313 maps = G_find_file("cell", params->elev, ""); 00314 if (maps == NULL) { 00315 fprintf(stderr, "file [%s] not found\n", params->elev); 00316 return -1; 00317 } 00318 00319 if (G_write_colors(params->elev, maps, &colors2) < 0) { 00320 fprintf(stderr, "Cannot write color table\n"); 00321 return -1; 00322 } 00323 G_quantize_fp_map_range(params->elev, mapset, 00324 zminac - 0.5, zmaxac + 0.5, 00325 (CELL) (zminac - 0.5), 00326 (CELL) (zmaxac + 0.5)); 00327 } 00328 else 00329 fprintf(stderr, 00330 "No color table for input file -- will not create color table\n"); 00331 } 00332 00333 /* colortable for slopes */ 00334 if (cond1 & (!params->deriv)) { 00335 G_init_colors(&colors); 00336 G_add_color_rule(0, 255, 255, 255, 2, 255, 255, 0, &colors); 00337 G_add_color_rule(2, 255, 255, 0, 5, 0, 255, 0, &colors); 00338 G_add_color_rule(5, 0, 255, 0, 10, 0, 255, 255, &colors); 00339 G_add_color_rule(10, 0, 255, 255, 15, 0, 0, 255, &colors); 00340 G_add_color_rule(15, 0, 0, 255, 30, 255, 0, 255, &colors); 00341 G_add_color_rule(30, 255, 0, 255, 50, 255, 0, 0, &colors); 00342 G_add_color_rule(50, 255, 0, 0, 90, 0, 0, 0, &colors); 00343 00344 if (params->slope != NULL) { 00345 maps = NULL; 00346 maps = G_find_file("cell", params->slope, ""); 00347 if (maps == NULL) { 00348 fprintf(stderr, "file [%s] not found\n", params->slope); 00349 return -1; 00350 } 00351 G_write_colors(params->slope, maps, &colors); 00352 G_quantize_fp_map_range(params->slope, mapset, 0., 90., 0, 90); 00353 00354 type = "raster"; 00355 G_short_history(params->slope, type, &hist1); 00356 if (params->elev != NULL) 00357 sprintf(hist1.edhist[0], "The elevation map is %s", 00358 params->elev); 00359 00360 sprintf(hist1.datsrc_1, "raster map %s", input); 00361 hist1.edlinecnt = 1; 00362 00363 G_write_history(params->slope, &hist1); 00364 } 00365 00366 /* colortable for aspect */ 00367 G_init_colors(&colors); 00368 G_add_color_rule(0, 255, 255, 255, 0, 255, 255, 255, &colors); 00369 G_add_color_rule(1, 255, 255, 0, 90, 0, 255, 0, &colors); 00370 G_add_color_rule(90, 0, 255, 0, 180, 0, 255, 255, &colors); 00371 G_add_color_rule(180, 0, 255, 255, 270, 255, 0, 0, &colors); 00372 G_add_color_rule(270, 255, 0, 0, 360, 255, 255, 0, &colors); 00373 00374 if (params->aspect != NULL) { 00375 maps = NULL; 00376 maps = G_find_file("cell", params->aspect, ""); 00377 if (maps == NULL) { 00378 fprintf(stderr, "file [%s] not found\n", params->aspect); 00379 return -1; 00380 } 00381 G_write_colors(params->aspect, maps, &colors); 00382 G_quantize_fp_map_range(params->aspect, mapset, 0., 360., 0, 360); 00383 00384 type = "raster"; 00385 G_short_history(params->aspect, type, &hist2); 00386 if (params->elev != NULL) 00387 sprintf(hist2.edhist[0], "The elevation map is %s", 00388 params->elev); 00389 00390 sprintf(hist2.datsrc_1, "raster map %s", input); 00391 hist2.edlinecnt = 1; 00392 00393 G_write_history(params->aspect, &hist2); 00394 } 00395 00396 /* colortable for curvatures */ 00397 if (cond2) { 00398 G_init_colors(&colors); 00399 00400 dat1 = (FCELL) amin1(c1min, c2min); 00401 dat2 = (FCELL) - 0.01; 00402 00403 G_add_f_raster_color_rule(&dat1, 50, 0, 155, 00404 &dat2, 0, 0, 255, &colors); 00405 dat1 = dat2; 00406 dat2 = (FCELL) - 0.001; 00407 G_add_f_raster_color_rule(&dat1, 0, 0, 255, 00408 &dat2, 0, 127, 255, &colors); 00409 dat1 = dat2; 00410 dat2 = (FCELL) - 0.00001; 00411 G_add_f_raster_color_rule(&dat1, 0, 127, 255, 00412 &dat2, 0, 255, 255, &colors); 00413 dat1 = dat2; 00414 dat2 = (FCELL) 0.00; 00415 G_add_f_raster_color_rule(&dat1, 0, 255, 255, 00416 &dat2, 200, 255, 200, &colors); 00417 dat1 = dat2; 00418 dat2 = (FCELL) 0.00001; 00419 G_add_f_raster_color_rule(&dat1, 200, 255, 200, 00420 &dat2, 255, 255, 0, &colors); 00421 dat1 = dat2; 00422 dat2 = (FCELL) 0.001; 00423 G_add_f_raster_color_rule(&dat1, 255, 255, 0, 00424 &dat2, 255, 127, 0, &colors); 00425 dat1 = dat2; 00426 dat2 = (FCELL) 0.01; 00427 G_add_f_raster_color_rule(&dat1, 255, 127, 0, 00428 &dat2, 255, 0, 0, &colors); 00429 dat1 = dat2; 00430 dat2 = (FCELL) amax1(c1max, c2max); 00431 G_add_f_raster_color_rule(&dat1, 255, 0, 0, 00432 &dat2, 155, 0, 20, &colors); 00433 maps = NULL; 00434 if (params->pcurv != NULL) { 00435 maps = G_find_file("cell", params->pcurv, ""); 00436 if (maps == NULL) { 00437 fprintf(stderr, "file [%s] not found\n", params->pcurv); 00438 return -1; 00439 } 00440 G_write_colors(params->pcurv, maps, &colors); 00441 00442 fprintf(stderr, "color map written\n"); 00443 00444 G_quantize_fp_map_range(params->pcurv, mapset, 00445 dat1, dat2, 00446 (CELL) (dat1 * MULT), 00447 (CELL) (dat2 * MULT)); 00448 type = "raster"; 00449 G_short_history(params->pcurv, type, &hist3); 00450 if (params->elev != NULL) 00451 sprintf(hist3.edhist[0], "The elevation map is %s", 00452 params->elev); 00453 00454 sprintf(hist3.datsrc_1, "raster map %s", input); 00455 hist3.edlinecnt = 1; 00456 00457 G_write_history(params->pcurv, &hist3); 00458 } 00459 00460 if (params->tcurv != NULL) { 00461 maps = NULL; 00462 maps = G_find_file("cell", params->tcurv, ""); 00463 if (maps == NULL) { 00464 fprintf(stderr, "file [%s] not found\n", params->tcurv); 00465 return -1; 00466 } 00467 G_write_colors(params->tcurv, maps, &colors); 00468 G_quantize_fp_map_range(params->tcurv, mapset, 00469 dat1, dat2, (CELL) (dat1 * MULT), 00470 (CELL) (dat2 * MULT)); 00471 00472 type = "raster"; 00473 G_short_history(params->tcurv, type, &hist4); 00474 if (params->elev != NULL) 00475 sprintf(hist4.edhist[0], "The elevation map is %s", 00476 params->elev); 00477 00478 sprintf(hist4.datsrc_1, "raster map %s", input); 00479 hist4.edlinecnt = 1; 00480 00481 G_write_history(params->tcurv, &hist4); 00482 } 00483 00484 if (params->mcurv != NULL) { 00485 maps = NULL; 00486 maps = G_find_file("cell", params->mcurv, ""); 00487 if (maps == NULL) { 00488 fprintf(stderr, "file [%s] not found\n", params->mcurv); 00489 return -1; 00490 } 00491 G_write_colors(params->mcurv, maps, &colors); 00492 G_quantize_fp_map_range(params->mcurv, mapset, 00493 dat1, dat2, 00494 (CELL) (dat1 * MULT), 00495 (CELL) (dat2 * MULT)); 00496 00497 type = "raster"; 00498 G_short_history(params->mcurv, type, &hist5); 00499 if (params->elev != NULL) 00500 sprintf(hist5.edhist[0], "The elevation map is %s", 00501 params->elev); 00502 00503 sprintf(hist5.datsrc_1, "raster map %s", input); 00504 hist5.edlinecnt = 1; 00505 00506 G_write_history(params->mcurv, &hist5); 00507 } 00508 } 00509 } 00510 00511 if (params->elev != NULL) { 00512 maps = G_find_file("cell", params->elev, ""); 00513 if (maps == NULL) { 00514 fprintf(stderr, "file [%s] not found \n", params->elev); 00515 return -1; 00516 } 00517 G_short_history(params->elev, "raster", &hist); 00518 00519 if (smooth != NULL) 00520 sprintf(hist.edhist[0], "tension=%f, smoothing=%s", 00521 params->fi * 1000. / (*dnorm), smooth); 00522 else 00523 sprintf(hist.edhist[0], "tension=%f", 00524 params->fi * 1000. / (*dnorm)); 00525 sprintf(hist.edhist[1], "dnorm=%f, zmult=%f", *dnorm, params->zmult); 00526 sprintf(hist.edhist[2], "KMAX=%d, KMIN=%d, errtotal=%f", params->kmax, 00527 params->kmin, sqrt(ertot / n_points)); 00528 sprintf(hist.edhist[3], "zmin_data=%f, zmax_data=%f", zmin, zmax); 00529 sprintf(hist.edhist[4], "zmin_int=%f, zmax_int=%f", zminac, zmaxac); 00530 00531 sprintf(hist.datsrc_1, "raster map %s", input); 00532 00533 hist.edlinecnt = 5; 00534 00535 G_write_history(params->elev, &hist); 00536 } 00537 00538 /* change region to initial region */ 00539 fprintf(stderr, "Changing the region back to initial...\n"); 00540 if (G_set_window(winhd) < 0) { 00541 fprintf(stderr, "Cannot set region to back to initial region!\n"); 00542 return -1; 00543 } 00544 00545 return 1; 00546 }