Skip to content

Commit 8567d8d

Browse files
mjrenomjreno
mjreno
authored and
mjreno
committed
enhanced grid support for netcdf
1 parent 6f3b043 commit 8567d8d

5 files changed

Lines changed: 136 additions & 59 deletions

File tree

‎flopy/discretization/grid.py‎

Lines changed: 7 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1185,12 +1185,14 @@ def read_usgs_model_reference_file(self, reffile="usgs.model.reference"):
11851185
yul = None
11861186
if os.path.exists(reffile):
11871187
with open(reffile) as input:
1188+
print(f"Updating grid based on reference: {reffile}")
11881189
for line in input:
11891190
if len(line) > 1:
11901191
if line.strip()[0] != "#":
11911192
info = line.strip().split("#")[0].split()
11921193
if len(info) > 1:
11931194
data = " ".join(info[1:]).strip("'").strip('"')
1195+
print(f"Grid update on reference: {info}")
11941196
if info[0] == "xll":
11951197
self._xoff = float(data)
11961198
elif info[0] == "yll":
@@ -1201,12 +1203,14 @@ def read_usgs_model_reference_file(self, reffile="usgs.model.reference"):
12011203
yul = float(data)
12021204
elif info[0] == "rotation":
12031205
self._angrot = float(data)
1204-
elif info[0] == "epsg":
1206+
elif info[0] == "length_units":
1207+
self._units = data.lower()
1208+
elif info[0].lower() == "epsg":
12051209
self.epsg = int(data)
12061210
elif info[0] == "proj4":
12071211
self.crs = data
1208-
elif info[0] == "start_date":
1209-
start_datetime = data
1212+
else:
1213+
print(" ->warn: update not applied.")
12101214

12111215
# model must be rotated first, before setting xoff and yoff
12121216
# when xul and yul are provided.

‎flopy/discretization/modeltime.py‎

Lines changed: 37 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,5 @@
11
import calendar
2+
import os
23
from dataclasses import dataclass, field
34
from datetime import datetime, timedelta
45
from difflib import SequenceMatcher
@@ -802,3 +803,39 @@ def reverse(self) -> "ModelTime":
802803
self.start_datetime,
803804
self.steady_state[::-1] if self.steady_state is not None else None,
804805
)
806+
807+
def read_usgs_model_reference_file(self, reffile="usgs.model.reference"):
808+
"""read spatial reference info from the usgs.model.reference file
809+
https://water.usgs.gov/ogw/policy/gw-model/modelers-setup.html"""
810+
if os.path.exists(reffile):
811+
start_date_time = ""
812+
with open(reffile) as input:
813+
print(f"Updating modeltime based on reference: {reffile}")
814+
for line in input:
815+
if len(line) > 1:
816+
if line.strip()[0] != "#":
817+
info = line.strip().split("#")[0].split()
818+
if len(info) > 1:
819+
data = " ".join(info[1:]).strip("'").strip('"')
820+
print(f"ModelTime update on reference: {info}")
821+
if info[0].lower() == "time_units":
822+
self.time_units = data.lower()
823+
elif info[0] == "start_date":
824+
if len(start_date_time) > 0:
825+
start_date_time = f"{data} {start_date_time}"
826+
else:
827+
start_date_time = data
828+
elif info[0] == "start_time":
829+
if len(start_date_time) > 0:
830+
start_date_time = f"{start_date_time} {data}"
831+
else:
832+
start_date_time = data
833+
else:
834+
print(" ->warn: update not applied.")
835+
836+
if len(start_date_time) > 0:
837+
self.start_datetime = start_date_time
838+
839+
return True
840+
else:
841+
return False

‎flopy/discretization/structuredgrid.py‎

Lines changed: 52 additions & 41 deletions
Original file line numberDiff line numberDiff line change
@@ -1770,6 +1770,46 @@ def get_plottable_layer_array(self, a, layer):
17701770
assert plotarray.shape == required_shape, msg
17711771
return plotarray
17721772

