GRASS Programmer's Manual  6.4.2(2012)
resout2d.c
Go to the documentation of this file.
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 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines