Exporting to netCDF

[ed76]

Hi all,

I want to export heads and specific discharge from a modflow-nwt-model (created with Groundwater Vistas) to netCDF. To keep it simple, I limited myself at first to heads and wrote a small py-script, closely following the flopy-tutorial (‘Exporting to netCDF’, ‘Freyberg model’). The netCDF file is created, but the content (viewed in QGIS) is not as expected:

1. I would have expected only heads to be exported. In addition to the heads, the nc file also contains: Elevation (all layers), Latitude, Longitude, x/y coordinate of projection.

2. All values are displayed incorrectly with incorrect coordinates. xoff=479750.0, yoff=5758550.0 was set but the coordinates in qgis are 0,0. furthermore the model is in reality > 6000 m in x- and y-direction. the content of the netcdf file however is 205m in x-direction and 225m in y-direction. this is exactly the number of columns/rows! the model has uneven cell sizes which leads to a distortion of the image.

in summary, each cell is exported with the dimensions 1x1m although in reality the cells are much larger and have unequal dimensions. xoff and yoff are not taken into account.i would be very happy if someone could help me. below you find first the input script and second the output.

py script:

# # Exporting to netCDF
 
# +
import os
import sys
from tempfile import TemporaryDirectory
 
import flopy
 
print(sys.version)
print(f"flopy version: {flopy.__version__}") # -
 
# Load model
 
nam_file = "gv66.nam"
 
model_ws = os.getcwd()
 
ml = flopy.modflow.Modflow.load(nam_file, model_ws=model_ws, check=False)
 
# We can see the ``Modelgrid`` instance has generic entries, as does ``start_datetime``
 
ml.modelgrid
 
ml.modeltime.start_datetime
 
# Setting the attributes of the ``ml.modelgrid`` is easy:
 
ml.modelgrid.set_coord_info(
    xoff=479750.0, yoff=5758550.0, angrot=0.0, crs=25832
)
ml.dis.start_datetime = "7/4/1776"
 
ml.modeltime.start_datetime
 
# load binary head
fhead = os.path.join(model_ws, "gv66.hds")
 
hds = flopy.utils.HeadFile(fhead)
 
export_dict = {"hds": hds}
 
# Define the path where the NetCDF file should be saved
pth = model_ws  # or specify another directory path
 
# export head output to netcdf
fnc = flopy.export.utils.output_helper(
    os.path.join(pth, "output.nc"), ml, export_dict
)
# -
 
try:
    # ignore PermissionError on Windows
    temp_dir.cleanup()
except:
    pass

 
Output:

3.12.7 (tags/v3.12.7:0b05ead, Oct  1 2024, 03:06:41) [MSC v.1941 64 bit (AMD64)]
flopy version: 3.8.2
initialize_geometry::
model crs: EPSG:25832
initialize_geometry::nc_crs = EPSG:4326
transforming coordinates using = proj=pipeline step inv proj=utm zone=32 ellps=GRS80 step proj=unitconvert xy_in=rad xy_out=deg
 
Warning (from warnings module):
  File "C:\Users\xxx\AppData\Roaming\Python\Python312\site-packages\flopy\export\netcdf.py", line 755
    "date_created", datetime.utcnow().strftime("%Y-%m-%dT%H:%M:%SZ")
DeprecationWarning: datetime.datetime.utcnow() is deprecated and scheduled for removal in a future version. Use timezone-aware objects to represent datetimes in UTC: datetime.datetime.now(datetime.UTC).

[ed76]

