GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-809ca1ee84
raster/window_map.c
Go to the documentation of this file.
1 /*!
2  * \file lib/raster/window_map.c
3  *
4  * \brief Raster Library - Window mapping functions.
5  *
6  * (C) 2001-2009 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 
14 #include <stdlib.h>
15 #include <grass/gis.h>
16 #include <grass/raster.h>
17 
18 #include "R.h"
19 
20 #define alloc_index(n) (COLUMN_MAPPING *)G_malloc((n) * sizeof(COLUMN_MAPPING))
21 
22 /*!
23  * \brief Create window mapping.
24  *
25  * Creates mapping from cell header into window. The boundaries and
26  * resolution of the two spaces do not have to be the same or aligned in
27  * any way.
28  *
29  * \param fd file descriptor
30  */
32 {
33  struct fileinfo *fcb = &R__.fileinfo[fd];
34  COLUMN_MAPPING *col;
35  int i;
36  int x;
37  double C1, C2;
38  double west, east;
39 
40  if (fcb->open_mode >= 0 && fcb->open_mode != OPEN_OLD) /* open for write? */
41  return;
42  if (fcb->open_mode == OPEN_OLD) /* already open ? */
43  G_free(fcb->col_map);
44 
45  col = fcb->col_map = alloc_index(R__.rd_window.cols);
46 
47  /*
48  * for each column in the window, go to center of the cell,
49  * compute nearest column in the data file
50  * if column is not in data file, set column to 0
51  *
52  * for lat/lon move window so that west is bigger than
53  * cellhd west.
54  */
55  west = R__.rd_window.west;
56  east = R__.rd_window.east;
58  while (west > fcb->cellhd.west + 360.0) {
59  west -= 360.0;
60  east -= 360.0;
61  }
62  while (west < fcb->cellhd.west) {
63  west += 360.0;
64  east += 360.0;
65  }
66  }
67 
69  C2 = (west - fcb->cellhd.west + R__.rd_window.ew_res / 2.0) /
70  fcb->cellhd.ew_res;
71  for (i = 0; i < R__.rd_window.cols; i++) {
72  x = C2;
73  if (C2 < x) /* adjust for rounding of negatives */
74  x--;
75  if (x < 0 || x >= fcb->cellhd.cols) /* not in data file */
76  x = -1;
77  *col++ = x + 1;
78  C2 += C1;
79  }
80 
81  /* do wrap around for lat/lon */
83 
84  while (east - 360.0 > fcb->cellhd.west) {
85  east -= 360.0;
86  west -= 360.0;
87 
88  col = fcb->col_map;
89  C2 = (west - fcb->cellhd.west + R__.rd_window.ew_res / 2.0) /
90  fcb->cellhd.ew_res;
91  for (i = 0; i < R__.rd_window.cols; i++) {
92  x = C2;
93  if (C2 < x) /* adjust for rounding of negatives */
94  x--;
95  if (x < 0 || x >= fcb->cellhd.cols) /* not in data file */
96  x = -1;
97  if (*col == 0) /* only change those not already set */
98  *col = x + 1;
99  col++;
100  C2 += C1;
101  }
102  }
103  }
104 
105  G_debug(3, "create window mapping (%d columns)", R__.rd_window.cols);
106  /* for (i = 0; i < R__.rd_window.cols; i++)
107  fprintf(stderr, "%s%ld", i % 15 ? " " : "\n", (long)fcb->col_map[i]);
108  fprintf(stderr, "\n");
109  */
110 
111  /* compute C1,C2 for row window mapping */
112  fcb->C1 = R__.rd_window.ns_res / fcb->cellhd.ns_res;
113  fcb->C2 =
114  (fcb->cellhd.north - R__.rd_window.north + R__.rd_window.ns_res / 2.0) /
115  fcb->cellhd.ns_res;
116 }
117 
118 /*!
119  * \brief Loops rows until mismatch?.
120  *
121  * This routine works fine if the mask is not set. It may give
122  * incorrect results with a mask, since the mask row may have a
123  * different repeat value. The issue can be fixed by doing it for the
124  * mask as well and using the smaller value.
125  *
126  * \param fd file descriptor
127  * \param row starting row
128  *
129  * \return number of rows completed
130  */
131 int Rast_row_repeat_nomask(int fd, int row)
132 {
133  struct fileinfo *fcb = &R__.fileinfo[fd];
134  double f;
135  int r1, r2;
136  int count;
137 
138  count = 1;
139 
140  /* r1 is the row in the raster map itself.
141  * r2 is the next row(s) in the raster map
142  * see get_row.c for details on this calculation
143  */
144  f = row * fcb->C1 + fcb->C2;
145  r1 = f;
146  if (f < r1)
147  r1--;
148 
149  while (++row < R__.rd_window.rows) {
150  f = row * fcb->C1 + fcb->C2;
151  r2 = f;
152  if (f < r2)
153  r2--;
154  if (r1 != r2)
155  break;
156 
157  count++;
158  }
159 
160  return count;
161 }
int COLUMN_MAPPING
Definition: R.h:17
#define OPEN_OLD
Definition: R.h:106
void G_free(void *)
Free allocated memory.
Definition: gis/alloc.c:150
int G_debug(int, const char *,...) __attribute__((format(printf
#define PROJECTION_LL
Projection code - Latitude-Longitude.
Definition: gis.h:130
int count
void Rast__create_window_mapping(int fd)
Create window mapping.
int Rast_row_repeat_nomask(int fd, int row)
Loops rows until mismatch?.
#define alloc_index(n)
double ew_res
Resolution - east to west cell size for 2D data.
Definition: gis.h:473
double north
Extent coordinates (north)
Definition: gis.h:483
double east
Extent coordinates (east)
Definition: gis.h:487
double ns_res
Resolution - north to south cell size for 2D data.
Definition: gis.h:477
int rows
Number of rows for 2D data.
Definition: gis.h:452
int cols
Number of columns for 2D data.
Definition: gis.h:456
int proj
Projection code.
Definition: gis.h:469
double west
Extent coordinates (west)
Definition: gis.h:489
Definition: R.h:87
struct fileinfo * fileinfo
Definition: R.h:101
struct Cell_head rd_window
Definition: R.h:97
Definition: R.h:53
struct Cell_head cellhd
Definition: R.h:55
COLUMN_MAPPING * col_map
Definition: R.h:63
int open_mode
Definition: R.h:54
double C2
Definition: R.h:64
double C1
Definition: R.h:64
#define x