1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
|
#include <unistd.h>
#include <string.h>
#include <grass/raster.h>
#include <grass/glocale.h>
#include "global.h"
/* Modified to support Grass 5.0 fp format 11 april 2000
*
* Pierre de Mouveaux - pmx@audiovu.com
*
*/
int rectify(char *name, char *mapset, char *result, int order, char *interp_method)
{
struct Cell_head cellhd;
int ncols, nrows;
int row, col;
double row_idx, col_idx;
int infd, cell_size, outfd;
void *trast, *tptr;
double n1, e1, nx, ex;
struct cache *ibuffer;
select_current_env();
Rast_get_cellhd(name, mapset, &cellhd);
/* open the file to be rectified
* set window to cellhd first to be able to read file exactly
*/
Rast_set_input_window(&cellhd);
infd = Rast_open_old(name, mapset);
map_type = Rast_get_map_type(infd);
cell_size = Rast_cell_size(map_type);
ibuffer = readcell(infd, seg_mb);
Rast_close(infd); /* (pmx) 17 april 2000 */
G_message(_("Rectify <%s@%s> (location <%s>)"),
name, mapset, G_location());
select_target_env();
G_message(_("into <%s@%s> (location <%s>) ..."),
result, G_mapset(), G_location());
nrows = target_window.rows;
ncols = target_window.cols;
if (strcmp(interp_method, "nearest") != 0) {
map_type = DCELL_TYPE;
cell_size = Rast_cell_size(map_type);
}
/* open the result file into target window
* this open must be first since we change the window later
* raster maps open for writing are not affected by window changes
* but those open for reading are
*/
outfd = Rast_open_new(result, map_type);
trast = Rast_allocate_output_buf(map_type);
for (row = 0; row < nrows; row++) {
n1 = target_window.north - (row + 0.5) * target_window.ns_res;
G_percent(row, nrows, 2);
Rast_set_null_value(trast, ncols, map_type);
tptr = trast;
for (col = 0; col < ncols; col++) {
e1 = target_window.west + (col + 0.5) * target_window.ew_res;
/* backwards transformation of target cell center */
if (order == 0)
I_georef_tps(e1, n1, &ex, &nx, E21_t, N21_t, &cp, 0);
else
I_georef(e1, n1, &ex, &nx, E21, N21, order);
/* convert to row/column indices of source raster */
row_idx = (cellhd.north - nx) / cellhd.ns_res;
col_idx = (ex - cellhd.west) / cellhd.ew_res;
/* resample data point */
interpolate(ibuffer, tptr, map_type, &row_idx, &col_idx, &cellhd);
tptr = G_incr_void_ptr(tptr, cell_size);
}
Rast_put_row(outfd, trast, map_type);
}
G_percent(1, 1, 1);
Rast_close(outfd); /* (pmx) 17 april 2000 */
G_free(trast);
close(ibuffer->fd);
G_free(ibuffer);
Rast_get_cellhd(result, G_mapset(), &cellhd);
if (cellhd.proj == 0) { /* x,y imagery */
cellhd.proj = target_window.proj;
cellhd.zone = target_window.zone;
}
if (target_window.proj != cellhd.proj) {
cellhd.proj = target_window.proj;
G_warning(_("Raster map <%s@%s>: projection don't match current settings"),
name, mapset);
}
if (target_window.zone != cellhd.zone) {
cellhd.zone = target_window.zone;
G_warning(_("Raster map <%s@%s>: zone don't match current settings"),
name, mapset);
}
select_current_env();
return 1;
}
|