From a0f1812c7b0492d1d5a3d8eb3e50f30f56b6253a Mon Sep 17 00:00:00 2001
From: Mike Taves <mwtoews@gmail.com>
Date: Tue, 31 Mar 2020 22:36:56 +1300
Subject: [PATCH] feat(netcdf): use modern features from pyproj>=2.2.0

* For newer pyproj, use CRS and Transformer classes, to take advantage
  of modern transformation pipelines to convert coordinates
* For older pyproj, use Proj class, which has limitations
* Remove a few unnecessary try/except blocks and repeated log messages
---
 flopy/export/netcdf.py | 57 ++++++++++++++++++++++++------------------
 1 file changed, 32 insertions(+), 25 deletions(-)

diff --git a/flopy/export/netcdf.py b/flopy/export/netcdf.py
index 042855ffc0..66e211842d 100644
--- a/flopy/export/netcdf.py
+++ b/flopy/export/netcdf.py
@@ -624,25 +624,26 @@ def initialize_geometry(self):
             needed for the netcdf file
         """
         try:
-            from pyproj import Proj, transform
-        except Exception as e:
-            raise Exception("NetCdf error importing pyproj module:\n" + str(e))
+            import pyproj
+        except ImportError as e:
+            raise ImportError(
+                "NetCdf error importing pyproj module:\n" + str(e))
+        from distutils.version import LooseVersion
+
+        # Check if using newer pyproj version conventions
+        pyproj220 = LooseVersion(pyproj.__version__) >= LooseVersion('2.2.0')
 
         proj4_str = self.proj4_str
         print('initialize_geometry::proj4_str = {}'.format(proj4_str))
 
         self.log("building grid crs using proj4 string: {}".format(proj4_str))
-        try:
-            self.grid_crs = Proj(proj4_str, preserve_units=True)
-
-        except Exception as e:
-            self.log("error building grid crs:\n{0}".format(str(e)))
-            raise Exception("error building grid crs:\n{0}".format(str(e)))
+        if pyproj220:
+            self.grid_crs = pyproj.CRS(proj4_str)
+        else:
+            self.grid_crs = pyproj.Proj(proj4_str, preserve_units=True)
 
         print('initialize_geometry::self.grid_crs = {}'.format(self.grid_crs))
 
-        self.log("building grid crs using proj4 string: {}".format(proj4_str))
-
         vmin, vmax = self.model_grid.botm.min(), \
                      self.model_grid.top.max()
         if self.z_positive == 'down':
@@ -654,21 +655,25 @@ def initialize_geometry(self):
         xs = self.model_grid.xyzcellcenters[0].copy()
 
         # Transform to a known CRS
-        nc_crs = Proj(self.nc_epsg_str)
+        if pyproj220:
+            nc_crs = pyproj.CRS(self.nc_epsg_str)
+            self.transformer = pyproj.Transformer.from_crs(
+                self.grid_crs, nc_crs, always_xy=True)
+        else:
+            nc_crs = pyproj.Proj(self.nc_epsg_str)
+            self.transformer = None
+
         print('initialize_geometry::nc_crs = {}'.format(nc_crs))
 
-        self.log("projecting grid cell center arrays " + \
-                 "from {} to {}".format(str(self.grid_crs.srs),
-                                        str(nc_crs.srs)))
-        try:
-            self.xs, self.ys = transform(self.grid_crs, nc_crs, xs, ys)
-        except Exception as e:
-            self.log("error projecting:\n{0}".format(str(e)))
-            raise Exception("error projecting:\n{0}".format(str(e)))
+        if pyproj220:
+            print('transforming coordinates using = {}'
+                  .format(self.transformer))
 
-        self.log("projecting grid cell center arrays " + \
-                 "from {0} to {1}".format(str(self.grid_crs),
-                                          str(nc_crs)))
+        self.log("projecting grid cell center arrays")
+        if pyproj220:
+            self.xs, self.ys = self.transformer.transform(xs, ys)
+        else:
+            self.xs, self.ys = pyproj.transform(self.grid_crs, nc_crs, xs, ys)
 
         # get transformed bounds and record to check against ScienceBase later
         xmin, xmax, ymin, ymax = self.model_grid.extent
@@ -676,10 +681,12 @@ def initialize_geometry(self):
                          [xmin, ymax],
                          [xmax, ymax],
                          [xmax, ymin]])
-        x, y = transform(self.grid_crs, nc_crs, *bbox.transpose())
+        if pyproj220:
+            x, y = self.transformer.transform(*bbox.transpose())
+        else:
+            x, y = pyproj.transform(self.grid_crs, nc_crs, *bbox.transpose())
         self.bounds = x.min(), y.min(), x.max(), y.max()
         self.vbounds = vmin, vmax
-        pass
 
     def initialize_file(self, time_values=None):
         """