← Back to branch summary

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

« back to all changes in this revision

Viewing changes to lib/gis/window_map.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/window_map.c
1
 
 
2
 
/**
3
 
 * \file window_map.c
4
 
 *
5
 
 * \brief GIS Library - Window mapping functions.
6
 
 *
7
 
 * (C) 2001-2008 by the GRASS Development Team
8
 
 *
9
 
 * This program is free software under the GNU General Public License
10
 
 * (>=v2). Read the file COPYING that comes with GRASS for details.
11
 
 *
12
 
 * \author GRASS GIS Development Team
13
 
 *
14
 
 * \date 1999-2008
15
 
 */
16
 
 
17
 
#include <stdlib.h>
 
1
/*!
 
2
  \file lib/gis/window_map.c
 
3
  
 
4
  \brief GIS Library - Window mapping functions.
 
5
  
 
6
  (C) 2001-2009, 2011 by the GRASS Development Team
 
7
  
 
8
  This program is free software under the GNU General Public License
 
9
  (>=v2). Read the file COPYING that comes with GRASS for details.
 
10
  
 
11
  \author Original author CERL
 
12
*/
 
13
 
18
14
#include <grass/gis.h>
 
15
 
19
16
#include "G.h"
20
17
 
21
 
 
22
 
#define alloc_index(n) (COLUMN_MAPPING *) G_malloc((n)*sizeof(COLUMN_MAPPING))
23
 
 
24
 
 
25
 
/**
26
 
 * \brief Create window mapping.
27
 
 *
28
 
 * Creates mapping from cell header into window. The boundaries and 
29
 
 * resolution of the two spaces do not have to be the same or aligned in 
30
 
 * any way.
31
 
 *
32
 
 * \param[in] fd file descriptor
33
 
 * \return always returns 0
34
 
 */
35
 
 
36
 
int G__create_window_mapping(int fd)
37
 
{
38
 
    struct fileinfo *fcb = &G__.fileinfo[fd];
39
 
    COLUMN_MAPPING *col;
40
 
    int i;
41
 
    int x;
42
 
    double C1, C2;
43
 
    double west;
44
 
 
45
 
    G__init_window();
46
 
 
47
 
    if (fcb->open_mode >= 0 && fcb->open_mode != OPEN_OLD)      /* open for write? */
48
 
        return 0;
49
 
    if (fcb->open_mode == OPEN_OLD)     /* already open ? */
50
 
        G_free(fcb->col_map);
51
 
 
52
 
    col = fcb->col_map = alloc_index(G__.window.cols);
53
 
 
54
 
    /*
55
 
     * for each column in the window, go to center of the cell,
56
 
     * compute nearest column in the data file
57
 
     * if column is not in data file, set column to 0
58
 
     *
59
 
     * for lat/lon move window so that west is bigger than
60
 
     * cellhd west.
61
 
     */
62
 
    west = G__.window.west;
63
 
    if (G__.window.proj == PROJECTION_LL) {
64
 
        while (west > fcb->cellhd.west + 360.0)
65
 
            west -= 360.0;
66
 
        while (west < fcb->cellhd.west)
67
 
            west += 360.0;
68
 
    }
69
 
 
70
 
    C1 = G__.window.ew_res / fcb->cellhd.ew_res;
71
 
    C2 = (west - fcb->cellhd.west +
72
 
          G__.window.ew_res / 2.0) / fcb->cellhd.ew_res;
73
 
    for (i = 0; i < G__.window.cols; i++) {
74
 
        x = C2;
75
 
        if (C2 < x)             /* adjust for rounding of negatives */
76
 
            x--;
77
 
        if (x < 0 || x >= fcb->cellhd.cols)     /* not in data file */
78
 
            x = -1;
79
 
        *col++ = x + 1;
80
 
        C2 += C1;
81
 
    }
82
 
 
83
 
    /* do wrap around for lat/lon */
84
 
    if (G__.window.proj == PROJECTION_LL) {
85
 
        col = fcb->col_map;
86
 
        C2 = (west - 360.0 - fcb->cellhd.west +
87
 
              G__.window.ew_res / 2.0) / fcb->cellhd.ew_res;
88
 
        for (i = 0; i < G__.window.cols; i++) {
89
 
            x = C2;
90
 
            if (C2 < x)         /* adjust for rounding of negatives */
91
 
                x--;
92
 
            if (x < 0 || x >= fcb->cellhd.cols) /* not in data file */
93
 
                x = -1;
94
 
            if (*col == 0)      /* only change those not already set */
95
 
                *col = x + 1;
96
 
            col++;
97
 
            C2 += C1;
98
 
        }
99
 
    }
100
 
 
101
 
    G_debug(3, "create window mapping (%d columns)", G__.window.cols);
102
 
/*  for (i = 0; i < G__.window.cols; i++)
103
 
        fprintf(stderr, "%s%ld", i % 15 ? " " : "\n", (long)fcb->col_map[i]);
104
 
    fprintf(stderr, "\n");
 
18
/*!
 
19
  \brief Adjust east longitude.
 
20
 
 
21
  This routine returns an equivalent <i>east</i> that is larger, but
 
22
  no more than 360 larger than the <i>west</i> coordinate.
 
23
 
 
24
  <b>Note:</b> This routine should be used only with
 
25
  latitude-longitude coordinates.
 
26
 
 
27
  \param east east coordinate
 
28
  \param west west coordinate
 
29
  
 
30
  \return east coordinate
105
31
*/
106
 
 
107
 
    /* compute C1,C2 for row window mapping */
108
 
    fcb->C1 = G__.window.ns_res / fcb->cellhd.ns_res;
109
 
    fcb->C2 =
110
 
        (fcb->cellhd.north - G__.window.north +
111
 
         G__.window.ns_res / 2.0) / fcb->cellhd.ns_res;
112
 
 
113
 
    return 0;
114
 
}
115
 
 
116
 
 
117
 
