Flopy生成的shapefile不考虑MODFLOW DIS lenuni值？

Question

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

Flopy生成的shapefile中的网格原点和旋转看起来正常。下面是用于生成shapefile的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链接，你会得到以下有效的url：

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')

这应该能让你上路了。