0

我使用 Flopy 生成了代表 MODFLOW River Package 特征的多边形特征的 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 是否不正确？我的代码中有什么不正确的吗？

非常感谢。

4

1 回答 1

0

所以相关代码在位于 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/

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

那应该让你上路。

于 2020-04-08T16:54:43.900 回答