← Back to branch summary

~ubuntu-branches/ubuntu/wily/grass/wily

« back to all changes in this revision

Viewing changes to raster/r.surf.idw2/main.c

Tags: 7.0.0~rc1+ds1-1~exp1
http://bugs.debian.org/776349
* New upstream release candidate.
* Repack upstream tarball, remove precompiled Python objects.
* Add upstream metadata.
* Update gbp.conf and Vcs-Git URL to use the experimental branch.
* Update watch file for GRASS 7.0.
* Drop build dependencies for Tcl/Tk, add build dependencies:
  python-numpy, libnetcdf-dev, netcdf-bin, libblas-dev, liblapack-dev
* Update Vcs-Browser URL to use cgit instead of gitweb.
* Update paths to use grass70.
* Add configure options: --with-netcdf, --with-blas, --with-lapack,
  remove --with-tcltk-includes.
* Update patches for GRASS 7.
* Update copyright file, changes:
  - Update copyright years
  - Group files by license
  - Remove unused license sections
* Add patches for various typos.
* Fix desktop file with patch instead of d/rules.
* Use minimal dh rules.
* Bump Standards-Version to 3.9.6, no changes.
* Use dpkg-maintscript-helper to replace directories with symlinks.
  (closes: #776349)
* Update my email to use @debian.org address.

Show diffs side-by-side

added added

removed removed

Lines of Context:
raster/r.surf.idw2/main.c
1
 
 
2
 
/****************************************************************************
3
 
 *
4
 
 * MODULE:       r.surf.idw2
5
 
 * AUTHOR(S):    Michael Shapiro, CERL (original contributor)
6
 
 *               Roberto Flor <flor itc.it>, Markus Neteler <neteler itc.it>,
7
 
 *               Glynn Clements <glynn gclements.plus.com>, Jachym Cepicky <jachym les-ejk.cz>,
8
 
 *               Jan-Oliver Wagner <jan intevation.de>, Radim Blazek <radim.blazek gmail.com>
9
 
 * PURPOSE:      
10
 
 * COPYRIGHT:    (C) 1999-2006 by the GRASS Development Team
11
 
 *
12
 
 *               This program is free software under the GNU General Public
13
 
 *               License (>=v2). Read the file COPYING that comes with GRASS
14
 
 *               for details.
15
 
 *
16
 
 *****************************************************************************/
17
 
#include <stdlib.h>
18
 
#include <unistd.h>
19
 
#include <grass/gis.h>
20
 
#include "local_proto.h"
21
 
#include <grass/glocale.h>
22
 
 
23
 
int search_points = 12;
24
 
 
25
 
int npoints = 0;
26
 
int npoints_alloc = 0;
27
 
int nsearch;
28
 
 
29
 
struct Point
30
 
{
31
 
    double north, east;
32
 
    double z;
33
 
    double dist;
34
 
};
35
 
struct Point *points = NULL;
36
 
struct Point *list;
37
 
 
38
 
 
39
 
int main(int argc, char *argv[])
40
 
{
41
 
    int fd, maskfd;
42
 
    CELL *cell, *mask;
43
 
    struct Cell_head window;
44
 
    int row, col;
45
 
    double north, east;
46
 
    double dx, dy;
47
 
    double maxdist, dist;
48
 
    double sum1, sum2;
49
 
    int i, n, max;
50
 
    struct GModule *module;
51
 
    struct History history;
52
 
    struct
53
 
    {
54
 
        struct Option *input, *npoints, *output;
55
 
    } parm;
56
 
 
57
 
    G_gisinit(argv[0]);
58
 
 
59
 
    module = G_define_module();
60
 
    module->keywords = _("raster");
61
 
    module->description = _("Surface generation program.");
62
 
 
63
 
    parm.input = G_define_standard_option(G_OPT_R_INPUT);
64
 
 
65
 
    parm.output = G_define_standard_option(G_OPT_R_OUTPUT);
66
 
 
67
 
    parm.npoints = G_define_option();
68
 
    parm.npoints->key = "npoints";
69
 
    parm.npoints->key_desc = "count";
70
 
    parm.npoints->type = TYPE_INTEGER;
71
 
    parm.npoints->required = NO;
72
 
    parm.npoints->description = _("Number of interpolation points");
73
 
    parm.npoints->answer = "12";
74
 
 
75
 
    if (G_parser(argc, argv))
76
 
        exit(EXIT_FAILURE);
77
 
 
78
 
    /* Make sure that the current projection is not lat/long */
79
 
    if ((G_projection() == PROJECTION_LL))
80
 
        G_fatal_error(_("Lat/long databases not supported by r.surf.idw2. Use r.surf.idw instead!"));
81
 
 
82
 
    if (G_legal_filename(parm.output->answer) < 0)
83
 
        G_fatal_error(_("<%s> is an illegal file name"), parm.output->answer);
84
 
 
85
 
    if (sscanf(parm.npoints->answer, "%d", &search_points) != 1 ||
86
 
        search_points < 1)
87
 
        G_fatal_error(_("%s=%s - illegal number of interpolation points"),
88
 
                      parm.npoints->key, parm.npoints->answer);
89
 
 
90
 
    list = (struct Point *)G_calloc(search_points, sizeof(struct Point));
91
 
 
92
 
    /* read the elevation points from the input raster map */
93
 
    read_cell(parm.input->answer);
94
 
 
95
 
    if (npoints == 0)
96
 
        G_fatal_error(_("%s: no data points found"), G_program_name());
97
 
    nsearch = npoints < search_points ? npoints : search_points;
98
 
 
99
 
    /* get the window, allocate buffers, etc. */
100
 
    G_get_set_window(&window);
101
 
 
102
 
    cell = G_allocate_cell_buf();
103
 
 
104
 
    if ((maskfd = G_maskfd()) >= 0)
105
 
        mask = G_allocate_cell_buf();
106
 
    else
107
 
        mask = NULL;
108
 
 
109
 
    fd = G_open_cell_new(parm.output->answer);
110
 
    if (fd < 0)
111
 
        G_fatal_error(_("Unable to create raster map <%s>"),
112
 
                      parm.output->answer);
113
 
 
114
 
    G_message(_("Interpolating raster map <%s>... %d rows... "),
115
 
              parm.output->answer, window.rows);
116
 
 
117
 
    north = window.north - window.ns_res / 2.0;
118
 
    for (row = 0; row < window.rows; row++) {
119
 
        G_percent(row, window.rows, 2);
120
 
 
121
 
        if (mask) {
122
 
            if (G_get_map_row(maskfd, mask, row) < 0)
123
 
                G_fatal_error(_("Cannot get row"));
124
 
        }
125
 
        north += window.ns_res;
126
 
        east = window.west - window.ew_res / 2.0;
127
 
        for (col = 0; col < window.cols; col++) {
128
 
            east += window.ew_res;
129
 
            /* don't interpolate outside of the mask */
130
 
            if (mask && mask[col] == 0) {
131
 
                cell[col] = 0;
132
 
                continue;
133
 
            }
134
 
            /* fill list with first nsearch points */
135
 
            for (i = 0; i < nsearch; i++) {
136
 
                dy = points[i].north - north;
137
 
                dx = points[i].east - east;
138
 
                list[i].dist = dy * dy + dx * dx;
139
 
                list[i].z = points[i].z;
140
 
            }
141
 
            /* find the maximum distance */
142
 
            maxdist = list[max = 0].dist;
143
 
            for (n = 1; n < nsearch; n++) {
144
 
                if (maxdist < list[n].dist)
145
 
                    maxdist = list[max = n].dist;
146
 
            }
147
 
            /* go thru rest of the points now */
148
 
            for (; i < npoints; i++) {
149
 
                dy = points[i].north - north;
150
 
                dx = points[i].east - east;
151
 
                dist = dy * dy + dx * dx;
152
 
 
153
 
                if (dist < maxdist) {
154
 
                    /* replace the largest dist */
155
 
                    list[max].z = points[i].z;
156
 
                    list[max].dist = dist;
157
 
                    maxdist = list[max = 0].dist;
158
 
                    for (n = 1; n < nsearch; n++) {
159
 
                        if (maxdist < list[n].dist)
160
 
                            maxdist = list[max = n].dist;
161
 
                    }
162
 
                }
163
 
            }
164
 
 
165
 
            /* interpolate */
166
 
            sum1 = 0.0;
167
 
            sum2 = 0.0;
168
 
            for (n = 0; n < nsearch; n++) {
169
 
                if ((dist = list[n].dist)) {
170
 
                    sum1 += list[n].z / dist;
171
 
                    sum2 += 1.0 / dist;
172
 
                }
173
 
                else {
174
 
                    sum1 = list[n].z;
175
 
                    sum2 = 1.0;
176
 
                    break;
177
 
                }
178
 
            }
179
 
            cell[col] = (CELL) (sum1 / sum2 + 0.5);
180
 
        }
181
 
 
182
 
        G_put_raster_row(fd, cell, CELL_TYPE);
183
 
    }
184
 
 
185
 
    G_free(points);
186
 
    G_free(cell);
187
 
    G_close_cell(fd);
188
 
 
189
 
    /* writing history file */
190
 
    G_short_history(parm.output->answer, "raster", &history);
191
 
    G_command_history(&history);
192
 
    G_write_history(parm.output->answer, &history);
193
 
    G_done_msg(" ");
194
 
 
195
 
    exit(EXIT_SUCCESS);
196
 
}
197
 
 
198
 
int newpoint(double z, double east, double north)
199
 
{
200
 
    if (npoints_alloc <= npoints) {
201
 
        npoints_alloc += 1024;
202
 
        points = (struct Point *)G_realloc(points,
203
 
                                           npoints_alloc *
204
 
                                           sizeof(struct Point));
205
 
    }
206
 
    points[npoints].north = north;
207
 
    points[npoints].east = east;
208
 
    points[npoints].z = z;
209
 
    npoints++;
210
 
 
211
 
    return 0;
212
 
}