← Back to branch summary

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

« back to all changes in this revision

Viewing changes to raster/r.proj.seg/bordwalk.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.proj.seg/bordwalk.c
1
 
/*
2
 
 * Function added to r.proj
3
 
 * GNU GPL by Morten Hulden <morten@ngb.se>, August 2000
4
 
 *
5
 
 * bordwalk.c - projects the border cell centers of a map or region
6
 
 * to check whether they are inside the borders of another map/region 
7
 
 * in a different location, and adjusts the cell header so
8
 
 * only overlapping areas are included. The function is called by main,
9
 
 * first to project the output region into the input map and trim the
10
 
 * borders of this in order to get a smaller map, and faster and less hungry 
11
 
 * memory allocation. Then main calls the function again, but reversed,
12
 
 * to project the input map on the output region, trimming this down to 
13
 
 * the smallest possible rectangular region.
14
 
 * 
15
 
 * Simply using corner and midpoints (original r.proj) will only work 
16
 
 * between cylindrical projections. In other projections, though he input 
17
 
 * map is always a rectangular area, the projected output can be of almost 
18
 
 * any shape and its position can be rotated any way. It can even be a 
19
 
 * discontinous area.
20
 
 *
21
 
 * In many projections, especially when large areas are displayed, the edges 
22
 
 * of rectangular GRASS regions do not necessarily represent east, west, north 
23
 
 * and south. Naming the region edges accordingly (as is regions and cellhd) can be 
24
 
 * misleading. (Well, invite code readers/writers to make assumptions anyway. Don't 
25
 
 * assume north is really direction north in this code ;)
26
 
 */
27
 
 
28
 
#include <math.h>
29
 
#include <stdio.h>
30
 
#include <grass/gis.h>
31
 
#include <grass/gprojects.h>
32
 
#include <grass/glocale.h>
33
 
 
34
 
void bordwalk(struct Cell_head *from_hd, struct Cell_head *to_hd,
35
 
              struct pj_info *from_pj, struct pj_info *to_pj)
36
 
