前言
由于工程上需要进行各种矢量的逻辑运算，类似arcgis的merge，union等，在开源环境下目前还没有人做过特别好的总结，这里由于涉及的逻辑运算处理很多，所以写一篇文章总结一下。（封面图:wallhaven，侵删）
多部分到单部分
这个挺有意思，英文是‘multipart to singleparts’，这个出现的原因是，由于使用分割（segmentation）以后做polygonized会导致‘shapefile polygon to polyline’锯齿化严重，，在做Multipolyline简化的时候，由于shapefile特征之间的离散导致简化会各个特征呈现不同状态，导致简化后反而更复杂。
那这个解决办法首先是做一次级联（Cascade），但是基于ogr的级联是通过迭代完成的，python在多节点wkt的polyline中循环显得十分缓慢，最后导致循环失败，不得已得找找别的方案。持能不造轮就不造的想法，首先当然是拥向矢量处理的业界大拿——geopandas。二在geopandas中，也轻而易举找到了相关解决方案，如上文所述，使用‘multipart to singleparts’以后，多feature会转为单feature，也就是说他们联合了，变成了单属性，然后通过生成矢量范围的Extend，得到区域的polygon，最后通过求交，得到wkt不重复的multipolyline。技术路线如下：
由于栅格无法直接转线，所以这里前两个步骤是一种折中方案，第三步求研究区范围和矢量线的交集，实际就是求一个整体的feature，这也是一种变换的思路。
单特征的polyline简化后，不会出现杂乱的情况。这里由于简化容差（tolerance）设置的数值比较大，所以造成很多转折点都被简化了。实际情况的容差值应该根据其分辨率数值进行设置和调整。
说了基本原理，其实应该有了大概的一个思路。那实现的代码是通过geopandas的simplify函数直接实现，前处理函数这里也一起给出。
1
2
3
4
5
6
7
8
9
10
11
|
def merge_all_feature_in_one(input,output):
print("开始合并feature")
gdf = gpd.read_file(input)
geom = gdf['geometry']
new_geom = gpd.tools.collect(geom)
df = {'id': [0], 'geometry': [new_geom]}
new_gdf = gpd.GeoDataFrame(df, crs="EPSG:4326")
new_gdf.to_file(output)
print("合并完成")
|
这个主要也是通过geopandas的collect函数实现。
这个是简化函数，其实也能用作面简化
1
2
3
4
|
def simplifyshp(input, output, tolerance):
gdf_main = gp.read_file(input)
simp = gdf_main.simplify(tolerance=tolerance , preserve_topology=True)
simp.to_file(output)
|
影像有效范围提取
经过正射校正的影像，往往在投影空间里是有一定旋转角度的，这点在ENVI中能清晰看到，有些实验需要提取影像有效范围，一般思路应该是通过gdal的六参数（Geotransform）取到四个角点，但实际上这样准确吗？或许在没经过正射校正之前是准确的，退一步说其实正射校正是最后的处理步骤，在该步骤之前进行有效区域提取，理论上确实没问题。但是工程上有些地方需要先进行正射校正，这里我就以正射校正为先的思路继续讨论。
那进行了正射校正以后，六参数不是对应的有效区域的角点，该怎么办呢，我这里就设计了一个二值提取思路，首先将影像二值化，以1为分段点，1一下为0，以上为1，二值化以后进行polygonized，然后根据polygon的属性值选取，就能得到影像有效区域。
那么先给出二值化函数：
1
2
3
4
5
6
7
8
9
10
|
def raster_binary(input_raster,out_raster):
dataset = gdal.Open(input_raster)
im_width = dataset.RasterXSize
im_height = dataset.RasterYSize
im_geotrans = dataset.GetGeoTransform()
im_proj = dataset.GetProjection()
im_data = dataset.ReadAsArray(0, 0, im_width, im_height)
ret, border0 = cv2.threshold(im_data, 0, 1, cv2.THRESH_BINARY)
write_img(out_raster, im_proj, im_geotrans, border0)
del dataset
|
其中辅助函数write_img如下：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
|
def write_img(filename, im_proj, im_geotrans, im_data):
if 'int8' in im_data.dtype.name:
datatype = gdal.GDT_Byte
elif 'int16' in im_data.dtype.name:
datatype = gdal.GDT_UInt16
else:
datatype = gdal.GDT_Float32
if len(im_data.shape) == 3:
im_bands, im_height, im_width = im_data.shape
else:
im_bands, (im_height, im_width) = 1, im_data.shape
driver = gdal.GetDriverByName("GTiff")
dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)
dataset.SetGeoTransform(im_geotrans)
dataset.SetProjection(im_proj)
if im_bands == 1:
dataset.GetRasterBand(1).WriteArray(im_data)
else:
for i in range(im_bands):
dataset.GetRasterBand(i + 1).WriteArray(im_data[i])
del dataset
|
通过二值化影像，取栅格转矢量polygon：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
|
def PolygonizeTheRaster_bina(inputfile,outputfile):
dataset = gdal.Open(inputfile, gdal.GA_ReadOnly)
srcband=dataset.GetRasterBand(1)
im_proj = dataset.GetProjection()
prj = osr.SpatialReference()
prj.ImportFromWkt(im_proj)
drv = ogr.GetDriverByName('ESRI Shapefile')
dst_ds = drv.CreateDataSource(outputfile)
dst_layername = 'out'
dst_layer = dst_ds.CreateLayer(dst_layername, srs=prj)
dst_fieldname = 'DN'
fd = ogr.FieldDefn(dst_fieldname, ogr.OFTInteger)
dst_layer.CreateField(fd)
dst_field = 0
gdal.Polygonize(srcband, None, dst_layer, dst_field)
|
通过选取shapefile polygon的’DN’值字段为1的feature，创建新shapefile，完成提取有效区域：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
|
def SelectByAttribute(InShp, outShp):
ds = ogr.Open(InShp,0)
layer = ds.GetLayer(0)
driver = ogr.GetDriverByName("ESRI Shapefile")
ds_new = driver.CreateDataSource(outShp)
spatialref_new = osr.SpatialReference()
spatialref_new.ImportFromEPSG(4326)
geomtype = ogr.wkbPolygon
layer_new = ds_new.CreateLayer(outShp, srs=spatialref_new, geom_type=geomtype,)
layer_new.CreateField(ogr.FieldDefn('DN', ogr.OFTInteger))
attriNames = ogr.OFTInteger
features = [1]
for i in range(layer.GetFeatureCount()):
feature = layer.GetFeature(i)
geom_polygon = feature.GetGeometryRef()
name = feature.GetField('DN')
feat = ogr.Feature(layer_new.GetLayerDefn())
if name in attriNames:
feat.SetGeometry(geom_polygon)
feat.SetField('DN', name)
features.append(feat)
for f in features:
layer_new.CreateFeature(f)
ds.Destroy()
ds_new.Destroy()
|
求取两幅影像的共同区域
这个就更简单了，其实放在矢量逻辑里就是求交，通过geopandas或者其他shp处理包都能轻易实现。但是geopandas是基于其他shp处理包的汇聚，理论上应该处理性能更好。
1
2
3
4
5
|
def get_intersect_shp(out_shp_ring1 , out_shp_ring2, outshp_intersect):
gdf_left = gpd.read_file(out_shp_ring1)
gdf_right = gpd.read_file(out_shp_ring2)
gdf_intersect = gpd.overlay(gdf_left,gdf_right,'intersection')
gdf_intersect.to_file(outshp_intersect)
|
这里的两个输入，应该根据文章上节求影像有效范围后，作为输入的两个shp。
其他函数
除了以上几个主要函数以外，还做了一些很鸡肋的函数，不过在这其中有些可能会对大家有启发，所以还是逐个介绍一下。
裁剪就不用说了，基于gdal就能轻松实现：
1
2
3
4
5
6
7
8
9
10
|
def clip_raster_from_intersect(input_main_raster, input_shp, output_raster):
input_raster=gdal.Open(input_main_raster)
ds = gdal.Warp(output_raster,
input_raster,
format = 'GTiff',
cutlineDSName = input_shp,
cutlineWhere="FIELD = 'whatever'",
dstNodata = -9999)
ds=None
|
如果想拿到影像有效范围的四个角点，该怎么做呢。实际上通过上文求出有效区域的shp以后，利用shp的extend就可以轻松拿到角点坐标，这里给出角点函数
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
|
def get_extends(intersect_shp, out_shp_point):
inDriver = ogr.GetDriverByName("ESRI Shapefile")
inDataSource = inDriver.Open(intersect_shp, 0)
inLayer = inDataSource.GetLayer()
corners = inLayer.GetExtent()
point1 = geometry.Point(corners[0,0],corners[1,0])
point2 = geometry.Point(corners[0,1],corners[1,1])
point3 = geometry.Point(corners[0,3],corners[1,3])
point4 = geometry.Point(corners[0,2],corners[1,2])
pointshp = gpd.GeoSeries([point1, point2, point3, point4],
crs='EPSG:4326', # 指定坐标系为WGS 1984
index=['0', '1'] # 相关的索引
)
pointshp.to_file(out_shp_point,driver='ESRI Shapefile',encoding='utf-8')
inDataSource = None
out_shp_point = None
|
这里其实可以写个字典，或者list迭代导入geometry，但是范围角点就四个，懒得写循环了🤣
其实还有个ogr版本的简化，或者说平滑，原理是通过缓冲区构建缓冲线，完成的平滑曲线，但是这个在多特征中行不通，也就是会出现简化后的线叠合，在前文中有提到，这里也给出来做个参考，这个函数是别人写的，我就不做讨论了。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
|
def smoothing(inShp, fname, bdistance=2):
"""
:param inShp: 输入的矢量路径
:param fname: 输出的矢量路径
:param bdistance: 缓冲区距离
:return:
"""
ogr.UseExceptions()
in_ds = ogr.Open(inShp)
in_lyr = in_ds.GetLayer()
# 创建输出Buffer文件
driver = ogr.GetDriverByName('ESRI Shapefile')
if Path(fname).exists():
driver.DeleteDataSource(fname)
# 新建DataSource，Layer
out_ds = driver.CreateDataSource(fname)
out_lyr = out_ds.CreateLayer(fname, in_lyr.GetSpatialRef(), ogr.wkbMultiLineString)
def_feature = out_lyr.GetLayerDefn()
# 遍历原始的Shapefile文件给每个Geometry做Buffer操作
for feature in tqdm.tqdm(in_lyr):
geometry = feature.GetGeometryRef()
buffer = geometry.Buffer(bdistance,1).Buffer(-bdistance,1)
out_feature = ogr.Feature(def_feature)
out_feature.SetGeometry(buffer)
out_lyr.CreateFeature(out_feature)
out_feature = None
out_ds.FlushCache()
del in_ds, out_ds
|
还有个逻辑求和，不过这个基本没啥用，对多特征且特征之间严格接触的矢量文件，使用它会导致全部特征融成了一个大的shp块，但是其中的wkt点全部都没了，是真正意义的融合，具体的实验shp文件截图我这忘记截了，有兴趣的读者可以自己实验看下。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
|
def uni(shpPath, fname):
"""
:param shpPath: 输入的矢量路径
:param fname: 输出的矢量路径
:return:
"""
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(shpPath, 1)
layer = dataSource.GetLayer()
# 新建DataSource，Layer
out_ds = driver.CreateDataSource(fname)
out_lyr = out_ds.CreateLayer(fname, layer.GetSpatialRef(), ogr.wkbLineStringZM)
def_feature = out_lyr.GetLayerDefn()
# 遍历原始的Shapefile文件给每个Geometry做Buffer操作
current_union = layer[0].Clone()
print('the length of layer:', len(layer))
if len(layer) == 0:
return
for i, feature in tqdm.tqdm(enumerate(layer)):
geometry = feature.GetGeometryRef()
if i == 0:
current_union = geometry.Clone()
current_union = current_union.Union(geometry).Clone()
if i == len(layer) - 1:
out_feature = ogr.Feature(def_feature)
out_feature.SetGeometry(current_union)
out_lyr.ResetReading()
out_lyr.CreateFeature(out_feature)
|
另外说一下，gdal有两个求和函数，一个是普通的Union，一个是级联UnionCascaded，不过两者间的区别我看不出有什么，实际上级联也是一位群友推荐我使用的，最后没有拿到理想结果，这里也给出级联使用的代码
1
2
3
4
5
6
7
|
def dedup(geometries):
"""Return a geometry that is the union of all geometries."""
if not geometries: return None
multi = ogr.Geometry(ogr.wkbMultiPolygon)
for g in geometries:
multi.AddGeometry(g.geometry())
return multi.UnionCascaded()
|
其实上文的角点坐标我是用来确定有效区域的最远点的，为了实现这个目的，我前期还自己造了一个轮子，结果轮子不够严谨，可能在这幅影像能跑通，另一幅就不行，主要的思路是通过迭代四周的距离影像中心最远的行列，找到最近的有值点，以四个方向的初始有值点定义为可能最远点，然后通过求距离并判断实现最远两点求解。但是有的影像波段之间的有效点不一定符合我的判断定义，加上其他众多复杂的原因，这个函数就弃用了，不过这里放出来给读者，希望能给读者一些启发思路。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
|
def get_start_end_points(in_raster, out_shp_point):
ogr.RegisterAll()
dataset = gdal.Open(in_raster)
geocd = dataset.GetGeoTransform()
row = dataset.RasterXSize #hang
line = dataset.RasterYSize #lie
band1=dataset.GetRasterBand(1)
im_datas1 = band1.ReadAsArray(0, 0, row, line)
im_datas1 = im_datas1.transpose((1,0))
for i in range(row):
if im_datas1[0,i] != 0:
lefttopx = geocd[0]
lefttopy = geocd[3] + i * geocd[5]
print(im_datas1[0,i])
break
elif im_datas1[1,i] != 0:
lefttopx = geocd[0] + 1 * geocd[1]
lefttopy = geocd[3] + i * geocd[5]
print(im_datas1[1,i])
break
for i in range(row):
if im_datas1[row-1,i] != 0:
rightdownx = geocd[0] + row * geocd[1]
rightdowny = geocd[3] + i * geocd[5]
print(im_datas1[row-1,i])
break
elif im_datas1[row-2,i] != 0:
rightdownx = geocd[0] + (row-1) * geocd[1]
rightdowny = geocd[3] + i * geocd[5]
print(im_datas1[row-2,i])
break
for i in range(line):
if im_datas1[i,0] != 0:
leftdownx = geocd[0] + i * geocd[1]
leftdowny = geocd[3]
print(im_datas1[i,0])
break
elif im_datas1[i,1] != 0:
leftdownx = geocd[0] + i * geocd[1]
leftdowny = geocd[3] + 1 * geocd[5]
print(im_datas1[i,1])
break
for i in range(line):
if im_datas1[i,line-1] != 0:
righttopx = geocd[0] + i * geocd[1]
righttopy = geocd[3] + line * geocd[5]
print(im_datas1[i,line-1])
break
elif im_datas1[i,line-2] != 0:
righttopx = geocd[0] + i * geocd[1]
righttopy = geocd[3] + (line-1) * geocd[5]
print(im_datas1[i,line-2])
break
if row <= line:
start = geometry.Point(leftdownx, leftdowny)
end = geometry.Point(righttopx, righttopy)
else:
start = geometry.Point(lefttopx, lefttopy)
end = geometry.Point(rightdownx, rightdowny)
pointshp = gpd.GeoSeries([start, end],
crs='EPSG:4326', # 指定坐标系为WGS 1984
index=['0', '1'] # 相关的索引
)
pointshp.to_file(out_shp_point,driver='ESRI Shapefile',encoding='utf-8')
del dataset
|
不过这个函数写出来被弃用，也给了我一次深刻教训——能不造轮子就不造，这真的要牢记。
后记
很大的教训还是，拒绝造轮。然后就是geopandas真的强，以前很少用到，但是这次使用给了我很深刻的印象，果然站在巨人肩膀上，看的还是远。其实中间还发现了一些深入原理层的shapefile polyline简化算法，但是作为工程实现来将，稳定性还是排在我使用的第一位，geopandas的稳定性自然不用说。其他的，留在以后文章中写吧
