resout2d.c

Go to the documentation of this file.

/*-
 * Written by H. Mitasova, I. Kosinovsky, D. Gerdes Summer 1993
 * University of Illinois
 * US Army Construction Engineering Research Lab
 * Copyright 1993, H. Mitasova (University of Illinois),
 * I. Kosinovsky, (USA-CERL), and D.Gerdes (USA-CERL)
 *
 * modified by McCauley in August 1995
 * modified by Mitasova in August 1995
 *
 */
 
#define MULT 100000
 
#include <stdio.h>
#include <math.h>
#include <errno.h>
 
#include <grass/gis.h>
#include <grass/raster.h>
#include <grass/bitmap.h>
#include <grass/linkm.h>
#include <grass/interpf.h>
#include <grass/glocale.h>
 
/* output cell maps for elevation, aspect, slope and curvatures */
 
static void do_history(const char *name, const char *input,
                       const struct interp_params *params)
{
    struct History hist;
 
    Rast_short_history(name, "raster", &hist);
    if (params->elev)
        Rast_append_format_history(&hist, "The elevation map is %s",
                                   params->elev);
 
    Rast_format_history(&hist, HIST_DATSRC_1, "raster map %s", input);
 
    Rast_write_history(name, &hist);
 
    Rast_free_history(&hist);
}
 
int IL_resample_output_2d(
    struct interp_params *params, double zmin,
    double zmax,                  /* min,max input z-values */
    double zminac, double zmaxac, /* min,max interpolated values */
    double c1min, double c1max, double c2min, double c2max, double gmin UNUSED,
    double gmax UNUSED, double ertot,       /* total interplating func. error */
    char *input,                            /* input file name */
    double *dnorm, struct Cell_head *outhd, /* Region with desired resolution */
    struct Cell_head *winhd,                /* Current region */
    char *smooth, int n_points)