/**
118
 
 * \brief Northing to row.
119
 
 *
120
 
 * Converts a <b>north</b>ing relative to a <b>window</b> to a row.<br>
121
 
 * <b>Note:</b> The result is a double. Casting it to an integer will 
122
 
 * give the row number.
123
 
 *
124
 
 * \param[in] north
125
 
 * \param[in] window
126
 
 * \return row id
127
 
 */
128
 
 
129
 
double G_northing_to_row(double north, const struct Cell_head *window)
130
 
{
131
 
    return (window->north - north) / window->ns_res;
132
 
}
133
 
 
134
 
 
135
 
/**
136
 
 * \brief Adjust east longitude.
137
 
 *
138
 
 * This routine returns an equivalent <b>east</b> that is
139
 
 * larger, but no more than 360 larger than the <b>west</b> 
140
 
 * coordinate.<br>
141
 
 * <b>Note:</b> This routine should be used only with latitude-longitude 
142
 
 * coordinates.
143
 
 *
144
 
 * \param[in,out] east
145
 
 * \param[in] west
146
 
 * \return east coordinate
147
 
 */
148
 
 
149
32
double G_adjust_east_longitude(double east, double west)
150
33
{
151
34
    while (east > west + 360.0)
156
39
    return east;
157
40
}
158
41
 
159
 
 
160
 
/**
161
 
 * \brief Returns east larger than west.
162
 
 *
163
 
 * If the region projection is <i>PROJECTION_LL</i>, then this routine 
164
 
 * returns an equivalent <b>east</b> that is larger, but no more than 
165
 
 * 360 degrees larger, than the coordinate for the western edge of the 
166
 
 * region. Otherwise no adjustment is made and the original
167
 
 * <b>east</b> is returned.
168
 
 *
169
 
 * \param[in,out] east
170
 
 * \param[in] window
171
 
 * \return east coordinate
172
 
 */
173
 
 
 
