← Back to branch summary

~ubuntu-branches/ubuntu/vivid/grass/vivid-proposed

« back to all changes in this revision

Viewing changes to lib/gis/sample.c

  • Committer: Package Import Robot
  • Author(s): Bas Couwenberg
  • Date: 2015-02-20 23:12:08 UTC
  • mfrom: (8.2.6 experimental)
  • Revision ID: package-import@ubuntu.com-20150220231208-1u6qvqm84v430b10
Tags: 7.0.0-1~exp1
* New upstream release.
* Update python-ctypes-ternary.patch to use if/else instead of and/or.
* Drop check4dev patch, rely on upstream check.
* Add build dependency on libpq-dev to grass-dev for libpq-fe.h.
* Drop patches applied upstream, refresh remaining patches.
* Update symlinks for images switched from jpg to png.

Show diffs side-by-side

added added

removed removed

Lines of Context:
lib/gis/sample.c
1
 
/*!
2
 
   \file sample.c
3
 
 
4
 
   \brief GIS library - sampling methods (extract a cell value from
5
 
   raster map)
6
 
 
7
 
   1/2006: moved to libgis from v.sample/v.drape for clone removal
8
 
 
9
 
   (C) 2001-2008 by the GRASS Development Team
10
 
 
11
 
   This program is free software under the 
12
 
   GNU General Public License (>=v2). 
13
 
   Read the file COPYING that comes with GRASS
14
 
   for details.
15
 
 
16
 
   \author James Darrell McCauley <darrell mccauley-usa.com>, http://mccauley-usa.com/
17
 
 
18
 
   \date 1994
19
 
 */
20
 
 
21
 
#include <string.h>
22
 
#include <unistd.h>
23
 
#include <math.h>
24
 
#include <grass/gis.h>
25
 
#include <grass/glocale.h>
26
 
 
27
 
/* prototypes */
28
 
static double scancatlabel(const char *);
29
 
static void raster_row_error(const struct Cell_head *window, double north,
30
 
                             double east);
31
 
 
32
 
 
33
 
/*!
34
 
 *  \brief Extract a cell value from raster map.
35
 
 *
36
 
 *  Extract a cell value from raster map at given northing and easting
37
 
 *  with a sampled 3x3 window using a specified interpolation method.
38
 
 *
39
 
 *  \param fd file descriptor
40
 
 *  \param window region settings
41
 
 *  \param cats categories
42
 
 *  \param north northing position
43
 
 *  \param east easting position
44
 
 *  \param usedesc flag to scan category label
45
 
 *  \param itype interpolation method
46
 
 *
47
 
 *  \return cell value at given position
48
 
 */
49
 
 
50
 
DCELL G_get_raster_sample(
51
 
    int fd,
52
 
    const struct Cell_head *window,
53
 
    struct Categories *cats,
54
 
    double north, double east,
55
 
    int usedesc, INTERP_TYPE itype)
56
 
{
57
 
    double retval;
58
 
 
59
 
    switch (itype) {
60
 
    case NEAREST:
61
 
        retval = G_get_raster_sample_nearest(fd, window, cats, north, east, usedesc);
62
 
        break;
63
 
    case BILINEAR:
64
 
        retval = G_get_raster_sample_bilinear(fd, window, cats, north, east, usedesc);
65
 
        break;
66
 
    case CUBIC:
67
 
        retval = G_get_raster_sample_cubic(fd, window, cats, north, east, usedesc);
68
 
        break;
69
 
    default:
70
 
        G_fatal_error("G_get_raster_sample: %s",
71
 
                      _("Unknown interpolation type"));
72
 
    }
73
 
 
74
 
    return retval;
75
 
}
76
 
 
77
 
 
78
 
DCELL G_get_raster_sample_nearest(
79
 
    int fd,
80
 
    const struct Cell_head *window,
81
 
    struct Categories *cats,
82
 
    double north, double east, int usedesc)
83
 