/*
 * Creates output files as well as history files  and color tables for
 * them.
 */
{
    FCELL *cell1; /* cell buffer */
    int cf1 = 0, cf2 = 0, cf3 = 0, cf4 = 0, cf5 = 0,
        cf6 = 0;      /* cell file descriptors */
    int nrows, ncols; /* current region rows and columns */
    int i;            /* loop counter */
    const char *mapset;
    float dat1, dat2;
    struct Colors colors, colors2;
    double value1, value2;
    struct History hist;
    struct _Color_Rule_ *rule;
    const char *maps;
    int cond1, cond2;
    CELL val1, val2;
 
    cond2 = ((params->pcurv != NULL) || (params->tcurv != NULL) ||
             (params->mcurv != NULL));
    cond1 = ((params->slope != NULL) || (params->aspect != NULL) || cond2);
 
    /* change region to output cell file region */
    G_verbose_message(
        _("Temporarily changing the region to desired resolution..."));
    Rast_set_output_window(outhd);
    mapset = G_mapset();
 
    cell1 = Rast_allocate_f_output_buf();
 
    if (params->elev)
        cf1 = Rast_open_fp_new(params->elev);
 
    if (params->slope)
        cf2 = Rast_open_fp_new(params->slope);
 
    if (params->aspect)
        cf3 = Rast_open_fp_new(params->aspect);
 
    if (params->pcurv)
        cf4 = Rast_open_fp_new(params->pcurv);
 
    if (params->tcurv)
        cf5 = Rast_open_fp_new(params->tcurv);
 
    if (params->mcurv)
        cf6 = Rast_open_fp_new(params->mcurv);
 
    nrows = outhd->rows;
    if (nrows != params->nsizr) {
        G_warning(_("First change your rows number(%d) to %d"), nrows,
                  params->nsizr);
        return -1;
    }
 
    ncols = outhd->cols;
    if (ncols != params->nsizc) {
        G_warning(_("First change your columns number(%d) to %d"), ncols,
                  params->nsizr);
        return -1;
    }
 
    if (params->elev != NULL) {
        G_fseek(params->Tmp_fd_z, 0L, 0); /* seek to the beginning */
        for (i = 0; i < params->nsizr; i++) {
            /* seek to the right row */
            G_fseek(params->Tmp_fd_z,
                    (off_t)(params->nsizr - 1 - i) * params->nsizc *
                        sizeof(FCELL),
                    0);
            if (fread(cell1, sizeof(FCELL), params->nsizc, params->Tmp_fd_z) !=
                (size_t)params->nsizc)
                G_fatal_error(_("RST library temporary file reading error: %s"),
                              strerror(errno));
            Rast_put_f_row(cf1, cell1);
        }
    }
 
    if (params->slope != NULL) {
        G_fseek(params->Tmp_fd_dx, 0L, 0); /* seek to the beginning */
        for (i = 0; i < params->nsizr; i++) {
            /* seek to the right row */
            G_fseek(params->Tmp_fd_dx,
                    (off_t)(params->nsizr - 1 - i) * params->nsizc *
                        sizeof(FCELL),
                    0);
            if (fread(cell1, sizeof(FCELL), params->nsizc, params->Tmp_fd_dx) !=
                (size_t)params->nsizc)
                G_fatal_error(_("RST library temporary file reading error: %s"),
                              strerror(errno));
            Rast_put_f_row(cf2, cell1);
        }
    }
 
    if (params->aspect != NULL) {
        G_fseek(params->Tmp_fd_dy, 0L, 0); /* seek to the beginning */
        for (i = 0; i < params->nsizr; i++) {
            /* seek to the right row */
            G_fseek(params->Tmp_fd_dy,
                    (off_t)(params->nsizr - 1 - i) * params->nsizc *
                        sizeof(FCELL),
                    0);
            if (fread(cell1, sizeof(FCELL), params->nsizc, params->Tmp_fd_dy) !=
                (size_t)params->nsizc)
                G_fatal_error(_("RST library temporary file reading error: %s"),
                              strerror(errno));
            Rast_put_f_row(cf3, cell1);
        }
    }
 
    if (params->pcurv != NULL) {
        G_fseek(params->Tmp_fd_xx, 0L, 0); /* seek to the beginning */
        for (i = 0; i < params->nsizr; i++) {
            /* seek to the right row */
            G_fseek(params->Tmp_fd_xx,
                    (off_t)(params->nsizr - 1 - i) * params->nsizc *
                        sizeof(FCELL),
                    0);
            if (fread(cell1, sizeof(FCELL), params->nsizc, params->Tmp_fd_xx) !=
                (size_t)params->nsizc)
                G_fatal_error(_("RST library temporary file reading error: %s"),
                              strerror(errno));
            Rast_put_f_row(cf4, cell1);
        }
    }
 
    if (params->tcurv != NULL) {
        G_fseek(params->Tmp_fd_yy, 0L, 0); /* seek to the beginning */
        for (i = 0; i < params->nsizr; i++) {
            /* seek to the right row */
            G_fseek(params->Tmp_fd_yy,
                    (off_t)(params->nsizr - 1 - i) * params->nsizc *
                        sizeof(FCELL),
                    0);
            if (fread(cell1, sizeof(FCELL), params->nsizc, params->Tmp_fd_yy) !=
                (size_t)params->nsizc)
                G_fatal_error(_("RST library temporary file reading error: %s"),
                              strerror(errno));
            Rast_put_f_row(cf5, cell1);
        }
    }
 
    if (params->mcurv != NULL) {
        G_fseek(params->Tmp_fd_xy, 0L, 0); /* seek to the beginning */
        for (i = 0; i < params->nsizr; i++) {
            /* seek to the right row */
            G_fseek(params->Tmp_fd_xy,
                    (off_t)(params->nsizr - 1 - i) * params->nsizc *
                        sizeof(FCELL),
                    0);
            if (fread(cell1, sizeof(FCELL), params->nsizc, params->Tmp_fd_xy) !=
                (size_t)params->nsizc)
                G_fatal_error(_("RST library temporary file reading error: %s"),
                              strerror(errno));
            Rast_put_f_row(cf6, cell1);
        }
    }
 
    if (cf1)
        Rast_close(cf1);
    if (cf2)
        Rast_close(cf2);
    if (cf3)
        Rast_close(cf3);
    if (cf4)
        Rast_close(cf4);
    if (cf5)
        Rast_close(cf5);
    if (cf6)
        Rast_close(cf6);
 
    /* write colormaps and history for output cell files */
    /* colortable for elevations */
    maps = G_find_file("cell", input, "");
 
    if (params->elev != NULL) {
        if (maps == NULL) {
            G_warning(_("Raster map <%s> not found"), input);
            return -1;
        }
        Rast_init_colors(&colors2);
        /*
         * Rast_mark_colors_as_fp(&colors2);
         */
 
        if (Rast_read_colors(input, maps, &colors) >= 0) {
            if (colors.modular.rules) {
                rule = colors.modular.rules;
 
                while (rule->next)
                    rule = rule->next;
 
                for (; rule; rule = rule->prev) {
                    value1 = rule->low.value * params->zmult;
                    value2 = rule->high.value * params->zmult;
                    Rast_add_modular_d_color_rule(
                        &value1, rule->low.red, rule->low.grn, rule->low.blu,
                        &value2, rule->high.red, rule->high.grn, rule->high.blu,
                        &colors2);
                }
            }
 
            if (colors.fixed.rules) {
                rule = colors.fixed.rules;
 
                while (rule->next)
                    rule = rule->next;
 
                for (; rule; rule = rule->prev) {
                    value1 = rule->low.value * params->zmult;
                    value2 = rule->high.value * params->zmult;
                    Rast_add_d_color_rule(&value1, rule->low.red, rule->low.grn,
                                          rule->low.blu, &value2,
                                          rule->high.red, rule->high.grn,
                                          rule->high.blu, &colors2);
                }
            }
 
            maps = NULL;
            maps = G_find_file("cell", params->elev, "");
            if (maps == NULL) {
                G_warning(_("Raster map <%s> not found"), params->elev);
                return -1;
            }
 
            Rast_write_colors(params->elev, maps, &colors2);
            Rast_quantize_fp_map_range(params->elev, mapset, zminac - 0.5,
                                       zmaxac + 0.5, (CELL)(zminac - 0.5),
                                       (CELL)(zmaxac + 0.5));
        }
        else
            G_warning(_("No color table for input raster map -- will not "
                        "create color table"));
    }
 
    /* colortable for slopes */
    if (cond1 & (!params->deriv)) {
        Rast_init_colors(&colors);
        val1 = 0;
        val2 = 2;
        Rast_add_c_color_rule(&val1, 255, 255, 255, &val2, 255, 255, 0,
                              &colors);
        val1 = 2;
        val2 = 5;
        Rast_add_c_color_rule(&val1, 255, 255, 0, &val2, 0, 255, 0, &colors);
        val1 = 5;
        val2 = 10;
        Rast_add_c_color_rule(&val1, 0, 255, 0, &val2, 0, 255, 255, &colors);
        val1 = 10;
        val2 = 15;
        Rast_add_c_color_rule(&val1, 0, 255, 255, &val2, 0, 0, 255, &colors);
        val1 = 15;
        val2 = 30;
        Rast_add_c_color_rule(&val1, 0, 0, 255, &val2, 255, 0, 255, &colors);
        val1 = 30;
        val2 = 50;
        Rast_add_c_color_rule(&val1, 255, 0, 255, &val2, 255, 0, 0, &colors);
        val1 = 50;
        val2 = 90;
        Rast_add_c_color_rule(&val1, 255, 0, 0, &val2, 0, 0, 0, &colors);
 
        if (params->slope != NULL) {
            maps = NULL;
            maps = G_find_file("cell", params->slope, "");
            if (maps == NULL) {
                G_warning(_("Raster map <%s> not found"), params->slope);
                return -1;
            }
            Rast_write_colors(params->slope, maps, &colors);
            Rast_quantize_fp_map_range(params->slope, mapset, 0., 90., 0, 90);
 
            do_history(params->slope, input, params);
        }
 
        /* colortable for aspect */
        Rast_init_colors(&colors);
        val1 = 0;
        val2 = 0;
        Rast_add_c_color_rule(&val1, 255, 255, 255, &val2, 255, 255, 255,
                              &colors);
        val1 = 1;
        val2 = 90;
        Rast_add_c_color_rule(&val1, 255, 255, 0, &val2, 0, 255, 0, &colors);
        val1 = 90;
        val2 = 180;
        Rast_add_c_color_rule(&val1, 0, 255, 0, &val2, 0, 255, 255, &colors);
        val1 = 180;
        val2 = 270;
        Rast_add_c_color_rule(&val1, 0, 255, 255, &val2, 255, 0, 0, &colors);
        val1 = 270;
        val2 = 360;
        Rast_add_c_color_rule(&val1, 255, 0, 0, &val2, 255, 255, 0, &colors);
 
        if (params->aspect != NULL) {
            maps = NULL;
            maps = G_find_file("cell", params->aspect, "");
            if (maps == NULL) {
                G_warning(_("Raster map <%s> not found"), params->aspect);
                return -1;
            }
            Rast_write_colors(params->aspect, maps, &colors);
            Rast_quantize_fp_map_range(params->aspect, mapset, 0., 360., 0,
                                       360);
 
            do_history(params->aspect, input, params);
        }
 
        /* colortable for curvatures */
        if (cond2) {
            Rast_init_colors(&colors);
 
            dat1 = (FCELL)amin1(c1min, c2min);
            dat2 = (FCELL)-0.01;
 
            Rast_add_f_color_rule(&dat1, 50, 0, 155, &dat2, 0, 0, 255, &colors);
            dat1 = dat2;
            dat2 = (FCELL)-0.001;
            Rast_add_f_color_rule(&dat1, 0, 0, 255, &dat2, 0, 127, 255,
                                  &colors);
            dat1 = dat2;
            dat2 = (FCELL)-0.00001;
            Rast_add_f_color_rule(&dat1, 0, 127, 255, &dat2, 0, 255, 255,
                                  &colors);
            dat1 = dat2;
            dat2 = (FCELL)0.00;
            Rast_add_f_color_rule(&dat1, 0, 255, 255, &dat2, 200, 255, 200,
                                  &colors);
            dat1 = dat2;
            dat2 = (FCELL)0.00001;
            Rast_add_f_color_rule(&dat1, 200, 255, 200, &dat2, 255, 255, 0,
                                  &colors);
            dat1 = dat2;
            dat2 = (FCELL)0.001;
            Rast_add_f_color_rule(&dat1, 255, 255, 0, &dat2, 255, 127, 0,
                                  &colors);
            dat1 = dat2;
            dat2 = (FCELL)0.01;
            Rast_add_f_color_rule(&dat1, 255, 127, 0, &dat2, 255, 0, 0,
                                  &colors);
            dat1 = dat2;
            dat2 = (FCELL)amax1(c1max, c2max);
            Rast_add_f_color_rule(&dat1, 255, 0, 0, &dat2, 155, 0, 20, &colors);
            maps = NULL;
            if (params->pcurv != NULL) {
                maps = G_find_file("cell", params->pcurv, "");
                if (maps == NULL) {
                    G_warning(_("Raster map <%s> not found"), params->pcurv);
                    return -1;
                }
                Rast_write_colors(params->pcurv, maps, &colors);
 
                fprintf(stderr, "color map written\n");
 
                Rast_quantize_fp_map_range(params->pcurv, mapset, dat1, dat2,
                                           (CELL)(dat1 * MULT),
                                           (CELL)(dat2 * MULT));
                do_history(params->pcurv, input, params);
            }
 
            if (params->tcurv != NULL) {
                maps = NULL;
                maps = G_find_file("cell", params->tcurv, "");
                if (maps == NULL) {
                    G_warning(_("Raster map <%s> not found"), params->tcurv);
                    return -1;
                }
                Rast_write_colors(params->tcurv, maps, &colors);
                Rast_quantize_fp_map_range(params->tcurv, mapset, dat1, dat2,
                                           (CELL)(dat1 * MULT),
                                           (CELL)(dat2 * MULT));
 
                do_history(params->tcurv, input, params);
            }
 
            if (params->mcurv != NULL) {
                maps = NULL;
                maps = G_find_file("cell", params->mcurv, "");
                if (maps == NULL) {
                    G_warning(_("Raster map <%s> not found"), params->mcurv);
                    return -1;
                }
                Rast_write_colors(params->mcurv, maps, &colors);
                Rast_quantize_fp_map_range(params->mcurv, mapset, dat1, dat2,
                                           (CELL)(dat1 * MULT),
                                           (CELL)(dat2 * MULT));
 
                do_history(params->mcurv, input, params);
            }
        }
    }
 
    if (params->elev != NULL) {
        if (!G_find_file2("cell", params->elev, "")) {
            G_warning(_("Raster map <%s> not found"), params->elev);
            return -1;
        }
 
        Rast_short_history(params->elev, "raster", &hist);
 
        if (smooth != NULL)
            Rast_append_format_history(&hist, "tension=%f, smoothing=%s",
                                       params->fi * 1000. / (*dnorm), smooth);
        else
            Rast_append_format_history(&hist, "tension=%f",
                                       params->fi * 1000. / (*dnorm));
 
        Rast_append_format_history(&hist, "dnorm=%f, zmult=%f", *dnorm,
                                   params->zmult);
        Rast_append_format_history(&hist, "KMAX=%d, KMIN=%d, errtotal=%f",
                                   params->kmax, params->kmin,
                                   sqrt(ertot / n_points));
        Rast_append_format_history(&hist, "zmin_data=%f, zmax_data=%f", zmin,
                                   zmax);
        Rast_append_format_history(&hist, "zmin_int=%f, zmax_int=%f", zminac,
                                   zmaxac);
 
        Rast_format_history(&hist, HIST_DATSRC_1, "raster map %s", input);
 
        Rast_write_history(params->elev, &hist);
 
        Rast_free_history(&hist);
    }
 
    /* change region to initial region */
    G_verbose_message(_("Changing the region back to initial..."));
    Rast_set_output_window(winhd);
 
    return 1;
}