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

Question

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

Flopy 生成的 shapefile 中的网格原点和旋转看起来没问题。这是用于生成 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)

执行最后一行代码时，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

所以在get_spatialreference中的相关代码位于shapefile_utils.py中的CRS对象中。所以很明显这个 epsg 没有存储在本地 epsg.json 数据库中，所以它正在调用 spatialreference.org。我不知道为什么的全部细节，但看起来您正在寻找的投影在网站上的编码不同：

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

此外，如果您单击它的 esri wkt link，您将获得以下工作 url：

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

如果您想打开软盘代码，也许这应该作为一个问题在软盘代码中的某个时刻解决（尽管考虑到所有 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')

这应该会让你上路。

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

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

flopy