Skip to content

Commit

Permalink
feat(netcdf): use modern features from pyproj>=2.2.0 (#840)
Browse files Browse the repository at this point in the history
* 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
  • Loading branch information
mwtoews committed Mar 31, 2020
1 parent abef931 commit fb942b4
Showing 1 changed file with 32 additions and 25 deletions.
57 changes: 32 additions & 25 deletions flopy/export/netcdf.py
Original file line number Diff line number Diff line change
Expand Up @@ -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':
Expand All @@ -654,32 +655,38 @@ 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
bbox = np.array([[xmin, ymin],
[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):
"""
Expand Down

0 comments on commit fb942b4

Please sign in to comment.