Flopy-generated shapefile not accounting for MODFLOW DIS lenuni value?

Question

I've used Flopy to generate a shapefile of polygon features representing MODFLOW River Package features. However, the size of the grid cell polygon features in the shapefile are 3.28 times larger than they should be. The length units of my model are in feet (the LENUNI variable in the MODFLOW Discretization Package of my model is equal to 1), and I'm using NAD83 / UTM Zone 16N (length unit is meters, EPSG:26916). Therefore, it looks like the conversion between MODFLOW model units (in feet) and GIS coordinate reference system (in meters) isn't happening for some reason.

The grid origin and rotation in the Flopy-generated shapefile look okay. Here is the Flopy code used to generate the shapefile:

model_ws = os.getcwd()
m = flopy.modflow.Modflow.load("model.nam", model_ws=model_ws, verbose=False,
                               check=False, exe_name="MODFLOW-NWT_64.exe")

grid = m.modelgrid
delr = grid.delr
delc = grid.delc
xll =  660768.2212
yll = 3282397.889
rot = -16.92485016
model_epsg = 26916
m.modelgrid.set_coord_info(xoff=xll, yoff=yll, angrot=rot, epsg='EPSG:26916')

m.riv.stress_period_data.export('{0}/riv_features.shp'.format(model_ws), sparse=True)

When the last line of code is executed, the shapefile is written to disk, but the following error messages precede the message confirming that the shapefile was output:

(<class 'urllib.error.HTTPError'>, <HTTPError 404: 'NOT FOUND'>, <traceback object at 0x11646208>)
No internet connection or epsg code EPSG:26916 not found at https://spatialreference.org/ref/epsg/EPSG:26916/esriwkt
No internet connection or epsg code EPSG:26916 not found at https://spatialreference.org/ref/esri/EPSG:26916/esriwkt

The URL associated with the above error messages is https://spatialreference.org/ref/epsg/EPSG:26916/esriwkt. This URL displayed the following:

Not found, /ref/epsg/EPSG:26916/esriwkt.

So could the problem be that Flopy is not getting the information that it needs from spatialreference.org? If so, is the URL that Flopy is generating incorrect? Is there something in my code that is incorrect?

Thanks very much.

Answer 1

So the relevant code in in get_spatialreference in the CRS object located in shapefile_utils.py. So apparently this epsg is not stored in the local epsg.json database, so it is calling spatialreference.org. I don't know the full details why, but it looks as though the projection you are looking for is coded differently on the website:

https://spatialreference.org/ref/epsg/nad83-utm-zone-16n/

Furthermore, if you click on its esri wkt link you get the following working url:

https://spatialreference.org/ref/epsg/nad83-utm-zone-16n/esriwkt/

Maybe this should be addressed at some point in the flopy code as an issue if you want to open one up, (although accounting for all CRS edge cases is a daunting task). However, in the mean time I think you can get around this by assigning epsg the following value:

m.modelgrid.set_coord_info(xoff=xll, yoff=yll, angrot=rot, epsg='nad83-utm-zone-16n')

Or if that's not ideal you should be able add it to the export function:

.riv.stress_period_data.export('{0}/riv_features.shp'.format(model_ws), sparse=True, epsg='nad83-utm-zone-16n')

That way the proper url will be called when getting the CRS info. Let me know if you have any issues.

EDIT:

I will look in to if flopy can handle the unit length conversions. In the mean time, since you are a python user you can do the following to scale your vectors if you install geopandas:

import geopandas as gpd
from shapely.ops import unary_union
df = gpd.read_file(proj_path + '/polytest.shp') #replace with your file location
union = unary_union(df.geometry.values)
df.geometry = df.geometry.scale(xfact=.304878,yfact=.304878, origin=union.centroid)  #You may require another origin possibly
df.to_file(proj_path + '/new_shapes.shape')

That should get you on your way.