{
84
 
    int row, col;
85
 
    DCELL result;
86
 
    DCELL *maprow = G_allocate_d_raster_buf();
87
 
 
88
 
    /* convert northing and easting to row and col, resp */
89
 
    row = (int)floor(G_northing_to_row(north, window));
90
 
    col = (int)floor(G_easting_to_col(east, window));
91
 
 
92
 
    if (row < 0 || row >= G_window_rows() ||
93
 
        col < 0 || col >= G_window_cols()) {
94
 
        G_set_d_null_value(&result, 1);
95
 
        goto done;
96
 
    }
97
 
 
98
 
    if (G_get_d_raster_row(fd, maprow, row) < 0)
99
 
        raster_row_error(window, north, east);
100
 
 
101
 
    if (G_is_d_null_value(&maprow[col])) {
102
 
        G_set_d_null_value(&result, 1);
103
 
        goto done;
104
 
    }
105
 
 
106
 
    if (usedesc) {
107
 
        char *buf = G_get_cat(maprow[col], cats);
108
 
 
109
 
        G_squeeze(buf);
110
 
        result = scancatlabel(buf);
111
 
    }
112
 
    else
113
 
        result = maprow[col];
114
 
 
115
 
done:
116
 
    G_free(maprow);
117
 
 
118
 
    return result;
119
 
}
120
 
 
121
 
 
122
 
DCELL G_get_raster_sample_bilinear(
123
 
    int fd,
124
 
    const struct Cell_head *window,
125
 
    struct Categories *cats,
126
 
    double north, double east, int usedesc)
127
 
{
128
 
    int row, col;
129
 
    double grid[2][2];
130
 
    DCELL *arow = G_allocate_d_raster_buf();
131
 
    DCELL *brow = G_allocate_d_raster_buf();
132
 
    double frow, fcol, trow, tcol;
133
 
    DCELL result;
134
 
 
135
 
    frow = G_northing_to_row(north, window);
136
 
    fcol = G_easting_to_col(east, window);
137
 
 
138
 
    /* convert northing and easting to row and col, resp */
139
 
    row = (int)floor(frow - 0.5);
140
 
    col = (int)floor(fcol - 0.5);
141
 
 
142
 
    trow = frow - row - 0.5;
143
 
    tcol = fcol - col - 0.5;
144
 
 
145
 
    if (row < 0 || row + 1 >= G_window_rows() ||
146
 
        col < 0 || col + 1 >= G_window_cols()) {
147
 
        G_set_d_null_value(&result, 1);
148
 
        goto done;
149
 
    }
150
 
 
151
 
    if (G_get_d_raster_row(fd, arow, row) < 0)
152
 
        raster_row_error(window, north, east);
153
 
    if (G_get_d_raster_row(fd, brow, row + 1) < 0)
154
 
        raster_row_error(window, north, east);
155
 
 
156
 
    if (G_is_d_null_value(&arow[col]) || G_is_d_null_value(&arow[col + 1]) ||
157
 
        G_is_d_null_value(&brow[col]) || G_is_d_null_value(&brow[col + 1])) {
158
 
        G_set_d_null_value(&result, 1);
159
 
        goto done;
160
 
    }
161
 
 
162
 
    /*-
163
 
     * now were ready to do bilinear interpolation over
164
 
     * arow[col], arow[col+1],
165
 
     * brow[col], brow[col+1]
166
 
     */
167
 
 
168
 
    if (usedesc) {
169
 
        char *buf;
170
 
 
171
 
        G_squeeze(buf = G_get_cat((int)arow[col], cats));
172
 
        grid[0][0] = scancatlabel(buf);
173
 
        G_squeeze(buf = G_get_cat((int)arow[col + 1], cats));
174
 
        grid[0][1] = scancatlabel(buf);
175
 
        G_squeeze(buf = G_get_cat((int)brow[col], cats));
176
 
        grid[1][0] = scancatlabel(buf);
177
 
        G_squeeze(buf = G_get_cat((int)brow[col + 1], cats));
178
 
        grid[1][1] = scancatlabel(buf);
179
 
    }
180
 
    else {
181
 
        grid[0][0] = arow[col];
182
 
        grid[0][1] = arow[col + 1];
183
 
        grid[1][0] = brow[col];
184
 
        grid[1][1] = brow[col + 1];
185
 
    }
186
 
 
187
 
    result = G_interp_bilinear(tcol, trow,
188
 
                               grid[0][0], grid[0][1], grid[1][0], grid[1][1]);
189
 
 
190
 
done:
191
 
    G_free(arow);
192
 
    G_free(brow);
193
 
 
194
 
    return result;
195
 
}
196
 
 
197
 