42
/*!
 
43
  \brief Returns east larger than west.
 
44
  
 
45
  If the region projection is <tt>PROJECTION_LL</tt>, then this
 
46
  routine returns an equivalent <i>east</i> that is larger, but no
 
47
  more than 360 degrees larger, than the coordinate for the western
 
48
  edge of the region. Otherwise no adjustment is made and the original
 
49
  <i>east</i> is returned.
 
50
  
 
51
  \param east east coordinate
 
52
  \param window pointer to Cell_head
 
53
  
 
54
  \return east coordinate
 
55
*/
174
56
double G_adjust_easting(double east, const struct Cell_head *window)
175
57
{
176
58
    if (window->proj == PROJECTION_LL) {
182
64
    return east;
183
65
}
184
66
 
185
 
 
186
 
/**
187
 
 * \brief Easting to column.
188
 
 *
189
 
 * Converts <b>east</b> relative to a <b>window</b> to a column.<br>
190
 
 * <b>Note:</b> The result is a <i>double</i>. Casting it to an 
191
 
 * <i>int</i> will give the column number.
192
 
 *
193
 
 * \param[in] east
194
 
 * \param[in] window
195
 
 * \return column id
196
 
 */
197
 
 
198
 
double G_easting_to_col(double east, const struct Cell_head *window)
199
 
{
200
 
    east = G_adjust_easting(east, window);
201
 
 
202
 
    return (east - window->west) / window->ew_res;
203
 
}
204
 
 
205
 
 
206
 
/**
207
 
 * \brief Row to northing.
208
 
 *
209
 
 * Converts a <b>row</b> relative to a <b>window</b> to a
210
 
 * northing.<br>
211
 
 * <b>Note:</b> row is a double:
212
 
 *  - row+0.0 will return the northing for the northern edge of the row.<br>
213
 
 *  - row+0.5 will return the northing for the center of the row.<br>
214
 
 *  - row+1.0 will return the northing for the southern edge of the row.<br>
215
 
 * <b>Note:</b> The result is a <i>double</i>. Casting it to an 
216
 
 * <i>int</i> will give the column number.
217
 
 *
218
 
 * \param[in] row
219
 
 * \param[in] window
220
 
 * \return north coordinate
221
 
 */
222
 
 
223
 
double G_row_to_northing(double row, const struct Cell_head *window)
224
 
{
225
 
    return window->north - row * window->ns_res;
226
 
}
227
 
 
228
 
 
229
 
/**
230
 
 * \brief Column to easting.
231
 
 *
232
 
 * Converts a <b>col</b> relative to a <b>window</b> to an easting.<br>
233
 
 * <b>Note:</b> <b>col</b> is a <i>double</i>:<br>
234
 
 *  - col+0.0 will return the easting for the western edge of the column.<br>
235
 
 *  - col+0.5 will return the easting for the center of the column.<br>
236
 
 *  - col+1.0 will return the easting for the eastern edge of the column.<br>
237
 
 *
238
 
 * \param[in] col
239
 
 * \param[in] window
240
 
 * \return east coordinate
241
 
 */
242
 
 
243
 
double G_col_to_easting(double col, const struct Cell_head *window)
244
 
{
245
 
    return window->west + col * window->ew_res;
246
 
}
247
 
 
248
 
 
249
 
/**
250
 
 * \brief Number of rows in active window.
251
 
 *
252
 
 * This routine returns the number of rows in the active module window. 
253
 
 * Before raster files can be read or written, it is necessary to
254
 
 * known how many rows are in the active window. For example:<br>
255
 
 \code  
256
 
  int nrows, cols;
257
 
  int row, col; 
258
 
  nrows = G_window_rows( ); 
259
 
  ncols = G_window_cols( ); 
260
 
  for (row = 0; row < nrows; row++)
261
 
  {
262
 
  <i>read</i> row ...
263
 
  for (col = 0; col < ncols; col++)
264
 
  {
265
 
  process col ...
266
 
  }
267
 
  }
268
 
 \endcode 
269
 
 *
270
 
 *  \return number of rows
271
 
 */
272
 
 
273
 
int G_window_rows(void)
274
 
{
275
 
    G__init_window();
276
 
 
277
 
    return G__.window.rows;
278
 
}
279
 
 
280
 
 
281
 
/**
282
 
 * \brief Number of columns in active window.
283
 
 *
284
 
 * These
285
 
 * routines return the number of rows and columns (respectively) in the active
286
 
 * module region. Before raster maps can be read or written, it is necessary to
287
 
 * known how many rows and columns are in the active region. For example:<br>
288
 
 \code  
289
 
  int nrows, cols;
290
 
  int row, col; 
291
 
  nrows = G_window_rows( ); 
292
 
  ncols = G_window_cols( ); 
293
 
  for (row = 0; row < nrows; row++)
294
 
  {
295
 
  <i>read</i> row ...
296
 
  for (col = 0; col < ncols; col++)
297
 
  {
298
 
  process col ...
299
 
  }
300
 
  }
301
 
 \endcode 
302
 
 *
303
 
 * \return number of columns
304
 
 */
305
 
 
306
 
int G_window_cols(void)
307
 
{
308
 
    G__init_window();
309
 
 
310
 
    return G__.window.cols;
311
 
}
312
 
 
313
 
 
314
 
/**
315
 
 * \brief Initialize window.
316
 
 *
317
 
 * \return always returns 0
318
 
 */
319
 
 
320
 
int G__init_window(void)
321
 
{
322
 
    if (!G__.window_set) {
323
 
        G__.window_set = 1;
324
 
        G_get_window(&G__.window);
325
 
    }
326
 
 
327
 
    return 0;
328
 
}
329
 
 
330
 
 
331
 
/**
332
 
 * \brief Loops rows until mismatch?.
333
 
 *
334
 
 * This routine works fine if the mask is not set. It may give incorrect 
335
 
 * results with a mask, since the mask row may have a different repeat 
336
 
 * value. The issue can be fixed by doing it for the mask as well and 
337
 
 * using the smaller value.
338
 
 *
339
 
 * \param[in] fd file descriptor
340
 
 * \param[in] row starting row
341
 
 * \return number of rows completed
342
 
 */
343
 
 
344
 
int G_row_repeat_nomask(int fd, int row)
345
 
{
346
 
    struct fileinfo *fcb = &G__.fileinfo[fd];
347
 
    double f;
348
 
    int r1, r2;
349
 
    int count;
350
 
 
351
 
    count = 1;
352
 
 
353
 
    /* r1 is the row in the raster map itself.
354
 
     * r2 is the next row(s) in the raster map
355
 
     * see get_row.c for details on this calculation
356
 
     */
357
 
    f = row * fcb->C1 + fcb->C2;
358
 
    r1 = f;
359
 
    if (f < r1)
360
 
        r1--;
361
 
 
362
 
    while (++row < G__.window.rows) {
363
 
        f = row * fcb->C1 + fcb->C2;
364
 
        r2 = f;
365
 
        if (f < r2)
366
 
            r2--;
367
 
        if (r1 != r2)
368
 
            break;
369
 
 
370
 
        count++;
371
 
    }
372
 
 
373
 
    return count;
374
 
}
 
67
/*!
 
68
  \brief Initialize window (region).
 
69
*/
 
70
void G__init_window(void)
 
71
{
 
72
    if (G_is_initialized(&G__.window_set))
 
73
        return;
 
74
    
 
75
    G_get_window(&G__.window);
 
76
 
 
77
    G_initialize_done(&G__.window_set);
 
78
}
 
79