1773+
def latlon(self):
1774+
try:
1775+
import warnings
1776+
1777+
from pyproj import Proj
1778+
1779+
epsg = None
1780+
if self.crs is not None:
1781+
epsg = self.crs.to_epsg()
1782+
1783+
proj = Proj(
1784+
f"EPSG:{epsg}",
1785+
)
1786+
1787+
lats = []
1788+
lons = []
1789+
x_local = []
1790+
y_local = []
1791+
for y in self.xycenters[1]:
1792+
for x in self.xycenters[0]:
1793+
x_local.append(x)
1794+
y_local.append(y)
1795+
1796+
x_global, y_global = self.get_coords(x_local, y_local)
1797+
1798+
for i, x in enumerate(x_global):
1799+
lon, lat = proj(x, y_global[i], inverse=True)
1800+
lats.append(lat)
1801+
lons.append(lon)
1802+
1803+
return np.array(lats), np.array(lons)
1804+
1805+
except Exception as e:
1806+
warnings.warn(
1807+
f"Cannot create coordinates from CRS: {e}",
1808+
UserWarning,
1809+
)
1810+
1811+
return None, None
1812+
17731813
def dataset(self, modeltime=None, mesh=None, configuration=None):
17741814
"""
17751815
modeltime : FloPy ModelTime object
@@ -1803,7 +1843,9 @@ def _layered_mesh_dataset(self, ds, modeltime=None, configuration=None):
18031843
}
18041844
ds = ds.assign(var_d)
18051845
ds["time"].attrs["calendar"] = "standard"
1806-
ds["time"].attrs["units"] = f"days since {modeltime.start_datetime}"
1846+
ds["time"].attrs["units"] = (
1847+
f"{modeltime.time_units} since {modeltime.start_datetime}"
1848+
)
18071849
ds["time"].attrs["axis"] = "T"
18081850
ds["time"].attrs["standard_name"] = "time"
18091851
ds["time"].attrs["long_name"] = "time"
@@ -1955,7 +1997,9 @@ def _structured_dataset(self, ds, modeltime=None, configuration=None):
19551997
ds = ds.assign(var_d)
19561998

19571999
ds["time"].attrs["calendar"] = "standard"
1958-
ds["time"].attrs["units"] = f"days since {modeltime.start_datetime}"
2000+
ds["time"].attrs["units"] = (
2001+
f"{modeltime.time_units} since {modeltime.start_datetime}"
2002+
)
19592003
ds["time"].attrs["axis"] = "T"
19602004
ds["time"].attrs["standard_name"] = "time"
19612005
ds["time"].attrs["long_name"] = "time"
@@ -1980,46 +2024,13 @@ def _structured_dataset(self, ds, modeltime=None, configuration=None):
19802024
and configuration["longitude"] is not None
19812025
)
19822026

1983-
if latlon_cfg or self.crs is not None:
1984-
if latlon_cfg:
1985-
lats = configuration["latitude"]
1986-
lons = configuration["longitude"]
1987-
else:
1988-
try:
1989-
import warnings
1990-
1991-
from pyproj import Proj
1992-
1993-
epsg_code = self.crs.to_epsg(min_confidence=90)
1994-
proj = Proj(
1995-
f"EPSG:{epsg_code}",
1996-
)
1997-
1998-
lats = []
1999-
lons = []
2000-
x_local = []
2001-
y_local = []
2002-
for y in self.xycenters[1]:
2003-
for x in self.xycenters[0]:
2004-
x_local.append(x)
2005-
y_local.append(y)
2006-
2007-
x_global, y_global = self.get_coords(x_local, y_local)
2008-
2009-
for i, x in enumerate(x_global):
2010-
lon, lat = proj(x, y_global[i], inverse=True)
2011-
lats.append(lat)
2012-
lons.append(lon)
2013-
2014-
lats = np.array(lats)
2015-
lons = np.array(lons)
2016-
2017-
except Exception as e:
2018-
warnings.warn(
2019-
f"Cannot create coordinates from CRS: {e}",
2020-
UserWarning,
2021-
)
2027+
if latlon_cfg:
2028+
lats = configuration["latitude"]
2029+
lons = configuration["longitude"]
2030+
elif self.crs is not None:
2031+
lats, lons = self.latlon()
20222032

2033+
if lats is not None and lons is not None:
20232034
# create coordinate vars
20242035
var_d = {
20252036
"lat": (["y", "x"], lats.reshape(yc.size, xc.size)),

‎flopy/discretization/vertexgrid.py‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -631,7 +631,9 @@ def dataset(self, modeltime=None, mesh=None, configuration=None):
631631
}
632632
ds = ds.assign(var_d)
633633
ds["time"].attrs["calendar"] = "standard"
634-
ds["time"].attrs["units"] = f"days since {modeltime.start_datetime}"
634+
ds["time"].attrs["units"] = (
635+
f"{modeltime.time_units} since {modeltime.start_datetime}"
636+
)
635637
ds["time"].attrs["axis"] = "T"
636638
ds["time"].attrs["standard_name"] = "time"
637639
ds["time"].attrs["long_name"] = "time"

‎flopy/mf6/mfmodel.py‎

Lines changed: 37 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -458,14 +458,19 @@ def modelgrid(self):
458458
if idomain is None:
459459
force_resync = True
460460
idomain = self._resolve_idomain(idomain, botm)
461+
crs = self._modelgrid.crs
462+
if crs is None and hasattr(dis, "crs"):
463+
crs = dis.crs.get_data()
464+
if crs is not None:
465+
crs = crs[0][1]
461466
self._modelgrid = StructuredGrid(
462467
delc=dis.delc.array,
463468
delr=dis.delr.array,
464469
top=dis.top.array,
465470
botm=botm,
466471
idomain=idomain,
467472
lenuni=dis.length_units.array,
468-
crs=self._modelgrid.crs,
473+
crs=crs,
469474
xoff=self._modelgrid.xoffset,
470475
yoff=self._modelgrid.yoffset,
471476
angrot=self._modelgrid.angrot,
@@ -496,14 +501,19 @@ def modelgrid(self):
496501
if idomain is None:
497502
force_resync = True
498503
idomain = self._resolve_idomain(idomain, botm)
504+
crs = self._modelgrid.crs
505+
if crs is None and hasattr(dis, "crs"):
506+
crs = dis.crs.get_data()
507+
if crs is not None:
508+
crs = crs[0][1]
499509
self._modelgrid = VertexGrid(
500510
vertices=dis.vertices.array,
501511
cell2d=dis.cell2d.array,
502512
top=dis.top.array,
503513
botm=botm,
504514
idomain=idomain,
505515
lenuni=dis.length_units.array,
506-
crs=self._modelgrid.crs,
516+
crs=crs,
507517
xoff=self._modelgrid.xoffset,
508518
yoff=self._modelgrid.yoffset,
509519
angrot=self._modelgrid.angrot,
@@ -525,6 +535,11 @@ def modelgrid(self):
525535
idomain = dis.idomain.array
526536
if idomain is None:
527537
idomain = np.ones(dis.nodes.array, dtype=int)
538+
crs = self._modelgrid.crs
539+
if crs is None and hasattr(dis, "crs"):
540+
crs = dis.crs.get_data()
541+
if crs is not None:
542+
crs = crs[0][1]
528543
if cell2d is None:
529544
if (
530545
self.simulation.simulation_data.verbosity_level.value
@@ -557,7 +572,7 @@ def modelgrid(self):
557572
idomain=idomain,
558573
lenuni=dis.length_units.array,
559574
ncpl=ncpl,
560-
crs=self._modelgrid.crs,
575+
crs=crs,
561576
xoff=self._modelgrid.xoffset,
562577
yoff=self._modelgrid.yoffset,
563578
angrot=self._modelgrid.angrot,
@@ -590,14 +605,19 @@ def modelgrid(self):
590605
if idomain is None:
591606
force_resync = True
592607
idomain = self._resolve_idomain(idomain, botm)
608+
crs = self._modelgrid.crs
609+
if crs is None and hasattr(dis, "crs"):
610+
crs = dis.crs.get_data()
611+
if crs is not None:
612+
crs = crs[0][1]
593613
self._modelgrid = VertexGrid(
594614
vertices=dis.vertices.array,
595615
cell1d=dis.cell1d.array,
596616
top=None,
597617
botm=botm,
598618
idomain=idomain,
599619
lenuni=dis.length_units.array,
600-
crs=self._modelgrid.crs,
620+
crs=crs,
601621
xoff=self._modelgrid.xoffset,
602622
yoff=self._modelgrid.yoffset,
603623
angrot=self._modelgrid.angrot,
@@ -628,14 +648,19 @@ def modelgrid(self):
628648
if idomain is None:
629649
force_resync = True
630650
idomain = self._resolve_idomain(idomain, botm)
651+
crs = self._modelgrid.crs
652+
if crs is None and hasattr(dis, "crs"):
653+
crs = dis.crs.get_data()
654+
if crs is not None:
655+
crs = crs[0][1]
631656
self._modelgrid = StructuredGrid(
632657
delc=dis.delc.array,
633658
delr=dis.delr.array,
634659
top=None,
635660
botm=botm,
636661
idomain=idomain,
637662
lenuni=dis.length_units.array,
638-
crs=self._modelgrid.crs,
663+
crs=crs,
639664
xoff=self._modelgrid.xoffset,
640665
yoff=self._modelgrid.yoffset,
641666
angrot=self._modelgrid.angrot,
@@ -666,14 +691,19 @@ def modelgrid(self):
666691
if idomain is None:
667692
force_resync = True
668693
idomain = self._resolve_idomain(idomain, botm)
694+
crs = self._modelgrid.crs
695+
if crs is None and hasattr(dis, "crs"):
696+
crs = dis.crs.get_data()
697+
if crs is not None:
698+
crs = crs[0][1]
669699
self._modelgrid = VertexGrid(
670700
vertices=dis.vertices.array,
671701
cell2d=dis.cell2d.array,
672702
top=None,
673703
botm=botm,
674704
idomain=idomain,
675705
lenuni=dis.length_units.array,
676-
crs=self._modelgrid.crs,
706+
crs=crs,
677707
xoff=self._modelgrid.xoffset,
678708
yoff=self._modelgrid.yoffset,
679709
angrot=self._modelgrid.angrot,
@@ -1362,13 +1392,6 @@ def write(
13621392
if write_netcdf:
13631393
# set data storage to write ascii for netcdf
13641394
pp._set_netcdf_storage()
1365-
if (
1366-
pp.package_type.startswith("dis")
1367-
and hasattr(pp, "crs")
1368-
):
1369-
crs = pp.crs.get_data()
1370-
if crs is not None and self.modelgrid.crs is None:
1371-
self.modelgrid.crs = crs[0][1]
13721395
if (
13731396
self.simulation_data.verbosity_level.value
13741397
>= VerbosityLevel.normal.value
@@ -2379,7 +2402,7 @@ def _write_netcdf(self, mesh=None):
23792402
nc_meta["attrs"]["history"] = f"first created {timestamp}"
23802403
if mesh is None:
23812404
nc_meta["attrs"]["Conventions"] = "CF-1.11"
2382-
elif mesh.upper() is "LAYERED":
2405+
elif mesh.upper() == "LAYERED":
23832406
nc_meta["attrs"]["Conventions"] = "CF-1.11 UGRID-1.0"
23842407

23852408
ds = self.update_dataset(

0 commit comments

Comments
 (0)