DCELL G_get_raster_sample_cubic(
198
 
    int fd,
199
 
    const struct Cell_head *window,
200
 
    struct Categories *cats,
201
 
    double north, double east, int usedesc)
202
 
{
203
 
    int i, j, row, col;
204
 
    double grid[4][4];
205
 
    DCELL *rows[4];
206
 
    double frow, fcol, trow, tcol;
207
 
    DCELL result;
208
 
 
209
 
    for (i = 0; i < 4; i++)
210
 
        rows[i] = G_allocate_d_raster_buf();
211
 
 
212
 
    frow = G_northing_to_row(north, window);
213
 
    fcol = G_easting_to_col(east, window);
214
 
 
215
 
    /* convert northing and easting to row and col, resp */
216
 
    row = (int)floor(frow - 1.5);
217
 
    col = (int)floor(fcol - 1.5);
218
 
 
219
 
    trow = frow - row - 1.5;
220
 
    tcol = fcol - col - 1.5;
221
 
 
222
 
    if (row < 0 || row + 3 >= G_window_rows() ||
223
 
        col < 0 || col + 3 >= G_window_cols()) {
224
 
        G_set_d_null_value(&result, 1);
225
 
        goto done;
226
 
    }
227
 
 
228
 
    for (i = 0; i < 4; i++)
229
 
        if (G_get_d_raster_row(fd, rows[i], row + i) < 0)
230
 
            raster_row_error(window, north, east);
231
 
 
232
 
    for (i = 0; i < 4; i++)
233
 
        for (j = 0; j < 4; j++)
234
 
            if (G_is_d_null_value(&rows[i][col + j])) {
235
 
                G_set_d_null_value(&result, 1);
236
 
                goto done;
237
 
            }
238
 
 
239
 
    /*
240
 
     * now were ready to do cubic interpolation over
241
 
     * arow[col], arow[col+1], arow[col+2], arow[col+3],
242
 
     * brow[col], brow[col+1], brow[col+2], brow[col+3],
243
 
     * crow[col], crow[col+1], crow[col+2], crow[col+3],
244
 
     * drow[col], drow[col+1], drow[col+2], drow[col+3],
245
 
     */
246
 
 
247
 
    if (usedesc) {
248
 
        char *buf;
249
 
 
250
 
        for (i = 0; i < 4; i++) {
251
 
            for (j = 0; j < 4; j++) {
252
 
                G_squeeze(buf = G_get_cat(rows[i][col + j], cats));
253
 
                grid[i][j] = scancatlabel(buf);
254
 
            }
255
 
        }
256
 
    }
257
 
    else {
258
 
        for (i = 0; i < 4; i++)
259
 
            for (j = 0; j < 4; j++)
260
 
                grid[i][j] = rows[i][col + j];
261
 
    }
262
 
 
263
 
    result = G_interp_bicubic(tcol, trow,
264
 
                              grid[0][0], grid[0][1], grid[0][2], grid[0][3],
265
 
                              grid[1][0], grid[1][1], grid[1][2], grid[1][3],
266
 
                              grid[2][0], grid[2][1], grid[2][2], grid[2][3],
267
 
                              grid[3][0], grid[3][1], grid[3][2], grid[3][3]);
268
 
 
269
 
done:
270
 
    for (i = 0; i < 4; i++)
271
 
        G_free(rows[i]);
272
 
 
273
 
    return result;
274
 
}
275
 
 
276
 
 
277
 
static double scancatlabel(const char *str)
278
 
{
279
 
    double val;
280
 
 
281
 
    if (strcmp(str, "no data") != 0)
282
 
        sscanf(str, "%lf", &val);
283
 
    else {
284
 
        G_warning(_("\"no data\" label found; setting to zero"));
285
 
        val = 0.0;
286
 
    }
287
 
 
288
 
    return val;
289
 
}
290
 
 
291
 
 
292
 
static void raster_row_error(const struct Cell_head *window, double north,
293
 
                             double east)
294
 
{
295
 
    G_debug(3, "DIAG: \tRegion is: n=%g s=%g e=%g w=%g",
296
 
            window->north, window->south, window->east, window->west);
297
 
    G_debug(3, "      \tData point is north=%g east=%g", north, east);
298
 
 
299
 
    G_fatal_error(_("Problem reading raster map"));
300
 
}