Warning1: RingSelf-intersectionatornearpoint112 .4866642030000334 .830899357000078 ERROR1: Cutlinepolygonisinvalid.
使用矢量shp文件剪裁栅格影像时会出现上述问题。提示shp文件自相交。搜集了大量资料,多数回答都是基于PostGIS进行处理。在下载安装过程中很复杂,本人试过多次也没有成功。这里给出一种使用python+ogr函数库的解决方法。
不废话,直接上代码。
from osgeo import ogr, gdal
import os
gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES")
gdal.SetConfigOption("SHAPE_ENCODING", "GBK")
# 思路:https://github.com/whosonfirst-data/whosonfirst-data/issues/1071
# Solve shapefile self-intersection
def s_shp_si(shpFile):
out_shp = shpFile[:-4] + '_ssi' + '.shp'
print(out_shp)
# 打开数据
ds = ogr.Open(shpFile, 0)
if ds is None:
print("打开文件 %s 失败!" % shpFile)
return
print("打开文件%s成功!" % shpFile)
# 获取该数据源中的图层个数,一般shp数据图层只有一个,如果是mdb、dxf等图层就会有多个
m_layer_count = ds.GetLayerCount()
m_layer = ds.GetLayerByIndex(0)
if m_layer is None:
print("获取第%d个图层失败!\n", 0)
return
# 创建输出文件
driver = ogr.GetDriverByName('ESRI Shapefile')
if os.path.exists(out_shp):
driver.DeleteDataSource(out_shp)
outds = driver.CreateDataSource(out_shp)
outlayer = outds.CreateLayer(out_shp[:-4], geom_type=ogr.wkbPolygon, options=["ENCODING=GBK"])
# 获取输出层的要素定义
outLayerDefn = outlayer.GetLayerDefn()
# 对图层进行初始化,如果对图层进行了过滤操作,执行这句后,之前的过滤全部清空
m_layer.ResetReading()
# 获取投影
prosrs = m_layer.GetSpatialRef()
# 添加字段
inLayerDefn = m_layer.GetLayerDefn()
for i in range(0, inLayerDefn.GetFieldCount()):
fieldDefn = inLayerDefn.GetFieldDefn(i)
outlayer.CreateField(fieldDefn)
# loop through the input features
m_feature = m_layer.GetNextFeature()
while m_feature:
# print(m_feature)
o_geometry = m_feature.GetGeometryRef()
# 关键,合并几何
o_geometry = o_geometry.Union(o_geometry)
outfeature = ogr.Feature(outLayerDefn)
outfeature.SetGeometry(o_geometry)
# 遍历每个要素的字段,并设置字段属性
for i in range(0, outLayerDefn.GetFieldCount()):
print(outLayerDefn.GetFieldDefn(i).GetNameRef())
outfeature.SetField(outLayerDefn.GetFieldDefn(i).GetNameRef(), m_feature.GetField(i))
outlayer.CreateFeature(outfeature)
# dereference the features and get the next input feature
outfeature = None
m_feature = m_layer.GetNextFeature()
ds.Destroy()
outds.Destroy()
# 保存投影文件
prjfile = (out_shp[:-4] + '.prj')
fn = open(prjfile, 'w')
fn.write(prosrs.ExportToWkt())
fn.close()
shpFile = r'E:/hhly/hhly.shp' # 裁剪矩形
# # # 注册所有的驱动
ogr.RegisterAll()
s_shp_si(shpFile)
参考链接:
https://blog.csdn.net/gisuuser/article/details/103293038
https://github.com/whosonfirst-data/whosonfirst-data/issues/1071
https://blog.csdn.net/weixin_42924891/article/details/85469881