r.proj: use new GRASS API for coordinate transformation

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@72519 15284696-431f-4ddb-bdfa-cd5b030d7da7

raster/r.proj/bordwalk.c

@@ -84,14 +84,17 @@ static void invert(struct Cell_head *cur_hd, const struct Cell_head *ref_hd,
 }
 
 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,41 +102,44 @@ void bordwalk1(const struct pj_info *from_pj, const struct pj_info *to_pj,
     /* 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);
 
     /* 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);
 
     /* 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);
 
     /* 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,28 +147,28 @@ static void reverse_check(const struct pj_info *from_pj, const struct pj_info *t
     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;
     }
 
     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;
     }
 
     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;
     }
 
     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 @@ static void snap_to_grid(struct Cell_head *cur_hd, const struct Cell_head *ref_h
 }
 
 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 @@ void bordwalk(const struct Cell_head *from_hd, struct Cell_head *to_hd,
 
     /* 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(const struct Cell_head *from_hd, struct Cell_head *to_hd,
 }
 
 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 @@ void bordwalk_edge(const struct Cell_head *from_hd, struct Cell_head *to_hd,
     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,36 +258,36 @@ void bordwalk_edge(const struct Cell_head *from_hd, struct Cell_head *to_hd,
     /* 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);
 
     /* 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);
 
     /* 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);
 
     /* 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);
 }

raster/r.proj/main.c

@@ -112,7 +112,8 @@ int main(int argc, char **argv)
     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 @@ int main(int argc, char **argv)
     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 @@ int main(int argc, char **argv)
 	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 @@ int main(int argc, char **argv)
 
     /* 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 @@ int main(int argc, char **argv)
     /* 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 @@ int main(int argc, char **argv)
 	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 @@ int main(int argc, char **argv)
 
 	    /* 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 */
 

raster/r.proj/r.proj.h

@@ -36,9 +36,11 @@ struct menu
 };
 
 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 *);