Loading...
Searching...
No Matches
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 G_UNUSED, double gmax G_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;
}
NULL
#define NULL
Definition ccmath.h:32
G_find_file2
const char * G_find_file2(const char *element, const char *name, const char *mapset)
Searches for a file from the mapset search list or in a specified mapset. (look but don't touch).
Definition find_file.c:234
G_find_file
const char * G_find_file(const char *element, char *name, const char *mapset)
Searches for a file from the mapset search list or in a specified mapset.
Definition find_file.c:186
G_verbose_message
void G_verbose_message(const char *msg,...)
Print a message to stderr but only if module is in verbose mode.
Definition gis/error.c:108
G_fatal_error
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition gis/error.c:159
G_warning
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition gis/error.c:203
amin1
double amin1(double, double)
Definition minmax.c:65
amax1
double amax1(double, double)
Definition minmax.c:52
G_mapset
const char * G_mapset(void)
Get current mapset name.
Definition mapset.c:33
name
const char * name
Definition named_colr.c:6
MULT
#define MULT
Definition output2d.c:30
IL_resample_output_2d
int IL_resample_output_2d(struct interp_params *params, double zmin, double zmax, double zminac, double zmaxac, double c1min, double c1max, double c2min, double c2max, double gmin, double gmax, double ertot, char *input, double *dnorm, struct Cell_head *outhd, struct Cell_head *winhd, char *smooth, int n_points)
Definition resout2d.c:45
interp_params
Definition interpf.h:70
interp_params::zmult
double zmult
Definition interpf.h:72
interp_params::Tmp_fd_xx
FILE * Tmp_fd_xx
Definition interpf.h:123
interp_params::Tmp_fd_xy
FILE * Tmp_fd_xy
Definition interpf.h:123
interp_params::pcurv
char * pcurv
Definition interpf.h:107
interp_params::deriv
int deriv
Definition interpf.h:114
interp_params::Tmp_fd_yy
FILE * Tmp_fd_yy
Definition interpf.h:123
interp_params::fi
double fi
Definition interpf.h:97
interp_params::Tmp_fd_dx
FILE * Tmp_fd_dx
Definition interpf.h:123
interp_params::nsizr
int nsizr
Definition interpf.h:92
interp_params::Tmp_fd_z
FILE * Tmp_fd_z
Definition interpf.h:123
interp_params::nsizc
int nsizc
Definition interpf.h:92
interp_params::tcurv
char * tcurv
Definition interpf.h:107
interp_params::mcurv
char * mcurv
Definition interpf.h:107
interp_params::Tmp_fd_dy
FILE * Tmp_fd_dy
Definition interpf.h:123
interp_params::aspect
char * aspect
Definition interpf.h:107
interp_params::elev
char * elev
Definition interpf.h:107
interp_params::kmin
int kmin
Definition interpf.h:84
interp_params::slope
char * slope
Definition interpf.h:107
interp_params::kmax
int kmax
Definition interpf.h:87