Issues with grid refinement "level 1" only with some .shp files

[Larkiie]

Hi all,
I am quite new to both python and flopy/modflow, so don't hesitate to ask if the issue isn't clear.

As part of my masters thesis I'm working with my supervisor to set up modflow simulations with flopy.
I've run into a peculiar issue, that we couldn't find a reason for, so my supervisor suggested asking on github.

The thing is - one of the elements I'm varying as part of my project is grid refinement. It is set up with refining with "levels", meaning level = 1 means that the cells in the specified area gets halved, creating a finer grid (i.e. in my case from 100x100m to 50x50). If you use a higher "level" is refines to 25 etc.

I do it with specifying an area with a shapefile created in qgis. However, I've encountered an issue where if I refine on level 1, it creates the grid fine, but when I run the simulation it fails to truncate.
It is only when running level 1 - I can use the same shapefile to run level 0, 2, 3, 4 without issue.
Some shapefiles I've manage to succeed with - but I cannot see why they would be different than other random shapefiles. 1 shapefile is a multilinestring consisting of 25 parts, 1 is linestring consisting of 1 part. I did make it work with a polygon as well. But other shp files it just creates the grid fine, but the simulation fails to truncate - it hasn't mattered if they were a multilinestring of few or many parts. I do not know why some works.

It is so strange. I've tested creating different shapefiles. I've checked if it's a naming convention - no. It doesn't seem to be an issue with the setup or my pc's limitations, since it works on all the other levels. It has the correct number of cells.

Has anybody encountered something similar, or have an Idea what could cause it?

I'll insert part of the code below:

%%time

Specify refinement level

level = 1 # (int) – The level of refinement for this features

Specify grid resolution etc

xoff = 528500.0 # bottom left # xoff = 530800.0
yoff = 6251600.0 # bottom left # yoff = 6251900.0
Lx = 541500-xoff # 5900
Ly = 6260300-yoff # 5900
extent = (0, Lx, 0, Ly)
vmin, vmax = 0.0, 100.0

specify a polygon of the refinement area for the quadtree grid

shapefile_path = '../../source_data/shapefiles/extents/ProjectArea.shp'

===========================

MF6 simulation to make grid

===========================

sim = flopy.mf6.MFSimulation()
gwf = flopy.mf6.ModflowGwf(sim)
dx = dy = 100.0
nr = int(Ly / dy)
nc = int(Lx / dx)
dis6 = flopy.mf6.ModflowGwfdis(
gwf,
nrow=nr, ncol=nc,
delr=dy, delc=dx,
xorigin=xoff, yorigin=yoff,
length_units='meters'
)

g = Gridgen(gwf.modelgrid, model_ws=temp_path)

g.add_refinement_features(
features=shapefile_path,
featuretype='polygon',
level=level,
layers=range(1))

g.build(verbose=False)
gridprops_vg = g.get_gridprops_vertexgrid()
quadtree_grid = flopy.discretization.VertexGrid(**gridprops_vg)