Export raster from binary head file

[Groundwater-modeler]

Hi,

I have a hard time figuring out how to export my heads file as a raster (.tif). I'm thinking i could use the build in Raster-module something like this:

mod_name = "my_model"
grb_name = f"{mod_name}.disv.grb"
m = flopy.mf6.MFSimulation.load(sim_ws=os.getcwd(), sim_name="mfsim.nam", verbosity_level=0)
ml = m.get_model()
modelgrid = ml.modelgrid
headfile = os.path.join( f"{mod_name}.bhd")
timestep1 = 1

# Read the binary head file
headobj = flopy.utils.HeadFile(headfile)

# Extract heads for the specified timesteps
heads = headobj.get_data(totim=timestep1)
    
# Create a PlotMapView for contouring
mapview = flopy.plot.PlotMapView(modelgrid=modelgrid, layer=1)

# Export raster file
raster = flopy.utils.Raster(modelgrid, heads_t1, crs='EPSG:3010',nodataval=-9999,transform=modelgrid)

output_rasterfile_name = f"my_run_heads.tif"
output_raster_path = os.path.join('..', 'output', output_rasterfile_name)

raster.write(output_raster_path)

I'ts just that the "transform" parameter seems to be giving me trouble... perhaps because this is a quad tree grid, but I'm not sure.
Does anyone have ideas of how I could go about? I didn't see a parameter for raster resolution in the module either which i believe must be included somewhere.

Thanks!

David

[Kotarum]

Hi!

I'm not sure if this helps, but when I need to get the head data as a raster, I use this approach:

head = gwf.output.head().get_data()
layer = 0
head_layer = head[layer, :, :]
export_utils.export_array(modelgrid=gwf.modelgrid, a=head_layer, filename='head_layer.tif', nodata=head.max())

This should generate a raster file that QGIS can read.

Hope it helps!

Javi

[Groundwater-modeler]

Thanks for your reply! In this case I'm looking for a method to export heads from a quadtree model. From the looks of it (and from what I can gather in the documentation if I'm not mistaken) this approach applies for structured and uniform models.

Any idea of how this could be handled?

/David

[jlarsen-usgs]

@Groundwater-modeler:

Here are some of the exporting capabilities that'd I'd recommend exploring for writing your heads to either raster or vector data formats.

1. Raster format: FloPy's Raster utility has a raster_from_array() class that can be used to generate a Raster object and to save raster data. This works well for structured grids, as raster data is traditionally structured (rectilinear) and has its geometry defined by an affine transformation. Here's an example of usage

from flopy.utils import Raster

rstr = Raster.from_array(my_3d_head_array, modelgrid=my_structured_grid)
rstr.write("my_head_raster.tif")

2. For unstructured and vertex-based grids I'd recommend creating a shapefile of the head data using the geopandas capabilities on the modelgrid class. Example below:

gdf = modelgrid.geo_dataframe

for ix, layer_arr in enumerate(my_3d_head_array):
    gdf[f"head_{ix}"] = layer_arr.ravel()

gdf.to_file("my_head_shapefile.shp")