[GRASS-SVN] r72519 - grass/trunk/raster/r.proj
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Mar 23 02:25:34 PDT 2018
Author: mmetz
Date: 2018-03-23 02:25:34 -0700 (Fri, 23 Mar 2018)
New Revision: 72519
Modified:
grass/trunk/raster/r.proj/bordwalk.c
grass/trunk/raster/r.proj/main.c
grass/trunk/raster/r.proj/r.proj.h
Log:
r.proj: use new GRASS API for coordinate transformation
Modified: grass/trunk/raster/r.proj/bordwalk.c
===================================================================
--- grass/trunk/raster/r.proj/bordwalk.c 2018-03-23 09:25:04 UTC (rev 72518)
+++ grass/trunk/raster/r.proj/bordwalk.c 2018-03-23 09:25:34 UTC (rev 72519)
@@ -84,14 +84,17 @@
}
static void proj_update(const struct pj_info *from_pj, const struct pj_info *to_pj,
+ const struct pj_info *trans_pj, int dir,
struct Cell_head *to_hd, double hx, double hy)
{
- if (pj_do_proj(&hx, &hy, from_pj, to_pj) < 0)
+ if (GPJ_transform(from_pj, to_pj, trans_pj, dir,
+ &hx, &hy, NULL) < 0)
return;
update(to_hd, hx, hy);
}
void bordwalk1(const struct pj_info *from_pj, const struct pj_info *to_pj,
+ const struct pj_info *trans_pj, int dir,
const struct Cell_head *from_hd, struct Cell_head *to_hd)
{
double idx;
@@ -99,7 +102,7 @@
/* Top */
for (idx = from_hd->west + from_hd->ew_res / 2; idx < from_hd->east;
idx += from_hd->ew_res)
- proj_update(from_pj, to_pj, to_hd, idx, from_hd->north - from_hd->ns_res / 2);
+ proj_update(from_pj, to_pj, trans_pj, dir, to_hd, idx, from_hd->north - from_hd->ns_res / 2);
debug("Top", to_hd);
@@ -106,7 +109,7 @@
/* Right */
for (idx = from_hd->north - from_hd->ns_res / 2; idx > from_hd->south;
idx -= from_hd->ns_res)
- proj_update(from_pj, to_pj, to_hd, from_hd->east - from_hd->ew_res / 2, idx);
+ proj_update(from_pj, to_pj, trans_pj, dir, to_hd, from_hd->east - from_hd->ew_res / 2, idx);
debug("Right", to_hd);
@@ -113,7 +116,7 @@
/* Bottom */
for (idx = from_hd->east - from_hd->ew_res / 2; idx > from_hd->west;
idx -= from_hd->ew_res)
- proj_update(from_pj, to_pj, to_hd, idx, from_hd->south + from_hd->ns_res / 2);
+ proj_update(from_pj, to_pj, trans_pj, dir, to_hd, idx, from_hd->south + from_hd->ns_res / 2);
debug("Bottom", to_hd);
@@ -120,20 +123,23 @@
/* Left */
for (idx = from_hd->south + from_hd->ns_res / 2; idx < from_hd->north;
idx += from_hd->ns_res)
- proj_update(from_pj, to_pj, to_hd, from_hd->west + from_hd->ew_res / 2, idx);
+ proj_update(from_pj, to_pj, trans_pj, dir, to_hd, from_hd->west + from_hd->ew_res / 2, idx);
debug("Left", to_hd);
}
static int proj_inside(const struct pj_info *from_pj, const struct pj_info *to_pj,
+ const struct pj_info *trans_pj, int dir,
const struct Cell_head *ref_hd, double hx, double hy)
{
- if (pj_do_proj(&hx, &hy, to_pj, from_pj) < 0)
+ if (GPJ_transform(from_pj, to_pj, trans_pj, -dir, &hx, &hy, NULL) < 0)
return 0;
return inside(ref_hd, hx, hy);
}
-static void reverse_check(const struct pj_info *from_pj, const struct pj_info *to_pj,
+static void reverse_check(const struct pj_info *from_pj,
+ const struct pj_info *to_pj,
+ const struct pj_info *trans_pj, int dir,
const struct Cell_head *from_hd,
const struct Cell_head *to_hd,
struct Cell_head *cur_hd)
@@ -141,7 +147,7 @@
if (cur_hd->west > to_hd->west) {
double hx = to_hd->west + to_hd->ew_res / 2;
double hy = to_hd->south + (to_hd->north - to_hd->south) / 2;
- if (proj_inside(from_pj, to_pj, from_hd, hx, hy))
+ if (proj_inside(from_pj, to_pj, trans_pj, dir, from_hd, hx, hy))
cur_hd->west = hx;
}
@@ -148,7 +154,7 @@
if (cur_hd->east < to_hd->east) {
double hx = to_hd->east - to_hd->ew_res / 2;
double hy = to_hd->south + (to_hd->north - to_hd->south) / 2;
- if (proj_inside(from_pj, to_pj, from_hd, hx, hy))
+ if (proj_inside(from_pj, to_pj, trans_pj, dir, from_hd, hx, hy))
cur_hd->east = hx;
}
@@ -155,7 +161,7 @@
if (cur_hd->south > to_hd->south) {
double hx = to_hd->west + (to_hd->east - to_hd->west) / 2;
double hy = to_hd->south + to_hd->ns_res / 2;
- if (proj_inside(from_pj, to_pj, from_hd, hx, hy))
+ if (proj_inside(from_pj, to_pj, trans_pj, dir, from_hd, hx, hy))
cur_hd->south = hy;
}
@@ -162,7 +168,7 @@
if (cur_hd->north < to_hd->north) {
double hx = to_hd->west + (to_hd->east - to_hd->west) / 2;
double hy = to_hd->north - to_hd->ns_res / 2;
- if (proj_inside(from_pj, to_pj, from_hd, hx, hy))
+ if (proj_inside(from_pj, to_pj, trans_pj, dir, from_hd, hx, hy))
cur_hd->north = hy;
}
}
@@ -190,7 +196,8 @@
}
void bordwalk(const struct Cell_head *from_hd, struct Cell_head *to_hd,
- const struct pj_info *from_pj, const struct pj_info *to_pj)
+ const struct pj_info *from_pj, const struct pj_info *to_pj,
+ const struct pj_info *trans_pj, int dir)
{
struct Cell_head cur_hd;
@@ -200,13 +207,13 @@
/* Start walking */
- bordwalk1(from_pj, to_pj, from_hd, &cur_hd);
+ bordwalk1(from_pj, to_pj, trans_pj, dir, from_hd, &cur_hd);
intersect(&cur_hd, to_hd);
/* check some special cases by reversing the projection */
- reverse_check(from_pj, to_pj, from_hd, to_hd, &cur_hd);
+ reverse_check(from_pj, to_pj, trans_pj, dir, from_hd, to_hd, &cur_hd);
debug("Extra check", &cur_hd);
@@ -227,7 +234,8 @@
}
void bordwalk_edge(const struct Cell_head *from_hd, struct Cell_head *to_hd,
- const struct pj_info *from_pj, const struct pj_info *to_pj)
+ const struct pj_info *from_pj, const struct pj_info *to_pj,
+ const struct pj_info *trans_pj, int dir)
{
double idx;
double hx, hy;
@@ -238,7 +246,8 @@
hx = (from_hd->west + from_hd->east) / 2.0;
hy = (from_hd->north + from_hd->south) / 2.0;
- if (pj_do_proj(&hx, &hy, from_pj, to_pj) < 0)
+ if (GPJ_transform(from_pj, to_pj, trans_pj, dir,
+ &hx, &hy, NULL) < 0)
G_fatal_error(_("Unable to reproject map center"));
to_hd->east = hx;
@@ -249,9 +258,9 @@
/* Top */
for (idx = from_hd->west; idx < from_hd->east;
idx += from_hd->ew_res)
- proj_update(from_pj, to_pj, to_hd, idx, from_hd->north);
+ proj_update(from_pj, to_pj, trans_pj, dir, to_hd, idx, from_hd->north);
idx = from_hd->east;
- proj_update(from_pj, to_pj, to_hd, idx, from_hd->north);
+ proj_update(from_pj, to_pj, trans_pj, dir, to_hd, idx, from_hd->north);
debug("Top", to_hd);
@@ -258,9 +267,9 @@
/* Right */
for (idx = from_hd->north; idx > from_hd->south;
idx -= from_hd->ns_res)
- proj_update(from_pj, to_pj, to_hd, from_hd->east, idx);
+ proj_update(from_pj, to_pj, trans_pj, dir, to_hd, from_hd->east, idx);
idx = from_hd->south;
- proj_update(from_pj, to_pj, to_hd, from_hd->east, idx);
+ proj_update(from_pj, to_pj, trans_pj, dir, to_hd, from_hd->east, idx);
debug("Right", to_hd);
@@ -267,9 +276,9 @@
/* Bottom */
for (idx = from_hd->east; idx > from_hd->west;
idx -= from_hd->ew_res)
- proj_update(from_pj, to_pj, to_hd, idx, from_hd->south);
+ proj_update(from_pj, to_pj, trans_pj, dir, to_hd, idx, from_hd->south);
idx = from_hd->west;
- proj_update(from_pj, to_pj, to_hd, idx, from_hd->south);
+ proj_update(from_pj, to_pj, trans_pj, dir, to_hd, idx, from_hd->south);
debug("Bottom", to_hd);
@@ -276,9 +285,9 @@
/* Left */
for (idx = from_hd->south; idx < from_hd->north;
idx += from_hd->ns_res)
- proj_update(from_pj, to_pj, to_hd, from_hd->west, idx);
+ proj_update(from_pj, to_pj, trans_pj, dir, to_hd, from_hd->west, idx);
idx = from_hd->north;
- proj_update(from_pj, to_pj, to_hd, from_hd->west, idx);
+ proj_update(from_pj, to_pj, trans_pj, dir, to_hd, from_hd->west, idx);
debug("Left", to_hd);
}
Modified: grass/trunk/raster/r.proj/main.c
===================================================================
--- grass/trunk/raster/r.proj/main.c 2018-03-23 09:25:04 UTC (rev 72518)
+++ grass/trunk/raster/r.proj/main.c 2018-03-23 09:25:34 UTC (rev 72519)
@@ -112,7 +112,8 @@
struct History history;
struct pj_info iproj, /* input map proj parameters */
- oproj; /* output map proj parameters */
+ oproj, /* output map proj parameters */
+ tproj; /* transformation parameters */
struct Key_Value *in_proj_info, /* projection information of */
*in_unit_info, /* input and output mapsets */
@@ -316,6 +317,9 @@
if (pj_get_kv(&iproj, in_proj_info, in_unit_info) < 0)
G_fatal_error(_("Unable to get projection key values of input map"));
+ if (GPJ_init_transform(&iproj, &oproj, &tproj) < 0)
+ G_fatal_error(_("Unable to initialize coordinate transformation"));
+
G_free_key_value(in_proj_info);
G_free_key_value(in_unit_info);
G_free_key_value(out_proj_info);
@@ -355,7 +359,7 @@
outcellhd.south = 1e9;
outcellhd.east = -1e9;
outcellhd.west = 1e9;
- bordwalk_edge(&incellhd, &outcellhd, &iproj, &oproj);
+ bordwalk_edge(&incellhd, &outcellhd, &iproj, &oproj, &tproj, PJ_FWD);
inorth = outcellhd.north;
isouth = outcellhd.south;
ieast = outcellhd.east;
@@ -385,7 +389,7 @@
/* Cut non-overlapping parts of input map */
if (!nocrop->answer)
- bordwalk(&outcellhd, &incellhd, &oproj, &iproj);
+ bordwalk(&outcellhd, &incellhd, &iproj, &oproj, &tproj, PJ_INV);
/* Add 2 cells on each side for bilinear/cubic & future interpolation methods */
/* (should probably be a factor based on input and output resolution) */
@@ -411,7 +415,7 @@
/* Adjust borders of output map */
if (!nocrop->answer)
- bordwalk(&incellhd, &outcellhd, &iproj, &oproj);
+ bordwalk(&incellhd, &outcellhd, &iproj, &oproj, &tproj, PJ_FWD);
#if 0
outcellhd.west = outcellhd.south = HUGE_VAL;
@@ -420,7 +424,9 @@
ycoord1 = Rast_row_to_northing((double)(row + 0.5), &incellhd);
for (col = 0; col < incellhd.cols; col++) {
xcoord1 = Rast_col_to_easting((double)(col + 0.5), &incellhd);
- pj_do_proj(&xcoord1, &ycoord1, &iproj, &oproj);
+ if (GPJ_transform(&iproj, &oproj, &tproj, PJ_FWD,
+ &xcoord1, &ycoord1, NULL) < 0)
+ G_fatal_error(_("Error in %s"), "GPJ_transform()");
if (xcoord1 > outcellhd.east)
outcellhd.east = xcoord1;
if (ycoord1 > outcellhd.north)
@@ -509,8 +515,11 @@
/* project coordinates in output matrix to */
/* coordinates in input matrix */
- if (pj_do_proj(&xcoord1, &ycoord1, &oproj, &iproj) < 0)
+ if (GPJ_transform(&iproj, &oproj, &tproj, PJ_INV,
+ &xcoord1, &ycoord1, NULL) < 0) {
+ G_warning(_("Error in %s"), "GPJ_transform()");
Rast_set_null_value(obufptr, 1, cell_type);
+ }
else {
/* convert to row/column indices of input matrix */
Modified: grass/trunk/raster/r.proj/r.proj.h
===================================================================
--- grass/trunk/raster/r.proj/r.proj.h 2018-03-23 09:25:04 UTC (rev 72518)
+++ grass/trunk/raster/r.proj/r.proj.h 2018-03-23 09:25:34 UTC (rev 72519)
@@ -36,9 +36,11 @@
};
extern void bordwalk(const struct Cell_head *, struct Cell_head *,
- const struct pj_info *, const struct pj_info *);
+ const struct pj_info *, const struct pj_info *,
+ const struct pj_info *, int);
extern void bordwalk_edge(const struct Cell_head *, struct Cell_head *,
- const struct pj_info *, const struct pj_info *);
+ const struct pj_info *, const struct pj_info *,
+ const struct pj_info *, int);
extern struct cache *readcell(int, const char *);
extern block *get_block(struct cache *, int);
extern void release_cache(struct cache *);
More information about the grass-commit
mailing list