Flopy生成的形状文件没有考虑MODFLOW DIS lenuni值

Question

我使用Flopy生成了一个多边形特征的形状文件，表示MODFLOW River Package特征。然而，shapefile中网格单元多边形特征的大小是其应有大小的3.28倍。我的模型的长度单位为英尺(我的模型MODFLOW离散化包中的LENNUI变量等于1(，我使用的是NAD83/UTM Zone 16N(长度单位为米，EPSG:26916(。因此，由于某种原因，MODFLOW模型单位(以英尺为单位(和GIS坐标参考系统(以米为单位(之间的转换似乎没有发生。

Flopy生成的形状文件中的栅格原点和旋转看起来不错。以下是用于生成形状文件的Flopy代码：

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)

当执行最后一行代码时，shapefile被写入磁盘，但在确认shapefile已输出的消息之前会出现以下错误消息：

(<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

与上述错误消息关联的URL为https://spatialreference.org/ref/epsg/EPSG:26916/esriwkt.此URL显示以下内容：

未找到，/ref/epsg/epsg:26916/esriwkt。

那么问题可能是Flopy没有从spatialreference.org获得所需的信息吗？如果是，Flopy生成的URL是否不正确？我的代码中有什么不正确的地方吗？

非常感谢。

Answer 1

因此，位于shapefile_utils.py中CRS对象的get_spatialreference中的相关代码。因此，显然这个epsg没有存储在本地epsg.json数据库中，因此它正在调用spatialreference.org。我不知道原因的全部细节，但看起来你要查找的投影在网站上的编码不同：

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

此外，如果你点击它的esri wkt链接，你会得到以下工作网址：

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

如果你想打开一个问题，也许这应该在flopy代码的某个时候作为一个问题来解决(尽管考虑所有CRS边缘情况是一项艰巨的任务(。然而，与此同时，我认为你可以通过给epsg分配以下值来解决这个问题：

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

或者，如果这不理想，你应该可以将其添加到导出功能中：

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

这样，当获取CRS信息时，将调用正确的url。如果你有任何问题，请告诉我。

编辑：

我会看看flopy是否能处理单位长度的转换。同时，由于您是python用户，如果您安装了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')

这应该会让你走上正轨。