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