Tests fail with GDAL 3.0.2

[phargogh]

In a new virtualenv, I pip installed pygeoprocessing==1.9.1 using dependencies (including GDAL 3.0.2 from Gohlke) from our PyPI. 3 of our tests fail. I've attached the log from this test run here: gdal3_error_log.txt

Tests work great on the latest version of GDAL v2 (currently 2.4.2, I believe)

[phargogh]

It's also worth noting that Gohlke isn't producing Python 3.8 wheels for GDAL 2: https://www.lfd.uci.edu/~gohlke/pythonlibs/#gdal, so GDAL3 should really be in our future.

[dcdenu4]

Found some interesting things going on with difference between GDAL 3 and GDAL 2.
When transforming bounding box from align_and_resize_raster_stack:
GDAL 2 [Python 3.7]:
target_bbox : [845079.5702749295, -99734.04483739965, 956619.8421206995, 11081.556155365919]
GDAL 3 [Python 3.8]:
target_bbox : [733729.7578060223, 11060.483077897718, 845058.6446509969, 121762.53706604337]

Wanted to document this before signing off for the day.

[dcdenu4]

There was some pretty big changes to GDAL 3.0 relating to PROJ 6. Here are some notes: proj6-wkt2. This will take awhile to digest.

[dcdenu4]

It looks like the change with the axis order is the issue. Switching the points in pygeo _transform_point from x, y to y, x gives the tested expected results. Supposedly setting wgs84_sr.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) is supposed to keep the previous order, but I have not gotten it to work yet.

This post sums it up quite well: OSGeo/gdal#1546

[dcdenu4]

Came across a fresh error in tests for GDAL 3+

__________________________________________________________________________________ UFRMTests.test_ufrm_regression ___________________________________________________________________________________

self = <test_ufrm.UFRMTests testMethod=test_ufrm_regression>

    def test_ufrm_regression(self):
        """UFRM: regression test."""
        from natcap.invest import urban_flood_risk_mitigation
        args = self._make_args()
        urban_flood_risk_mitigation.execute(args)

        result_vector = gdal.OpenEx(os.path.join(
            args['workspace_dir'], 'flood_risk_service_Test1.shp'),
            gdal.OF_VECTOR)
        result_layer = result_vector.GetLayer()
>       result_feature = next(result_layer)
E       TypeError: 'Layer' object is not an iterator

tests\test_ufrm.py:63: TypeError

[dcdenu4]

Instead of using next it might be better to call: GetNextFeature() from : https://gdal.org/python/osgeo.ogr.Layer-class.html#GetNextFeature. This allows the test to pass in GDAL 3.

[dcdenu4]

GDAL 3 expected test outputs will vary slightly in some cases from GDAL 2. Should we support both GDAL 2 and GDAL 3 test cases?

[dcdenu4]

Below is a post I made to the premature PR for this issue:

I just spent a significant amount of time in the WKT2 docs, PROJ 6 docs, and GDAL 3 Docs. I think I'm starting to develop a deeper understanding. I think the basics are as follows:

• GDAL 3 conforms to the authority defining coordinate reference systems (CRS).
• For Geographic CRS (WGS84), this means no longer ignoring Axis Order but abiding by LAT,LON order
• For Projected CRS (UTM), most EPSG in the registry use Easting, Northing order
• We think of things in general as x,y because it is convenient, but now a Transform expects Axis Order, where if long=x and lat=y, GDAL 3 would take Transform(y, x)
• Thus when transforming from Projected CRS to Geographic CRS, Transform(x,y) will return y,x (Lat,Lon)
• Thus when transforming from Geographic CRS to Projected CRS, Transform(y,x) will return x,y (Easting, Northing)
• Thus when transforming from Geographic CRS to Geographic CRS, Transform(y,x) will return y,x (Lat,Lon)
• Thus when transforming from Projected CRS to Projected CRS, Transform(x,y) will return x,y (Easting, Northing)

note: Most projected CRS use Easting first, Northing second convention, however, some defined in the EPSG registry use the reverse convention.

So, I think you're correct in that not.____.IsProjected() is probably wrong, but _____.IsGeographic() could be OK. Basically that Axis Order setting SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) is saying, if we've got a Geographic CRS, use Lon,Lat like the good old days (since we think of lon,lat as x,y. We certainly can, but maybe should, NOT do that.

Perhaps a better solution is to create a mapping of Axis Order to Geographic (Lat,Lon) and Projected (Easting, Northing), such that when handling the extents and bounding box points we can correctly return x, y.

[dcdenu4]

I'm closing this issue since the GDAL 3 PR has been merged and since there is another issue #44 that reflects why the tests are failing.