import os
import platform
import socket
import copy
import json
import numpy as np
from datetime import datetime
import time
from .metadata import acdd
import flopy
# globals
FILLVALUE = -99999.9
ITMUNI = {
0: "undefined",
1: "seconds",
2: "minutes",
3: "hours",
4: "days",
5: "years",
}
PRECISION_STRS = ["f4", "f8", "i4"]
STANDARD_VARS = ["longitude", "latitude", "layer", "elevation", "time"]
path = os.path.split(__file__)[0]
with open(path + "/longnames.json") as f:
NC_LONG_NAMES = json.load(f)
[docs]class Logger(object):
"""
Basic class for logging events during the linear analysis calculations
if filename is passed, then an file handle is opened
Parameters
----------
filename : bool or string
if string, it is the log file to write. If a bool, then log is
written to the screen. echo (bool): a flag to force screen output
Attributes
----------
items : dict
tracks when something is started. If a log entry is
not in items, then it is treated as a new entry with the string
being the key and the datetime as the value. If a log entry is
in items, then the end time and delta time are written and
the item is popped from the keys
"""
def __init__(self, filename, echo=False):
self.items = {}
self.echo = bool(echo)
if filename == True:
self.echo = True
self.filename = None
elif filename:
self.f = open(filename, "w", 0) # unbuffered
self.t = datetime.now()
self.log("opening " + str(filename) + " for logging")
else:
self.filename = None
[docs] def log(self, phrase):
"""
log something that happened
Parameters
----------
phrase : str
the thing that happened
"""
pass
t = datetime.now()
if phrase in self.items.keys():
s = (
str(t)
+ " finished: "
+ str(phrase)
+ ", took: "
+ str(t - self.items[phrase])
+ "\n"
)
if self.echo:
print(
s,
)
if self.filename:
self.f.write(s)
self.items.pop(phrase)
else:
s = str(t) + " starting: " + str(phrase) + "\n"
if self.echo:
print(
s,
)
if self.filename:
self.f.write(s)
self.items[phrase] = copy.deepcopy(t)
[docs] def warn(self, message):
"""
Write a warning to the log file
Parameters
----------
message : str
the warning text
"""
s = str(datetime.now()) + " WARNING: " + message + "\n"
if self.echo:
print(
s,
)
if self.filename:
self.f.write(s)
return
[docs]class NetCdf(object):
"""
Support for writing a netCDF4 compliant file from a flopy model
Parameters
----------
output_filename : str
Name of the .nc file to write
model : flopy model instance
time_values : the entries for the time dimension
if not None, the constructor will initialize
the file. If None, the perlen array of ModflowDis
will be used
z_positive : str ('up' or 'down')
Positive direction of vertical coordinates written to NetCDF file.
(default 'down')
verbose : if True, stdout is verbose. If str, then a log file
is written to the verbose file
forgive : what to do if a duplicate variable name is being created. If
True, then the newly requested var is skipped. If False, then
an exception is raised.
**kwargs : keyword arguments
modelgrid : flopy.discretization.Grid instance
user supplied model grid which will be used in lieu of the model
object modelgrid for netcdf production
Notes
-----
This class relies heavily on the grid and modeltime objects,
including these attributes: lenuni, itmuni, start_datetime, and proj4.
Make sure these attributes have meaningful values.
"""
def __init__(
self,
output_filename,
model,
time_values=None,
z_positive="up",
verbose=None,
prj=None,
logger=None,
forgive=False,
**kwargs
):
assert output_filename.lower().endswith(".nc")
if verbose is None:
verbose = model.verbose
if logger is not None:
self.logger = logger
else:
self.logger = Logger(verbose)
self.var_attr_dict = {}
self.log = self.logger.log
if os.path.exists(output_filename):
self.logger.warn("removing existing nc file: " + output_filename)
os.remove(output_filename)
self.output_filename = output_filename
self.forgive = bool(forgive)
self.model = model
self.model_grid = model.modelgrid
if "modelgrid" in kwargs:
self.model_grid = kwargs.pop("modelgrid")
self.model_time = model.modeltime
if prj is not None:
self.model_grid.proj4 = prj
if self.model_grid.grid_type == "structured":
self.dimension_names = ("layer", "y", "x")
STANDARD_VARS.extend(["delc", "delr"])
# elif self.model_grid.grid_type == 'vertex':
# self.dimension_names = ('layer', 'ncpl')
else:
raise Exception(
"Grid type {} not supported.".format(self.model_grid.grid_type)
)
self.shape = self.model_grid.shape
try:
import dateutil.parser
except:
print(
"python-dateutil is not installed\n"
+ "try pip install python-dateutil"
)
return
self.start_datetime = self._dt_str(
dateutil.parser.parse(self.model_time.start_datetime)
)
self.logger.warn("start datetime:{0}".format(str(self.start_datetime)))
proj4_str = self.model_grid.proj4
if proj4_str is None:
proj4_str = "epsg:4326"
self.log(
"Warning: model has no coordinate reference system specified. "
"Using default proj4 string: {}".format(proj4_str)
)
self.proj4_str = proj4_str
self.grid_units = self.model_grid.units
self.z_positive = z_positive
if self.grid_units is None:
self.grid_units = "undefined"
assert self.grid_units in ["feet", "meters", "undefined"], (
"unsupported length units: " + self.grid_units
)
self.time_units = self.model_time.time_units
# this gives us confidence that every NetCdf instance
# has the same attributes
self.log("initializing attributes")
self._initialize_attributes()
self.log("initializing attributes")
self.time_values_arg = time_values
self.log("initializing file")
self.initialize_file(time_values=self.time_values_arg)
self.log("initializing file")
def __add__(self, other):
new_net = NetCdf.zeros_like(self)
if np.isscalar(other) or isinstance(other, np.ndarray):
for vname in self.var_attr_dict.keys():
new_net.nc.variables[vname][:] = (
self.nc.variables[vname][:] + other
)
elif isinstance(other, NetCdf):
for vname in self.var_attr_dict.keys():
new_net.nc.variables[vname][:] = (
self.nc.variables[vname][:] + other.nc.variables[vname][:]
)
else:
raise Exception(
"NetCdf.__add__(): unrecognized other:{0}".format(
str(type(other))
)
)
return new_net
def __sub__(self, other):
new_net = NetCdf.zeros_like(self)
if np.isscalar(other) or isinstance(other, np.ndarray):
for vname in self.var_attr_dict.keys():
new_net.nc.variables[vname][:] = (
self.nc.variables[vname][:] - other
)
elif isinstance(other, NetCdf):
for vname in self.var_attr_dict.keys():
new_net.nc.variables[vname][:] = (
self.nc.variables[vname][:] - other.nc.variables[vname][:]
)
else:
raise Exception(
"NetCdf.__sub__(): unrecognized other:{0}".format(
str(type(other))
)
)
return new_net
def __mul__(self, other):
new_net = NetCdf.zeros_like(self)
if np.isscalar(other) or isinstance(other, np.ndarray):
for vname in self.var_attr_dict.keys():
new_net.nc.variables[vname][:] = (
self.nc.variables[vname][:] * other
)
elif isinstance(other, NetCdf):
for vname in self.var_attr_dict.keys():
new_net.nc.variables[vname][:] = (
self.nc.variables[vname][:] * other.nc.variables[vname][:]
)
else:
raise Exception(
"NetCdf.__mul__(): unrecognized other:{0}".format(
str(type(other))
)
)
return new_net
def __div__(self, other):
return self.__truediv__(other)
def __truediv__(self, other):
new_net = NetCdf.zeros_like(self)
with np.errstate(invalid="ignore"):
if np.isscalar(other) or isinstance(other, np.ndarray):
for vname in self.var_attr_dict.keys():
new_net.nc.variables[vname][:] = (
self.nc.variables[vname][:] / other
)
elif isinstance(other, NetCdf):
for vname in self.var_attr_dict.keys():
new_net.nc.variables[vname][:] = (
self.nc.variables[vname][:]
/ other.nc.variables[vname][:]
)
else:
raise Exception(
"NetCdf.__sub__(): unrecognized other:{0}".format(
str(type(other))
)
)
return new_net
[docs] def append(self, other, suffix="_1"):
assert isinstance(other, NetCdf) or isinstance(other, dict)
if isinstance(other, NetCdf):
for vname in other.var_attr_dict.keys():
attrs = other.var_attr_dict[vname].copy()
var = other.nc.variables[vname]
new_vname = vname
if vname in self.nc.variables.keys():
if vname not in STANDARD_VARS:
new_vname = vname + suffix
if "long_name" in attrs:
attrs["long_name"] += " " + suffix
else:
continue
assert (
new_vname not in self.nc.variables.keys()
), "var already exists:{0} in {1}".format(
new_vname, ",".join(self.nc.variables.keys())
)
attrs["max"] = var[:].max()
attrs["min"] = var[:].min()
new_var = self.create_variable(
new_vname, attrs, var.dtype, dimensions=var.dimensions
)
new_var[:] = var[:]
else:
for vname, array in other.items():
vname_norm = self.normalize_name(vname)
assert (
vname_norm in self.nc.variables.keys()
), "dict var not in " "self.vars:{0}-->".format(
vname
) + ",".join(
self.nc.variables.keys()
)
new_vname = vname_norm + suffix
assert new_vname not in self.nc.variables.keys()
attrs = self.var_attr_dict[vname_norm].copy()
attrs["max"] = np.nanmax(array)
attrs["min"] = np.nanmin(array)
attrs["name"] = new_vname
attrs["long_name"] = attrs["long_name"] + " " + suffix
var = self.nc.variables[vname_norm]
# assert var.shape == array.shape,\
# "{0} shape ({1}) doesn't make array shape ({2})".\
# format(new_vname,str(var.shape),str(array.shape))
new_var = self.create_variable(
new_vname, attrs, var.dtype, dimensions=var.dimensions
)
try:
new_var[:] = array
except:
new_var[:, 0] = array
return
[docs] def copy(self, output_filename):
new_net = NetCdf.zeros_like(self, output_filename=output_filename)
for vname in self.var_attr_dict.keys():
new_net.nc.variables[vname][:] = self.nc.variables[vname][:]
return new_net
[docs] @classmethod
def zeros_like(
cls, other, output_filename=None, verbose=None, logger=None
):
new_net = NetCdf.empty_like(
other,
output_filename=output_filename,
verbose=verbose,
logger=logger,
)
# add the vars to the instance
for vname in other.var_attr_dict.keys():
if new_net.nc.variables.get(vname) is not None:
new_net.logger.warn(
"variable {0} already defined, skipping".format(vname)
)
continue
new_net.log("adding variable {0}".format(vname))
var = other.nc.variables[vname]
data = var[:]
try:
mask = data.mask
data = np.array(data)
except:
mask = None
new_data = np.zeros_like(data)
new_data[mask] = FILLVALUE
new_var = new_net.create_variable(
vname,
other.var_attr_dict[vname],
var.dtype,
dimensions=var.dimensions,
)
new_var[:] = new_data
new_net.log("adding variable {0}".format(vname))
global_attrs = {}
for attr in other.nc.ncattrs():
if attr not in new_net.nc.ncattrs():
global_attrs[attr] = other.nc[attr]
new_net.add_global_attributes(global_attrs)
return new_net
[docs] @classmethod
def empty_like(
cls, other, output_filename=None, verbose=None, logger=None
):
if output_filename is None:
output_filename = (
str(time.mktime(datetime.now().timetuple())) + ".nc"
)
while os.path.exists(output_filename):
print("{}...already exists".format(output_filename))
output_filename = (
str(time.mktime(datetime.now().timetuple())) + ".nc"
)
print(
"creating temporary netcdf file..."
+ "{}".format(output_filename)
)
new_net = cls(
output_filename,
other.model,
time_values=other.time_values_arg,
verbose=verbose,
logger=logger,
)
return new_net
[docs] def difference(
self, other, minuend="self", mask_zero_diff=True, onlydiff=True
):
"""
make a new NetCDF instance that is the difference with another
netcdf file
Parameters
----------
other : either an str filename of a netcdf file or
a netCDF4 instance
minuend : (optional) the order of the difference operation.
Default is self (e.g. self - other). Can be "self" or "other"
mask_zero_diff : bool flag to mask differences that are zero. If
True, positions in the difference array that are zero will be set
to self.fillvalue
only_diff : bool flag to only add non-zero diffs to output file
Returns
-------
net NetCDF instance
Notes
-----
assumes the current NetCDF instance has been populated. The
variable names and dimensions between the two files must match
exactly. The name of the new .nc file is
<self.output_filename>.diff.nc. The masks from both self and
other are carried through to the new instance
"""
assert self.nc is not None, (
"can't call difference() if nc " + "hasn't been populated"
)
try:
import netCDF4
except Exception as e:
mess = "error import netCDF4: {0}".format(str(e))
self.logger.warn(mess)
raise Exception(mess)
if isinstance(other, str):
assert os.path.exists(
other
), "filename 'other' not found:" + "{0}".format(other)
other = netCDF4.Dataset(other, "r")
assert isinstance(other, netCDF4.Dataset)
# check for similar variables
self_vars = set(self.nc.variables.keys())
other_vars = set(other.variables)
diff = self_vars.symmetric_difference(other_vars)
if len(diff) > 0:
self.logger.warn(
"variables are not the same between the two "
+ "nc files: "
+ ",".join(diff)
)
return
# check for similar dimensions
self_dimens = self.nc.dimensions
other_dimens = other.dimensions
for d in self_dimens.keys():
if d not in other_dimens:
self.logger.warn("missing dimension in other:{0}".format(d))
return
if len(self_dimens[d]) != len(other_dimens[d]):
self.logger.warn(
"dimension not consistent: "
+ "{0}:{1}".format(self_dimens[d], other_dimens[d])
)
return
# should be good to go
time_values = self.nc.variables.get("time")[:]
new_net = NetCdf(
self.output_filename.replace(".nc", ".diff.nc"),
self.model,
time_values=time_values,
)
# add the vars to the instance
for vname in self_vars:
if (
vname not in self.var_attr_dict
or new_net.nc.variables.get(vname) is not None
):
self.logger.warn("skipping variable: {0}".format(vname))
continue
self.log("processing variable {0}".format(vname))
s_var = self.nc.variables[vname]
o_var = other.variables[vname]
s_data = s_var[:]
o_data = o_var[:]
o_mask, s_mask = None, None
# keep the masks to apply later
if isinstance(s_data, np.ma.MaskedArray):
self.logger.warn("masked array for {0}".format(vname))
s_mask = s_data.mask
s_data = np.array(s_data)
s_data[s_mask] = 0.0
else:
np.nan_to_num(s_data)
if isinstance(o_data, np.ma.MaskedArray):
o_mask = o_data.mask
o_data = np.array(o_data)
o_data[o_mask] = 0.0
else:
np.nan_to_num(o_data)
# difference with self
if minuend.lower() == "self":
d_data = s_data - o_data
elif minuend.lower() == "other":
d_data = o_data - s_data
else:
mess = "unrecognized minuend {0}".format(minuend)
self.logger.warn(mess)
raise Exception(mess)
# check for non-zero diffs
if onlydiff and d_data.sum() == 0.0:
self.logger.warn(
"var {0} has zero differences, skipping...".format(vname)
)
continue
self.logger.warn(
"resetting diff attrs max,min:{0},{1}".format(
d_data.min(), d_data.max()
)
)
attrs = self.var_attr_dict[vname].copy()
attrs["max"] = np.nanmax(d_data)
attrs["min"] = np.nanmin(d_data)
# reapply masks
if s_mask is not None:
self.log("applying self mask")
s_mask[d_data != 0.0] = False
d_data[s_mask] = FILLVALUE
self.log("applying self mask")
if o_mask is not None:
self.log("applying other mask")
o_mask[d_data != 0.0] = False
d_data[o_mask] = FILLVALUE
self.log("applying other mask")
d_data[np.isnan(d_data)] = FILLVALUE
if mask_zero_diff:
d_data[np.where(d_data == 0.0)] = FILLVALUE
var = new_net.create_variable(
vname, attrs, s_var.dtype, dimensions=s_var.dimensions
)
var[:] = d_data
self.log("processing variable {0}".format(vname))
def _dt_str(self, dt):
"""for datetime to string for year < 1900"""
dt_str = "{0:04d}-{1:02d}-{2:02d}T{3:02d}:{4:02d}:{5:02}Z".format(
dt.year, dt.month, dt.day, dt.hour, dt.minute, dt.second
)
return dt_str
[docs] def write(self):
"""write the nc object to disk"""
self.log("writing nc file")
assert (
self.nc is not None
), "netcdf.write() error: nc file not initialized"
# write any new attributes that have been set since
# initializing the file
for k, v in self.global_attributes.items():
try:
if self.nc.attributes.get(k) is not None:
self.nc.setncattr(k, v)
except Exception:
self.logger.warn(
"error setting global attribute {0}".format(k)
)
self.nc.sync()
self.nc.close()
self.log("writing nc file")
def _initialize_attributes(self):
"""private method to initial the attributes
of the NetCdf instance
"""
assert (
"nc" not in self.__dict__.keys()
), "NetCdf._initialize_attributes() error: nc attribute already set"
self.nc_epsg_str = "epsg:4326"
self.nc_crs_longname = "http://www.opengis.net/def/crs/EPSG/0/4326"
self.nc_semi_major = float(6378137.0)
self.nc_inverse_flat = float(298.257223563)
self.global_attributes = {}
self.global_attributes["namefile"] = self.model.namefile
self.global_attributes["model_ws"] = self.model.model_ws
self.global_attributes["exe_name"] = self.model.exe_name
self.global_attributes["modflow_version"] = self.model.version
self.global_attributes["create_hostname"] = socket.gethostname()
self.global_attributes["create_platform"] = platform.system()
self.global_attributes["create_directory"] = os.getcwd()
htol, rtol = -999, -999
try:
htol, rtol = self.model.solver_tols()
except Exception as e:
self.logger.warn(
"unable to get solver tolerances:" + "{0}".format(str(e))
)
self.global_attributes["solver_head_tolerance"] = htol
self.global_attributes["solver_flux_tolerance"] = rtol
spatial_attribs = {
"xll": self.model_grid.xoffset,
"yll": self.model_grid.yoffset,
"rotation": self.model_grid.angrot,
"proj4_str": self.model_grid.proj4,
}
for n, v in spatial_attribs.items():
self.global_attributes["flopy_sr_" + n] = v
self.global_attributes[
"start_datetime"
] = self.model_time.start_datetime
self.fillvalue = FILLVALUE
# initialize attributes
self.grid_crs = None
self.zs = None
self.ys = None
self.xs = None
self.nc = None
[docs] def initialize_geometry(self):
"""initialize the geometric information
needed for the netcdf file
"""
try:
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))
if pyproj220:
self.grid_crs = pyproj.CRS(proj4_str)
else:
if proj4_str.lower().startswith("epsg:"):
proj4_str = "+init=" + proj4_str
self.grid_crs = pyproj.Proj(proj4_str, preserve_units=True)
print("initialize_geometry::self.grid_crs = {}".format(self.grid_crs))
vmin, vmax = self.model_grid.botm.min(), self.model_grid.top.max()
if self.z_positive == "down":
vmin, vmax = vmax, vmin
else:
self.zs = self.model_grid.xyzcellcenters[2].copy()
ys = self.model_grid.xyzcellcenters[1].copy()
xs = self.model_grid.xyzcellcenters[0].copy()
# Transform to a known CRS
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_epsg_str = self.nc_epsg_str
if nc_epsg_str.lower().startswith("epsg:"):
nc_epsg_str = "+init=" + nc_epsg_str
nc_crs = pyproj.Proj(nc_epsg_str)
self.transformer = None
print("initialize_geometry::nc_crs = {}".format(nc_crs))
if pyproj220:
print(
"transforming coordinates using = {}".format(self.transformer)
)
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]]
)
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
[docs] def initialize_file(self, time_values=None):
"""
initialize the netcdf instance, including global attributes,
dimensions, and grid information
Parameters
----------
time_values : list of times to use as time dimension
entries. If none, then use the times in
self.model.dis.perlen and self.start_datetime
"""
if self.nc is not None:
raise Exception("nc file already initialized")
if self.grid_crs is None:
self.log("initializing geometry")
self.initialize_geometry()
self.log("initializing geometry")
try:
import netCDF4
except Exception as e:
self.logger.warn("error importing netCDF module")
msg = "NetCdf error importing netCDF4 module:\n" + str(e)
raise Exception(msg)
# open the file for writing
try:
self.nc = netCDF4.Dataset(self.output_filename, "w")
except Exception as e:
msg = "error creating netcdf dataset:\n{}".format(str(e))
raise Exception(msg)
# write some attributes
self.log("setting standard attributes")
self.nc.setncattr(
"Conventions",
"CF-1.6, ACDD-1.3, flopy {}".format(flopy.__version__),
)
self.nc.setncattr(
"date_created", datetime.utcnow().strftime("%Y-%m-%dT%H:%M:00Z")
)
self.nc.setncattr(
"geospatial_vertical_positive", "{}".format(self.z_positive)
)
min_vertical = np.min(self.zs)
max_vertical = np.max(self.zs)
self.nc.setncattr("geospatial_vertical_min", min_vertical)
self.nc.setncattr("geospatial_vertical_max", max_vertical)
self.nc.setncattr("geospatial_vertical_resolution", "variable")
self.nc.setncattr("featureType", "Grid")
for k, v in self.global_attributes.items():
try:
self.nc.setncattr(k, v)
except:
self.logger.warn(
"error setting global attribute {0}".format(k)
)
self.global_attributes = {}
self.log("setting standard attributes")
# spatial dimensions
self.log("creating dimensions")
# time
if time_values is None:
time_values = np.cumsum(self.model_time.perlen)
self.nc.createDimension("time", len(time_values))
for name, length in zip(self.dimension_names, self.shape):
self.nc.createDimension(name, length)
self.log("creating dimensions")
self.log("setting CRS info")
# Metadata variables
crs = self.nc.createVariable("crs", "i4")
crs.long_name = self.nc_crs_longname
crs.epsg_code = self.nc_epsg_str
crs.semi_major_axis = self.nc_semi_major
crs.inverse_flattening = self.nc_inverse_flat
self.log("setting CRS info")
attribs = {
"units": "{} since {}".format(
self.time_units, self.start_datetime
),
"standard_name": "time",
"long_name": NC_LONG_NAMES.get("time", "time"),
"calendar": "gregorian",
"_CoordinateAxisType": "Time",
}
time = self.create_variable(
"time", attribs, precision_str="f8", dimensions=("time",)
)
self.logger.warn("time_values:{0}".format(str(time_values)))
time[:] = np.asarray(time_values)
# Elevation
attribs = {
"units": self.model_grid.units,
"standard_name": "elevation",
"long_name": NC_LONG_NAMES.get("elevation", "elevation"),
"axis": "Z",
"valid_min": min_vertical,
"valid_max": max_vertical,
"positive": self.z_positive,
}
elev = self.create_variable(
"elevation",
attribs,
precision_str="f8",
dimensions=self.dimension_names,
)
elev[:] = self.zs
# Longitude
attribs = {
"units": "degrees_east",
"standard_name": "longitude",
"long_name": NC_LONG_NAMES.get("longitude", "longitude"),
"axis": "X",
"_CoordinateAxisType": "Lon",
}
lon = self.create_variable(
"longitude",
attribs,
precision_str="f8",
dimensions=self.dimension_names[1:],
)
lon[:] = self.xs
self.log("creating longitude var")
# Latitude
self.log("creating latitude var")
attribs = {
"units": "degrees_north",
"standard_name": "latitude",
"long_name": NC_LONG_NAMES.get("latitude", "latitude"),
"axis": "Y",
"_CoordinateAxisType": "Lat",
}
lat = self.create_variable(
"latitude",
attribs,
precision_str="f8",
dimensions=self.dimension_names[1:],
)
lat[:] = self.ys
# x
self.log("creating x var")
attribs = {
"units": self.model_grid.units,
"standard_name": "projection_x_coordinate",
"long_name": NC_LONG_NAMES.get("x", "x coordinate of projection"),
"axis": "X",
}
x = self.create_variable(
"x_proj",
attribs,
precision_str="f8",
dimensions=self.dimension_names[1:],
)
x[:] = self.model_grid.xyzcellcenters[0]
# y
self.log("creating y var")
attribs = {
"units": self.model_grid.units,
"standard_name": "projection_y_coordinate",
"long_name": NC_LONG_NAMES.get("y", "y coordinate of projection"),
"axis": "Y",
}
y = self.create_variable(
"y_proj",
attribs,
precision_str="f8",
dimensions=self.dimension_names[1:],
)
y[:] = self.model_grid.xyzcellcenters[1]
# grid mapping variable
crs = flopy.utils.reference.crs(
prj=self.model_grid.prj, epsg=self.model_grid.epsg
)
attribs = crs.grid_mapping_attribs
if attribs is not None:
self.log("creating grid mapping variable")
self.create_variable(
attribs["grid_mapping_name"], attribs, precision_str="f8"
)
# layer
self.log("creating layer var")
attribs = {
"units": "",
"standard_name": "layer",
"long_name": NC_LONG_NAMES.get("layer", "layer"),
"positive": "down",
"axis": "Z",
}
lay = self.create_variable("layer", attribs, dimensions=("layer",))
lay[:] = np.arange(0, self.shape[0])
self.log("creating layer var")
if self.model_grid.grid_type == "structured":
# delc
attribs = {
"units": self.model_grid.units.strip("s"),
"long_name": NC_LONG_NAMES.get(
"delc", "Model grid cell spacing along a column"
),
}
delc = self.create_variable("delc", attribs, dimensions=("y",))
delc[:] = self.model_grid.delc[::-1]
if self.model_grid.angrot != 0:
delc.comments = (
"This is the row spacing that applied to the UNROTATED grid. "
+ "This grid HAS been rotated before being saved to NetCDF. "
+ "To compute the unrotated grid, use the origin point and this array."
)
# delr
attribs = {
"units": self.model_grid.units.strip("s"),
"long_name": NC_LONG_NAMES.get(
"delr", "Model grid cell spacing along a row"
),
}
delr = self.create_variable("delr", attribs, dimensions=("x",))
delr[:] = self.model_grid.delr[::-1]
if self.model_grid.angrot != 0:
delr.comments = (
"This is the col spacing that applied to the UNROTATED grid. "
+ "This grid HAS been rotated before being saved to NetCDF. "
+ "To compute the unrotated grid, use the origin point and this array."
)
# else:
# vertices
# attribs = {"units": self.model_grid.lenuni.strip('s'),
# "long_name": NC_LONG_NAMES.get("vertices",
# "List of vertices used in the model by cell"),
# }
# vertices = self.create_variable('vertices', attribs, dimensions=('ncpl',))
# vertices[:] = self.model_grid.vertices
# Workaround for CF/CDM.
# http://www.unidata.ucar.edu/software/thredds/current/netcdf-java/
# reference/StandardCoordinateTransforms.html
# "explicit_field"
exp = self.nc.createVariable("VerticalTransform", "S1")
exp.transform_name = "explicit_field"
exp.existingDataField = "elevation"
exp._CoordinateTransformType = "vertical"
exp._CoordinateAxes = "layer"
return
[docs] def initialize_group(
self,
group="timeseries",
dimensions=("time",),
attributes=None,
dimension_data=None,
):
"""
Method to initialize a new group within a netcdf file. This group
can have independent dimensions from the global dimensions
Parameters:
----------
name : str
name of the netcdf group
dimensions : tuple
data dimension names for group
dimension_shape : tuple
tuple of data dimension lengths
attributes : dict
nested dictionary of {dimension : {attributes}} for each netcdf
group dimension
dimension_data : dict
dictionary of {dimension : [data]} for each netcdf group dimension
"""
if attributes is None:
attributes = {}
if dimension_data is None:
dimension_data = {}
if self.nc is None:
self.initialize_file()
if group in self.nc.groups:
raise AttributeError("{} group already initialized".format(group))
self.log("creating netcdf group {}".format(group))
self.nc.createGroup(group)
self.log("{} group created".format(group))
self.log("creating {} group dimensions".format(group))
for dim in dimensions:
if dim == "time":
if "time" not in dimension_data:
time_values = np.cumsum(self.model_time.perlen)
else:
time_values = dimension_data["time"]
self.nc.groups[group].createDimension(dim, len(time_values))
else:
if dim not in dimension_data:
raise AssertionError(
"{} information must be supplied "
"to dimension data".format(dim)
)
else:
self.nc.groups[group].createDimension(
dim, len(dimension_data[dim])
)
self.log("created {} group dimensions".format(group))
dim_names = tuple([i for i in dimensions if i != "time"])
for dim in dimensions:
if dim.lower() == "time":
if "time" not in attributes:
unit_value = "{} since {}".format(
self.time_units, self.start_datetime
)
attribs = {
"units": unit_value,
"standard_name": "time",
"long_name": NC_LONG_NAMES.get("time", "time"),
"calendar": "gregorian",
"Axis": "Y",
"_CoordinateAxisType": "Time",
}
else:
attribs = attributes["time"]
time = self.create_group_variable(
group,
"time",
attribs,
precision_str="f8",
dimensions=("time",),
)
time[:] = np.asarray(time_values)
elif dim.lower() == "zone":
if "zone" not in attributes:
attribs = {
"units": "N/A",
"standard_name": "zone",
"long_name": "zonebudget zone",
"Axis": "X",
"_CoordinateAxisType": "Zone",
}
else:
attribs = attributes["zone"]
zone = self.create_group_variable(
group,
"zone",
attribs,
precision_str="i4",
dimensions=("zone",),
)
zone[:] = np.asarray(dimension_data["zone"])
else:
attribs = attributes[dim]
var = self.create_group_variable(
group,
dim,
attribs,
precision_str="f8",
dimensions=dim_names,
)
var[:] = np.asarray(dimension_data[dim])
[docs] @staticmethod
def normalize_name(name):
return name.replace(".", "_").replace(" ", "_").replace("-", "_")
[docs] def create_group_variable(
self, group, name, attributes, precision_str, dimensions=("time",)
):
"""
Create a new group variable in the netcdf object
Parameters
----------
name : str
the name of the variable
attributes : dict
attributes to add to the new variable
precision_str : str
netcdf-compliant string. e.g. f4
dimensions : tuple
which dimensions the variable applies to
default : ("time","layer","x","y")
group : str
which netcdf group the variable goes in
default : None which creates the variable in root
Returns
-------
nc variable
Raises
------
AssertionError if precision_str not right
AssertionError if variable name already in netcdf object
AssertionError if one of more dimensions do not exist
"""
name = self.normalize_name(name)
if (
name in STANDARD_VARS
and name in self.nc.groups[group].variables.keys()
):
return
if name in self.nc.groups[group].variables.keys():
if self.forgive:
self.logger.warn(
"skipping duplicate {} group variable: {}".format(
group, name
)
)
return
else:
raise Exception(
"duplicate {} group variable name: {}".format(group, name)
)
self.log("creating group {} variable: {}".format(group, name))
if precision_str not in PRECISION_STRS:
raise AssertionError(
"netcdf.create_variable() error: precision "
"string {} not in {}".format(precision_str, PRECISION_STRS)
)
if group not in self.nc.groups:
raise AssertionError(
"netcdf group `{}` must be created before "
"variables can be added to it".format(group)
)
self.var_attr_dict["{}/{}".format(group, name)] = attributes
var = self.nc.groups[group].createVariable(
name,
precision_str,
dimensions,
fill_value=self.fillvalue,
zlib=True,
)
for k, v in attributes.items():
try:
var.setncattr(k, v)
except:
self.logger.warn(
"error setting attribute"
+ "{} for group {} variable {}".format(k, group, name)
)
self.log("creating group {} variable: {}".format(group, name))
return var
[docs] def create_variable(
self,
name,
attributes,
precision_str="f4",
dimensions=("time", "layer"),
group=None,
):
"""
Create a new variable in the netcdf object
Parameters
----------
name : str
the name of the variable
attributes : dict
attributes to add to the new variable
precision_str : str
netcdf-compliant string. e.g. f4
dimensions : tuple
which dimensions the variable applies to
default : ("time","layer","x","y")
group : str
which netcdf group the variable goes in
default : None which creates the variable in root
Returns
-------
nc variable
Raises
------
AssertionError if precision_str not right
AssertionError if variable name already in netcdf object
AssertionError if one of more dimensions do not exist
"""
# Normalize variable name
name = self.normalize_name(name)
# if this is a core var like a dimension...
# long_name = attributes.pop("long_name",name)
if name in STANDARD_VARS and name in self.nc.variables.keys():
return
if (
name not in self.var_attr_dict.keys()
and name in self.nc.variables.keys()
):
if self.forgive:
self.logger.warn(
"skipping duplicate variable: {0}".format(name)
)
return
else:
raise Exception("duplicate variable name: {0}".format(name))
if name in self.nc.variables.keys():
raise Exception("duplicate variable name: {0}".format(name))
self.log("creating variable: " + str(name))
assert (
precision_str in PRECISION_STRS
), "netcdf.create_variable() error: precision string {0} not in {1}".format(
precision_str, PRECISION_STRS
)
if self.nc is None:
self.initialize_file()
# check that the requested dimension exists and
# build up the chuck sizes
# chunks = []
# for dimension in dimensions:
# assert self.nc.dimensions.get(dimension) is not None, \
# "netcdf.create_variable() dimension not found:" + dimension
# chunk = self.chunks[dimension]
# assert chunk is not None, \
# "netcdf.create_variable() chunk size of {0} is None in self.chunks". \
# format(dimension)
# chunks.append(chunk)
self.var_attr_dict[name] = attributes
var = self.nc.createVariable(
name,
precision_str,
dimensions,
fill_value=self.fillvalue,
zlib=True,
) # ,
# chunksizes=tuple(chunks))
for k, v in attributes.items():
try:
var.setncattr(k, v)
except:
self.logger.warn(
"error setting attribute"
+ "{0} for variable {1}".format(k, name)
)
self.log("creating variable: " + str(name))
return var
[docs] def add_global_attributes(self, attr_dict):
"""add global attribute to an initialized file
Parameters
----------
attr_dict : dict(attribute name, attribute value)
Returns
-------
None
Raises
------
Exception of self.nc is None (initialize_file()
has not been called)
"""
if self.nc is None:
# self.initialize_file()
mess = (
"NetCDF.add_global_attributes() should only "
+ "be called after the file has been initialized"
)
self.logger.warn(mess)
raise Exception(mess)
self.log("setting global attributes")
self.nc.setncatts(attr_dict)
self.log("setting global attributes")
def _check_vs_sciencebase(self, md):
"""Check that model bounds read from flopy are consistent with those in ScienceBase."""
xmin, ymin, xmax, ymax = self.bounds
tol = 1e-5
assert md.geospatial_lon_min - xmin < tol
assert md.geospatial_lon_max - xmax < tol
assert md.geospatial_lat_min - ymin < tol
assert md.geospatial_lat_max - ymax < tol
assert md.geospatial_vertical_min - self.vbounds[0] < tol
assert md.geospatial_vertical_max - self.vbounds[1] < tol
[docs] def get_longnames_from_docstrings(self, outfile="longnames.json"):
"""
This is experimental.
Scrape Flopy module docstrings and return docstrings for parameters
included in the list of variables added to NetCdf object. Create
a dictionary of longnames keyed by the NetCdf variable names; make each
longname from the first sentence of the docstring for that parameter.
One major limitation is that variables from mflists often aren't described
in the docstrings.
"""
def startstop(ds):
"""Get just the Parameters section of the docstring."""
start, stop = 0, -1
for i, l in enumerate(ds):
if "Parameters" in l and "----" in ds[i + 1]:
start = i + 2
if l.strip() in ["Attributes", "Methods", "Returns", "Notes"]:
stop = i - 1
break
if i >= start and "----" in l:
stop = i - 2
break
return start, stop
def get_entries(ds):
"""Parse docstring entries into dictionary."""
stuff = {}
k = None
for line in ds:
if (
len(line) >= 5
and line[:4] == " " * 4
and line[4] != " "
and ":" in line
):
k = line.split(":")[0].strip()
stuff[k] = ""
# lines with parameter descriptions
elif k is not None and len(line) > 10: # avoid orphans
stuff[k] += line.strip() + " "
return stuff
# get a list of the flopy classes
# packages = inspect.getmembers(flopy.modflow, inspect.isclass)
packages = [(pp.name[0], pp) for pp in self.model.packagelist]
# get a list of the NetCDF variables
attr = [v.split("_")[-1] for v in self.nc.variables]
# parse docstrings to get long names
longnames = {}
for pkg in packages:
# parse the docstring
obj = pkg[-1]
ds = obj.__doc__.split("\n")
start, stop = startstop(ds)
txt = ds[start:stop]
if stop - start > 0:
params = get_entries(txt)
for k, v in params.items():
if k in attr:
longnames[k] = v.split(". ")[0]
# add in any variables that weren't found
for var in attr:
if var not in longnames.keys():
longnames[var] = ""
with open(outfile, "w") as output:
json.dump(longnames, output, sort_keys=True, indent=2)
return longnames