{
37
 
    double idx;
38
 
    double hx, hy;
39
 
    double xmin, xmax;
40
 
    double ymin, ymax;
41
 
 
42
 
    /* Set some (un)reasonable defaults before we walk the borders */
43
 
 
44
 
    xmax = to_hd->west - 0.000001;
45
 
    xmin = to_hd->east + 0.000001;
46
 
    ymin = to_hd->north + 0.000001;
47
 
    ymax = to_hd->south - 0.000001;
48
 
 
49
 
    /* Start walking */
50
 
 
51
 
    /* Top */
52
 
    for (idx = from_hd->west + from_hd->ew_res / 2; idx < from_hd->east;
53
 
         idx += from_hd->ew_res) {
54
 
        hx = idx;
55
 
        hy = from_hd->north - from_hd->ns_res / 2;
56
 
        if (pj_do_proj(&hx, &hy, from_pj, to_pj) < 0)
57
 
            continue;
58
 
        /* check if we are within the region, but allow for some 'almost inside' points */
59
 
        /* (should probably be a factor based on input and output resolutions) */
60
 
        if (!(hx < to_hd->west - to_hd->ew_res) &&
61
 
            !(hx > to_hd->east + to_hd->ew_res) &&
62
 
            !(hy < to_hd->south - to_hd->ns_res) &&
63
 
            !(hy > to_hd->north + to_hd->ns_res)) {
64
 
            xmin = !(hx > xmin) ? hx : xmin;
65
 
            xmax = !(hx < xmax) ? hx : xmax;
66
 
            ymin = !(hy > ymin) ? hy : ymin;
67
 
            ymax = !(hy < ymax) ? hy : ymax;
68
 
        }
69
 
    }
70
 
 
71
 
    G_debug(3, "Top: xmin: %f; xmax: %f; ymin: %f; ymax: %f", xmin, xmax,
72
 
            ymin, ymax);
73
 
 
74
 
    /* Right */
75
 
    for (idx = from_hd->north - from_hd->ns_res / 2; idx > from_hd->south;
76
 
         idx -= from_hd->ns_res) {
77
 
        hx = from_hd->east - from_hd->ew_res / 2;
78
 
        hy = idx;
79
 
        if (pj_do_proj(&hx, &hy, from_pj, to_pj) < 0)
80
 
            continue;
81
 
        if (!(hx < to_hd->west - to_hd->ew_res) &&
82
 
            !(hx > to_hd->east + to_hd->ew_res) &&
83
 
            !(hy < to_hd->south - to_hd->ns_res) &&
84
 
            !(hy > to_hd->north + to_hd->ns_res)) {
85
 
            xmin = !(hx > xmin) ? hx : xmin;
86
 
            xmax = !(hx < xmax) ? hx : xmax;
87
 
            ymin = !(hy > ymin) ? hy : ymin;
88
 
            ymax = !(hy < ymax) ? hy : ymax;
89
 
        }
90
 
    }
91
 
 
92
 
    G_debug(3, "Right: xmin: %f; xmax: %f; ymin: %f; ymax: %f", xmin, xmax,
93
 
            ymin, ymax);
94
 
 
95
 
    /* Bottom */
96
 
    for (idx = from_hd->east - from_hd->ew_res / 2; idx > from_hd->west;
97
 
         idx -= from_hd->ew_res) {
98
 
        hx = idx;
99
 
        hy = from_hd->south + from_hd->ns_res / 2;
100
 
        if (pj_do_proj(&hx, &hy, from_pj, to_pj) < 0)
101
 
            continue;
102
 
        if (!(hx < to_hd->west - to_hd->ew_res) &&
103
 
            !(hx > to_hd->east + to_hd->ew_res) &&
104
 
            !(hy < to_hd->south - to_hd->ns_res) &&
105
 
            !(hy > to_hd->north + to_hd->ns_res)) {
106
 
            xmin = !(hx > xmin) ? hx : xmin;
107
 
            xmax = !(hx < xmax) ? hx : xmax;
108
 
            ymin = !(hy > ymin) ? hy : ymin;
109
 
            ymax = !(hy < ymax) ? hy : ymax;
110
 
        }
111
 
    }
112
 
 
113
 
    G_debug(3, "Bottom: xmin: %f; xmax: %f; ymin: %f; ymax: %f", xmin, xmax,
114
 
            ymin, ymax);
115
 
 
116
 
    /* Left */
117
 
    for (idx = from_hd->south + from_hd->ns_res / 2; idx < from_hd->north;
118
 
         idx += from_hd->ns_res) {
119
 
        hx = from_hd->west + from_hd->ew_res / 2;
120
 
        hy = idx;
121
 
        if (pj_do_proj(&hx, &hy, from_pj, to_pj) < 0)
122
 
            continue;
123
 
        if (!(hx < to_hd->west - to_hd->ew_res) &&
124
 
            !(hx > to_hd->east + to_hd->ew_res) &&
125
 
            !(hy < to_hd->south - to_hd->ns_res) &&
126
 
            !(hy > to_hd->north + to_hd->ns_res)) {
127
 
            xmin = !(hx > xmin) ? hx : xmin;
128
 
            xmax = !(hx < xmax) ? hx : xmax;
129
 
            ymin = !(hy > ymin) ? hy : ymin;
130
 
            ymax = !(hy < ymax) ? hy : ymax;
131
 
        }
132
 
    }
133
 
 
134
 
    G_debug(3, "Left: xmin: %f; xmax: %f; ymin: %f; ymax: %f", xmin, xmax,
135
 
            ymin, ymax);
136
 
 
137
 
    /* check some special cases by reversing the projection */
138
 
 
139
 
    if (xmin > to_hd->west) {
140
 
        hx = to_hd->west + to_hd->ew_res / 2;
141
 
        hy = to_hd->south + (to_hd->north - to_hd->south) / 2;
142
 
        if (!(pj_do_proj(&hx, &hy, to_pj, from_pj) < 0) &&
143
 
            !(hx < from_hd->west) && !(hx > from_hd->east) &&
144
 
            !(hy < from_hd->south) && !(hy > from_hd->north))
145
 
            xmin = to_hd->west + to_hd->ew_res / 2;
146
 
    }
147
 
 
148
 
    if (xmax < to_hd->east) {
149
 
        hx = to_hd->east - to_hd->ew_res / 2;
150
 
        hy = to_hd->south + (to_hd->north - to_hd->south) / 2;
151
 
        if (!(pj_do_proj(&hx, &hy, to_pj, from_pj) < 0) &&
152
 
            !(hx < from_hd->west) && !(hx > from_hd->east) &&
153
 
            !(hy < from_hd->south) && !(hy > from_hd->north))
154
 
            xmax = to_hd->east - to_hd->ew_res / 2;
155
 
    }
156
 
 
157
 
    if (ymin > to_hd->south) {
158
 
        hx = to_hd->west + (to_hd->east - to_hd->west) / 2;
159
 
        hy = to_hd->south + to_hd->ns_res / 2;
160
 
        if (!(pj_do_proj(&hx, &hy, to_pj, from_pj) < 0) &&
161
 
            !(hx < from_hd->west) && !(hx > from_hd->east) &&
162
 
            !(hy < from_hd->south) && !(hy > from_hd->north))
163
 
            ymin = to_hd->south + to_hd->ns_res / 2;
164
 
    }
165
 
 
166
 
    if (ymax < to_hd->north) {
167
 
        hx = to_hd->west + (to_hd->east - to_hd->west) / 2;
168
 
        hy = to_hd->north - to_hd->ns_res / 2;
169
 
        if (!(pj_do_proj(&hx, &hy, to_pj, from_pj) < 0) &&
170
 
            !(hx < from_hd->west) && !(hx > from_hd->east) &&
171
 
            !(hy < from_hd->south) && !(hy > from_hd->north))
172
 
            ymax = to_hd->north - to_hd->ns_res / 2;
173
 
    }
174
 
 
175
 
    G_debug(3, "Extra check: xmin: %f; xmax: %f; ymin: %f; ymax: %f", xmin,
176
 
            xmax, ymin, ymax);
177
 
 
178
 
    /* if we still have some unresonable default minmax left, then abort */
179
 
 
180
 
    if ((xmin > to_hd->east) || (xmax < to_hd->west)
181
 
        || (ymin > to_hd->north) || (ymax < to_hd->south))
182
 
        G_fatal_error(_("Input raster map is outside current region"));
183
 
 
184
 
    if (xmin < to_hd->west + to_hd->ew_res / 2)
185
 
        xmin = to_hd->west + to_hd->ew_res / 2;
186
 
    if (xmax > to_hd->east - to_hd->ew_res / 2)
187
 
        xmax = to_hd->east - to_hd->ew_res / 2;
188
 
    if (ymin < to_hd->south + to_hd->ns_res / 2)
189
 
        ymin = to_hd->south + to_hd->ns_res / 2;
190
 
    if (ymax > to_hd->north - to_hd->ns_res / 2)
191
 
        ymax = to_hd->north - to_hd->ns_res / 2;
192
 
 
193
 
    /* adjust to edges */
194
 
 
195
 
    idx = (int)floor(G_easting_to_col(xmin, to_hd));
196
 
    xmin = G_col_to_easting(idx + 0.0, to_hd);
197
 
    idx = (int)floor(G_easting_to_col(xmax, to_hd));
198
 
    xmax = G_col_to_easting(idx + 1.0, to_hd);
199
 
    idx = (int)floor(G_northing_to_row(ymin, to_hd));
200
 
    ymin = G_row_to_northing(idx + 1.0, to_hd);
201
 
    idx = (int)floor(G_northing_to_row(ymax, to_hd));
202
 
    ymax = G_row_to_northing(idx + 0.0, to_hd);
203
 
 
204
 
    to_hd->west = (xmin < to_hd->west) ? to_hd->west : xmin;
205
 
    to_hd->east = (xmax > to_hd->east) ? to_hd->east : xmax;
206
 
    to_hd->south = (ymin < to_hd->south) ? to_hd->south : ymin;
207
 
    to_hd->north = (ymax > to_hd->north) ? to_hd->north : ymax;
208
 
 
209
 
    G_debug(3, "Final check: xmin: %f; xmax: %f; ymin: %f; ymax: %f", xmin,
210
 
            xmax, ymin, ymax);
211
 
}