libproj: fix area bbox for PROJ (#2467)

* fix area bbox for PROJ * add tests for valid area bbox * test if crossing the anti-meridian
OSGeo · Jun 30, 2022 · da18254 · da18254
1 parent 6c1c13c
commit da18254
Show file tree

Hide file tree

Showing 2 changed files with 50 additions and 6 deletions.
diff --git a/lib/proj/do_proj.c b/lib/proj/do_proj.c
@@ -17,6 +17,7 @@
 #include <stdio.h>
 #include <string.h>
 #include <ctype.h>
+#include <math.h>
 
 #include <grass/gis.h>
 #include <grass/gprojects.h>
@@ -118,8 +119,9 @@ int get_pj_area(const struct pj_info *iproj, double *xmin, double *xmax,
 	}
 	G_free(indef);
 
-	estep = (window.west + window.east) / 21.;
-	nstep = (window.north + window.south) / 21.;
+	/* inpired by gdal/ogr/ogrct.cpp OGRProjCT::ListCoordinateOperations() */
+	estep = (window.east - window.west) / 21.;
+	nstep = (window.north - window.south) / 21.;
 	for (i = 0; i < 20; i++) {
 	    x[i] = window.west + estep * (i + 1);
 	    y[i] = window.north;
@@ -162,6 +164,17 @@ int get_pj_area(const struct pj_info *iproj, double *xmin, double *xmax,
 		*ymax = y[i];
 	}
 
+	/* The west longitude is generally lower than the east longitude,
+	 * except for areas of interest that go across the anti-meridian.
+	 */
+	if (x[80] > x[82] && x[81] > x[83]) {
+	    /* west > east, must be crossing the anti-meridian */
+	    double tmpx = *xmin;
+
+	    *xmin = *xmax;
+	    *xmax = tmpx;
+	}
+
 	G_debug(1, "input window north: %.8f", window.north);
 	G_debug(1, "input window south: %.8f", window.south);
 	G_debug(1, "input window east: %.8f", window.east);
@@ -171,6 +184,30 @@ int get_pj_area(const struct pj_info *iproj, double *xmin, double *xmax,
 	G_debug(1, "transformed xmax: %.8f", *xmax);
 	G_debug(1, "transformed ymin: %.8f", *ymin);
 	G_debug(1, "transformed ymax: %.8f", *ymax);
+
+	/* test valid values, as in
+	 * gdal/ogr/ogrct.cpp OGRCoordinateTransformationOptions::SetAreaOfInterest()
+	 */
+	if (fabs(*xmin) > 180) {
+	    G_warning(_("Invalid west longitude %g"), *xmin);
+	    return 0;
+	}
+	if (fabs(*xmax) > 180) {
+	    G_warning(_("Invalid east longitude %g"), *xmax);
+	    return 0;
+	}
+	if (fabs(*ymin) > 90) {
+	    G_warning(_("Invalid south latitude %g"), *ymin);
+	    return 0;
+	}
+	if (fabs(*ymax) > 90) {
+	    G_warning(_("Invalid north latitude %g"), *ymax);
+	    return 0;
+	}
+	if (*ymin > *ymax) {
+	    G_warning(_("South %g is larger than north %g"), *ymin, *ymax);
+	    return 0;
+	}
     }
     G_debug(1, "get_pj_area(): xmin %g, xmax %g, ymin %g, ymax %g",
             *xmin, *xmax, *ymin, *ymax);
@@ -506,6 +543,11 @@ int GPJ_init_transform(const struct pj_info *info_in,
 	    pj_area = proj_area_create();
 	    proj_area_set_bbox(pj_area, xmin, ymin, xmax, ymax);
 	}
+	else {
+	    G_warning(_("Unable to determine area of interest for '%s'"), info_in->def);
+
+	    return -1;
+	}
 
 	G_debug(1, "source proj string: %s", info_in->def);
 	G_debug(1, "source type: %s", get_pj_type_string(info_in->pj));
@@ -692,7 +734,9 @@ int GPJ_init_transform(const struct pj_info *info_in,
 	    proj_list_destroy(op_list);
 
 	/* try proj_create_crs_to_crs() */
+	/*
 	G_debug(1, "trying %s to %s", indef, outdef);
+	*/
 
 	/* proj_create_crs_to_crs() does not work because it calls
 	 * proj_create_crs_to_crs_from_pj() which calls

diff --git a/vector/v.proj/main.c b/vector/v.proj/main.c
@@ -201,10 +201,6 @@ int main(int argc, char *argv[])
     info_out.srid = G_get_projsrid();
     info_out.wkt = G_get_projwkt();
 
-    if (G_verbose() == G_verbose_max()) {
-	pj_print_proj_params(&info_in, &info_out);
-    }
-
     info_trans.def = NULL;
 #ifdef HAVE_PROJ_H
     if (pipeline->answer) {
@@ -279,6 +275,10 @@ int main(int argc, char *argv[])
 	info_in.srid = G_get_projsrid();
 	info_in.wkt = G_get_projwkt();
 
+	if (G_verbose() == G_verbose_max()) {
+	    pj_print_proj_params(&info_in, &info_out);
+	}
+
 	Vect_set_open_level(1);
 	G_debug(1, "Open old: location: %s mapset : %s", G_location_path(),
 		G_mapset());