Update issue natcap#27 adding coordinate transform func

Adding coordinate transformation function to utils.py so that all InVEST
models that handle their own transformations can now have a single
entry point for it. This function is compatible with gdal 3 but
NOT gdal 2. This function sets unprojected spatial references to
be "GIS Friendly", meaning use lon,lat when transforming points.

Adding tests for this new function.

src/natcap/invest/utils.py

@@ -539,3 +539,33 @@ def mean_pixel_size_and_area(pixel_size_tuple):
                 pixel_size_tuple))
 
     return (x_size, x_size*y_size)
+
+
+def create_coordinate_transformer(base_ref, target_ref):
+    """Create a spatial reference coordinate transformation function.
+
+    The function creates a transformer that is compatable with gdal 3, 
+    that uses the gdal 2 to gdal 3 migration recommendation from: 
+    https://trac.osgeo.org/gdal/wiki/rfc73_proj6_wkt2_srsbarn.
+
+    Parameter:
+        base_ref (osr spatial reference): A defined spatial reference to 
+            transform FROM
+        target_ref (osr spatial reference): A defined spatial reference 
+            to transform TO
+
+    Returns:
+        An OSR Coordinate Transformation object
+
+    """
+    # GDAL 3 handles lat/lon transformations differently where the transformer
+    # expects lat,lon instead of lon,lat. GDAL migration help recommends 
+    # using SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) as a 
+    # workaround. Only need to worry about this for unprojected references.
+    if not base_ref.IsProjected():
+            base_ref.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) 
+    if not target_ref.IsProjected():
+            target_ref.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) 
+
+    transformer = osr.CreateCoordinateTransformation(base_ref, target_ref)
+    return transformer

tests/test_utils.py

@@ -14,6 +14,7 @@
 from pygeoprocessing.testing import scm
 import pygeoprocessing.testing
 from osgeo import gdal
+from osgeo import osr
 
 
 class SuffixUtilsTests(unittest.TestCase):
@@ -621,3 +622,103 @@ def test_csv_latin_1_encoding(self):
         self.assertEqual(lookup_dict[4]['header 2'], 5)
         self.assertEqual(lookup_dict[4]['header 3'], 'foo')
         self.assertEqual(lookup_dict[1]['header 1'], 1)
+
+
+class CreateCoordinateTransformationTests(unittest.TestCase):
+    """Tests for natcap.invest.utils.create_coordinate_transformer."""
+
+    def test_latlon_to_latlon_transformer(self):
+        """Utils: test transformer for lat/lon to lat/lon"""
+        from natcap.invest import utils
+
+        # Willamette valley in lat/lon for reference
+        bounding_box = [-123.587984, 44.415778, -123.397976, 44.725814]
+        lon = -124.525
+        lat = 44.525
+
+        base_srs = osr.SpatialReference()
+        base_srs.ImportFromEPSG(4326) # WSG84 EPSG
+
+        target_srs = osr.SpatialReference()
+        target_srs.ImportFromEPSG(4326)
+
+        transformer = utils.create_coordinate_transformer(base_srs, target_srs)
+        x, y, _ = transformer.TransformPoint(lon, lat)
+
+        expected_x = -124.525 
+        expected_y = 44.525
+
+        self.assertEqual(expected_x, x)
+        self.assertEqual(expected_y, y)
+
+    def test_latlon_to_projected_transformer(self):
+        """Utils: test transformer for lat/lon to projected."""
+        from natcap.invest import utils
+
+        # Willamette valley in lat/lon for reference
+        bounding_box = [-123.587984, 44.415778, -123.397976, 44.725814]
+        lon = -124.525
+        lat = 44.525
+
+        base_srs = osr.SpatialReference()
+        base_srs.ImportFromEPSG(4326) # WSG84 EPSG
+
+        target_srs = osr.SpatialReference()
+        target_srs.ImportFromEPSG(26910) # UTM10N EPSG
+
+        transformer = utils.create_coordinate_transformer(base_srs, target_srs)
+        x, y, _ = transformer.TransformPoint(lon, lat)
+
+        expected_x = 378816.2531852932 
+        expected_y = 4931317.807472325 
+
+        self.assertEqual(expected_x, x)
+        self.assertEqual(expected_y, y)
+
+    def test_projected_to_latlon_transformer(self):
+        """Utils: test transformer for projected to lat/lon."""
+        from natcap.invest import utils
+
+        # Willamette valley in lat/lon for reference
+        bounding_box = [-123.587984, 44.415778, -123.397976, 44.725814]
+        known_x = 378816.2531852932 
+        known_y = 4931317.807472325 
+
+        base_srs = osr.SpatialReference()
+        base_srs.ImportFromEPSG(26910) # UTM10N EPSG
+
+        target_srs = osr.SpatialReference()
+        target_srs.ImportFromEPSG(4326) # WSG84 EPSG
+
+        transformer = utils.create_coordinate_transformer(base_srs, target_srs)
+        x, y, _ = transformer.TransformPoint(known_x, known_y)
+
+        expected_x = -124.52500000000002
+        expected_y = 44.525 
+
+        self.assertEqual(expected_x, x)
+        self.assertEqual(expected_y, y)
+
+    def test_projected_to_projected_transformer(self):
+        """Utils: test transformer for projected to projected."""
+        from natcap.invest import utils
+
+        # Willamette valley in lat/lon for reference
+        bounding_box = [-123.587984, 44.415778, -123.397976, 44.725814]
+        known_x = 378816.2531852932 
+        known_y = 4931317.807472325 
+
+        base_srs = osr.SpatialReference()
+        base_srs.ImportFromEPSG(26910) # UTM10N EPSG
+
+        target_srs = osr.SpatialReference()
+        target_srs.ImportFromEPSG(26910) # UTM10N EPSG
+
+        transformer = utils.create_coordinate_transformer(base_srs, target_srs)
+        x, y, _ = transformer.TransformPoint(known_x, known_y)
+
+        expected_x = 378816.2531852932 
+        expected_y = 4931317.807472325 
+
+        self.assertEqual(expected_x, x)
+        self.assertEqual(expected_y, y)