GDAL python教程基础篇(7)OGR空间计算
创始人
2025-05-29 03:46:38
0

1.空间计算

地理数据处理(geoprocessing)计算函数:
多边形(Polygon):
1、交:poly3.Intersection(poly2)
2、并:poly3.Union(poly2)
3、差:poly3.Difference(poly2)
4、补:poly3.SymmetricDifference(poly2)

下面我们来做一段练习,首先我们有如下一份数据,
在这里插入图片描述
该图层内有两个几何体,我们就用这两个集合体来进行练习

from osgeo import ogr, osrshp = r'D:/work/全国行政区划/江苏省.shp'drivier = ogr.GetDriverByName('ESRI Shapefile')ds = drivier.Open(shp)layer = ds.GetLayer(0)srs = layer.GetSpatialRef()  # 获取投影print(layer.GetFeatureCount())feat_js = layer.GetNextFeature()
feat_rec = layer.GetNextFeature()geom_js = feat_js.GetGeometryRef()
geom_rec = feat_rec.GetGeometryRef()geom_In = geom_js.Intersection(geom_rec)   # 相交
geom_un = geom_js.Union(geom_rec)          # 并集
geom_diff = geom_js.Difference(geom_rec)   # 差集
geom_sy = geom_js.SymmetricDifference(geom_rec)   # 补充#  写入shp文件
ds2 = drivier.CreateDataSource('test2.shp')
layer = ds2.CreateLayer('test', srs, geom_type=ogr.wkbPolygon)#定义字段
fieldDefn = ogr.FieldDefn('id', ogr.OFTString)
fieldDefn.SetWidth(4)
layer.CreateField(fieldDefn)#定义要素
featureDefn = layer.GetLayerDefn()
feature1 = ogr.Feature(featureDefn)
feature2 = ogr.Feature(featureDefn)
feature3 = ogr.Feature(featureDefn)
feature4 = ogr.Feature(featureDefn)#设置几何形状
feature1.SetGeometry(geom_In)
feature2.SetGeometry(geom_un)
feature3.SetGeometry(geom_diff)
feature4.SetGeometry(geom_sy)#设定某字段的数值
feature1.SetField('id', 1)
feature2.SetField('id', 2)
feature3.SetField('id', 3)
feature4.SetField('id', 4)#将feature写入layer
layer.CreateFeature(feature1)
layer.CreateFeature(feature2)
layer.CreateFeature(feature3)
layer.CreateFeature(feature4)ds2.Destroy()

得到的结果
1、两个几何体的交集
交集
2、两个几何体的并集在这里插入图片描述
3、差集(geom_js去除相交的部分后剩余的部分)
在这里插入图片描述
4、补集,并集减去相交的部分
在这里插入图片描述
几何形状(Geometry):

1.2给geometry加buffer,效果就是把点、线变成多边形,线条加粗:

.Buffer()

1.3两个geometry是否相等:.Equal()

1.4两个geometry之间的最短距离:.Distance()

1.5用方框边界框住这个几何形状,并返回四个角的点坐标(minx,maxx, miny, maxy):.GetEnvelope()

geom_js.GetENvelope()
(116.35518309500014, 121.92747178000022, 30.760280365000085, 35.127189627500144)

2.空间分析的方法

2.1.获取面积

建议采用GetArea()获得整个多边形覆盖的面积。

geom_js.GetArea()`
>>> 9.896416947183278

不过此用法以“度”为量纲,已经失去它作为面积的意义了。若要较为准确的面积,可根据空间参考将度转换为米,再将面积计算为平方米。

判断两个对象的关系

2.2Intersect语句是判断两个要素是否相交:

poly2.Intersect(poly1)

返回0表示不相交,返回1表示相交。

2.3Disjoint是判断两个要素是否不相交:

poly2.Disjoint(poly1)

返回1表示不相交,返回0表示相交,跟Intersect正好相反。

2.4Touch表示相邻:

poly2.Touches(poly1)

返回0表示不相邻,返回1表示相邻。

2.5Crosses表示横穿,一般是指一条线穿过一个多边形:

poly2.Crosses(line)

返回0表示不穿过,返回1表示穿过。

2.6Within表示包含,指一个要素完全被另一个要素包含:

ptB.Within(poly1)

返回0表示点在多边形外,返回1表示点在多边形内。

2.7Contains也表示包含,跟Within正好相反:

poly1.Contains(ptB)

与Within的主调对象与参数对换。

2.8Overlaps重叠,只有两个多边形之间才能算作overlap

poly2.Overlaps(poly3)

返回0表示不重叠,返回1表示重叠。

相关内容

热门资讯

喜欢穿一身黑的男生性格(喜欢穿... 今天百科达人给各位分享喜欢穿一身黑的男生性格的知识,其中也会对喜欢穿一身黑衣服的男人人好相处吗进行解...
发春是什么意思(思春和发春是什... 本篇文章极速百科给大家谈谈发春是什么意思,以及思春和发春是什么意思对应的知识点,希望对各位有所帮助,...
网络用语zl是什么意思(zl是... 今天给各位分享网络用语zl是什么意思的知识,其中也会对zl是啥意思是什么网络用语进行解释,如果能碰巧...
为什么酷狗音乐自己唱的歌不能下... 本篇文章极速百科小编给大家谈谈为什么酷狗音乐自己唱的歌不能下载到本地?,以及为什么酷狗下载的歌曲不是...
华为下载未安装的文件去哪找(华... 今天百科达人给各位分享华为下载未安装的文件去哪找的知识,其中也会对华为下载未安装的文件去哪找到进行解...
家里可以做假山养金鱼吗(假山能... 今天百科达人给各位分享家里可以做假山养金鱼吗的知识,其中也会对假山能放鱼缸里吗进行解释,如果能碰巧解...
四分五裂是什么生肖什么动物(四... 本篇文章极速百科小编给大家谈谈四分五裂是什么生肖什么动物,以及四分五裂打一生肖是什么对应的知识点,希...
怎么往应用助手里添加应用(应用... 今天百科达人给各位分享怎么往应用助手里添加应用的知识,其中也会对应用助手怎么添加微信进行解释,如果能...
一帆风顺二龙腾飞三阳开泰祝福语... 本篇文章极速百科给大家谈谈一帆风顺二龙腾飞三阳开泰祝福语,以及一帆风顺二龙腾飞三阳开泰祝福语结婚对应...
美团联名卡审核成功待激活(美团... 今天百科达人给各位分享美团联名卡审核成功待激活的知识,其中也会对美团联名卡审核未通过进行解释,如果能...