i have spent a few more hours working on the problem described. i have tried other options from the tutorial (e.g. “# export the whole model (inputs and outputs)” etc.). i always succeed in creating a netcf file, but the problem with the cell dimensions and coordinates described above still persists. i therefore think that there is either an error/bug in flopy or in the tutorial. i would be very happy if someone could help, especially as the problem is easy to reproduce.

[langevin-usgs]

Hi @ed76, sorry to hear about the troubles with the netcdf export. The flopy capability that you're using probably doesn't get a lot of use, so it's not too surprising that you are running into issues. The one thing that might be causing the problem is the qgis support for the netcdf file that flopy is writing. I'm just speculating here, but that older netcdf support in flopy was primarily focused on geographic coordinates, so I'm thinking that the lat/lons are correct, but the projected coordinates are not fully supported (seems like delr/delc are not being used to construct projected coordinates).

We are working on more comprehensive support for netcdf in MODFLOW 6, including both structured and unstructured grids. The netcdf files can be written directly by MODFLOW 6 when it runs. Our testing so far shows that these netcdf files can be viewed in qgis. I know this probably doesn't help you at the moment, but it may be something to keep in mind for later.

Maybe @mjreno has something to add here.

[mjreno]

@langevin-usgs @ed76 apologies for not responding earlier. I was thinking about running through this exercise myself to see if there was anything I can add. @ed76 I'll plan to do this early this week and report back if that will still work for you.

[mjreno]

Hi @ed76 would you mind adding your generated NetCDF file to this discussion? I have generated a couple of exports with local tests and would like to compare to your file if possible.

[ed76]

thanks @langevin-usgs and @mjreno for your support.

attached is the netcdf file that was created with the script above.

output.zip

[mjreno]

Thanks for adding the export file @ed76. Just adding a few observations here- feel free to clarify if I'm misunderstanding or get something wrong.

In general, I believe the extra variables in the file are to assist with visualization in tools like QGIS. The head array itself is a timeseries data array and should share the same dimensions of the model grid (along with the time dimension). Here is the definition from the provided file:

   float head(time, layer, y, x) ;
                head:_FillValue = -99999.9f ;
                head:long_name = "head" ;
                head:coordinates = "time layer latitude longitude" ;
                head:min = 78.15503f ;
                head:max = 126.3569f ;
                head:units = "meters" ;

Note the "coordinates" attribute which references other variables in the file. These coordinate variables are used by visualization tools to associate the data array to a location so it can be correctly displayed.

It looks to me like "xoff" and "yoff" are correctly taken into account in the file- the "x_proj" and y_proj" arrays use these values and grid dimensions provided from delr/delc to write the projected coordinate values for the given CRS. These values are then used to generate latitude and longitude arrays with a pyproj transformer specific to the CRS.

Note the EPSG number in the file is 4326 but above in the example script it is 25832- not sure which is correct. I'm not seeing the grid display correctly in my QGIS either but I am not able to find either of these EPSG's in my QGIS version and I think that would be the first step.

[ed76]

the excerpt of the definitions you posted seems correct (head min/max, no-data values etc.) but they are not displayed correctly in qgis. i find your comment interesting

“Note the EPSG number in the file is 4326 but above in the example script it is 25832- not sure which is correct.” .

when I look at the output (first post above):

“initialize_geometry::
model crs: EPSG:25832
initialize_geometry::nc_crs = EPSG:4326
transforming coordinates using = proj=pipeline step inv proj=utm zone=32 ellps=GRS80 step proj=unitconvert xy_in=rad xy_out=deg”

this is exactly the point i don't understand. the xoff and yoff values refer to the coordinate system epsg 25832 (utm, i live in germany). this is definitely correct. i don't understand why everything is transformed to epsg 4326, maybe that's the problem?

if i understand you correctly, you can't set epsg 25832 in qgis? you could search for “25832”, then you should find it, see attached photo.

thanks for your efforts!

[mjreno]

In my reading of the code, 4326 is the hard coded CRS of the exported file. In this case, 25832 is provided as the "model CRS" and when any such model CRS is provided a transform is applied to convert to 4326.

What do you see if you don't provide a model CRS for the export and then just try to explicitly assign the CRS in QGIS?

[ed76]

"What do you see if you don't provide a model CRS for the export and then just try to explicitly assign the CRS in QGIS?"

in the end nothing changes, everything stays the same, attached is the newly generated nc-file without indication of crs

output.zip