gdal笔记--chapter7
使用GR进行矢量分析
1.计算相交区域面积的两种方法（对要素或图层操作）
from osgeo import gdal, ogr water_ds = ogr.Open(r'') #打开矢量数据，沼泽数据 water_lyr = water_ds.GetLayer(0) #打开图层 water_lyr.SetAttributeFilter('WaterbdyID=1011327') #属性选择 marsh_feat = water_lyr.GetNextFeature() #获取要素 marsh_geom = marsh_feat.geometry().Clone() #获取要素的几何体并clone nola_ds = ogr.Open(r'') #打开城市数据 nola_lyr = nola_ds.GetLayer(0) nola_feat = nola_lyr.GetNextFeature() nola_geom = nola_feat.geometry().Clone() #相交操作 intersection = marsh_geom.Intersection(nola_geom) #遍历要素计算相交区域面积 water_lyr.SetAttributeFilter("Feature != 'LaKe'") #限定为非湖泊 water_lyr.SetSpatialFilter(nola_geom) #市区内空间过滤，减少不必要要素，减少后续处理时间 wetlands_area = 0 for feat in water_lyr: intersect = feat.geometry().Intersection(nola_geom) #交集部分 wetlands_area += intersect.GetArea() pcnt = wetlands_area/nola_geom.GetArea() print('{:.1%} of New Orleans is wetland'.format(pcnt)) #图层相交计算相交区域面积 water_lyr.SetAttributeFilter("Feature!='Lake'") memory_driver = ogr.GetDriverByName('Memory') #内存数据为临时数据 temp_ds = memory_driver.CreateDataSource('temp') temp_lyr = temp_ds.CreateLayer('temp') nola_lyr.Intersection(water_lyr, temp_lyr)#图层相交执行空间过滤，结果存储在临时图层 sql = 'SELECT SUM(OGR_GEOM_AREA) AS area FROM temp' #查询area的sql语句 lyr = temp_ds.ExecuteSQL(sql) pcnt = lyr.GetFeature(0).GetField('area')/nola_geom.GetArea() #获取area属性 print('{:.1%} of New Orleans is wetland'.format(pcnt))
2.创建缓冲区确定火山附近城市
#确定火山附近城市数量的有缺陷方法 shp_ds = ogr.Open(r'') volcano_lyr = shp_ds.GetLayer('') #打开火山和城市shp cities_lyr = shp_ds.GetLayer('cities_albers') memory_driver = ogr.GetDriverByName('memory') memory_ds = memory_driver.CreateDataSource('temp') buff_lyr = memory_ds.CreateLayer('buffer') #创建临时数据作为缓冲区 buff_feat = ogr.Feature(buff_lyr.GetLayerDefn()) #创建缓冲区要素 #直接SetGeometry会让不同火山的缓冲区有叠加的可能，使得城市可能位于多个火山缓冲区导致重复。 #逐要素添加缓冲区 for volcano_feat in volcano_lyr: buff_geom = volcano_feat.geometry().Buffer(16000) #缓冲区操作 tmp = buff_feat.SetGeometry(buff_geom) #要素添加几何体 tmp = buff_lyr.CreateFeature(buff_feat) #图层添加要素 result_lyr = memory_ds.CreateLayer('result') buff_lyr.Intersection(cities_lyr, result_lyr) #求交集 print('Cities:{}'.format(result_lyr.GetFeatureCount())) #更好的方法 volcano_lyr.ResetReading() #回到图层起始位置 multipoly = ogr.Geometry(ogr.wkbMultiPolygon) #创建复合多边形 for volcano_feat in volcano_lyr: buff_geom = volcano_feat.geometry().Buffer(16000) multipoly.AddGeometry(buff_geom) #复合多边形添加缓冲区 #UnionCascaded方法把缓冲区合在一起，相同部分合并，避免了缓冲区的重复 cities_lyr.SetSpatialFilter(multipoly.UnionCascaded()) #UnionCascaded合并多边形 print('Cities:{}'.format(cities_lyr.GetFeatureCount()))
3.计算某火山和某城市距离
#计算某火山和某城市的距离（round计算几何体的距离） volcano_lyr.SetAttributeFilter("NAME='Rainier'") feat = volcano_lyr.GetNextFeature() rainier = feat.geometry().Clone() #属性过滤获取lyr范围，读取要素，读取几何体并Clone cities_lyr.SetAttributeFilter("NAME='Seattle'") feat = cities_lyr.GetNextFeature() seattle = feat.geometry().Clone() meters = round(rainier.Distance(seattle)) #计算几何体距离 miles = meters/1600 #英里 print('{} meters {} miles'.format(meters, miles))
4.各种操作忽略Z值
pt1_2d = ogr.Geometry(ogr.wkbPoint) pt1_2d.AddPoint(15, 15) pt2_2d = ogr.Geometry(ogr.wkbPoint) pt2_2d.AddPoint(15, 19) print(pt1_2d.Distance(pt2_2d)) pt1_25d = ogr.Geometry(ogr.wkbPoint) pt1_25d.AddPoint(15, 15, 0) pt2_25d = ogr.Geometry(ogr.wkbPoint) pt2_25d.AddPoint(15, 19, 3) print(pt1_25d.Distance(pt2_25d)) #结果相同，未对Z运算 ring = ogr.Geometry(ogr.wkbLinearRing) #创建环 ring.AddPoint(10,10) ring.AddPoint(10,20) ring.AddPoint(20,20) ring.AddPoint(20,10) poly_2d = ogr.Geometry(ogr.wkbPolygon) #创建多边形 poly_2d.AddGeometry(ring) poly_2d.CloseRings() #多边形闭合 print(poly_2d.GetArea()) #2.5d同忽略Z ring = ogr.Geometry(ogr.wkbLinearRing) ring.AddPoint(10, 10, 0) ring.AddPoint(10, 20, 0) ring.AddPoint(20, 20, 10) ring.AddPoint(20, 10, 10) poly_25d = ogr.Geometry(ogr.wkbPolygon25D) poly_25d.AddGeometry(ring) poly_25d.CloseRings() print(poly_25d.GetArea()) poly_2d.Contains(pt1_2d) poly_25d.Contains(pt1_2d) #叠加操作忽略Z
5.案例1：风力场选址
#7.3风力场选址 census_fn = r'' #人口普查数据 census_ds = ogr.Open(census_fn, True) census_lyr = census_ds.GetLayer() density_field = ogr.FieldDefn('popsqkm', ogr.OFTReal) #字段定义 census_lyr.CreateField(density_field) #给图层创建新字段 for row in census_lyr: #遍历图层要素 pop = row.GetField('HD01_S001') #读取pop字段值 sqkm = row.Geometry().GetArea()/1000000 #要素几何体获取面积 row.SetField('popsqkm', pop/sqkm) #设置popsqkm字段值 census_lyr.SetFeature(row) #图层添加要素 country_fn = r'' country_ds = ogr.Open(r'') country_lyr = country_ds.GetLayer() country_lyr.SetAttributeFilter("COUNTRY='Imperial Country'") country_row = next(country_lyr) #读取要素 country_geom = country_row.geometry().Clone() del country_ds country_geom.TransformTo(census_lyr.GetSpatialRef()) #投影转换 census_lyr.SetSpatialFilter(country_geom) census_lyr.SetAttributeFilter('popsqkm<0.5') #空间和属性过滤完的census_lyr wind_fn = r'' wind_ds = ogr.Open(wind_fn) wind_lyr = wind_ds.GetLayer() wind_lyr.SetAttributeFilter('WPC>=3') #打开风力数据属性过滤 out_fn = r'' out_ds = ogr.GetDriverByName('ESRI Shapefile').CreateDataSource(out_fn) #创建输出文件 out_lyr = out_ds.CreateLayer('wind_farm', wind_lyr.GetSpatialRef(), ogr.wkbPolygon) out_lyr.CreateField(ogr.FieldDefn('wind', ogr.OFTInteger)) #创建字段 out_lyr.CreateField(ogr.FieldDefn('popsqkm', ogr.OFTReal)) out_row = ogr.Feature(out_lyr.GetLayerDfen()) #新建空要素 #将人口普查数据与风力数据相交 for census_row in census_lyr: census_geom = census_row.geometry() census_geom = census_geom.Intersection(country_geom) #人口普查数据和县界相交 wind_lyr.SetSpatialFilter(census_geom)#空间过滤风力lyr print('Intersecting census tract with {0} wind polygons'.format(wind_lyr.GetFeatureCount())) if wind_lyr.GetFeatureCount()>0: out_row.SetField('popsqkm', census_row.GetField('popsqkm')) for wind_row in wind_lyr: wind_geom = wind_row.geometry() if census_geom.Intersect(wind_geom): #人口普查数据和风力多边形相交 new_geom = census_geom.Intersection(wind_geom) out_row.SetField('wind', wind_row.GetField('WPC')) out_row.SetGeometry(new_geom) out_lyr.CreateFeature(out_row) #将小的多边形合并成大的多边形 folder = '' ds = ogr.Open(folder, True) in_lyr = ds.GetLayerByName('wind_farm') #源图层 out_lyr = ds.CreateLayer('wind_farm2', in_lyr.GetSpatialRef(), ogr.wkbPolygon) #输出图层 out_row = ogr.Feature(out_lyr.GetLayerDefn()) #新建输出空要素 multipoly = ogr.Geometry(ogr.wkbMultiPolygon) #创建复合多边形 for in_row in in_lyr: #遍历图层每个要素 in_geom = in_row.geometry().Clone() in_geom_type = in_geom.GetGeometryType() #查询类型是多边形还是复合多边形 if in_geom_type == ogr.wkbPolygon: multipoly.AddGeometry(in_geom) #多边形则添加到新复合多边形 elif in_geom_type == ogr.wkbMultiPolygon: for i in range(in_geom.GetGeometryCount()): #复合多边形，则in_geom.GetGeometryCount()遍历 multipoly.AddGeometry(in_geom.GetGeometryRef(i)) #打断成多边形后添加到新复合多边形 multipoly = multipoly.UnionCascaded() #合并操作 for i in range(multipoly.GetGeometryCount()): #遍历新复合多边形 poly = multipoly.GetGeometryRef(i) if poly.GetArea() > 1000000: #保留面积大于...的多边形几何体 out_row.SetGeometry(poly) #输出要素添加多边形几何体 out_lyr.CreateFeature(out_row) #图层添加要素 del ds
6.案例2：动物跟踪数据
#7.4动物跟踪数据 #从csv文件新建shapefile文件 from osgeo import ogr, osr csv_fn = r'' shp_fn = r'' sr = osr.SpatialReference(osr.SRS_WKT_WGS84) #设置shpfile投影信息 #新建shpfile文件并创建字段存储时间、日期，新建要素 shp_ds = ogr.GetDriverByName('ESRI Shapefile').CreateDataSource(shp_fn) shp_lyr = shp_ds.CreateLayer('albatross_dd', sr, ogr.wkbPoint) shp_lyr.CreateField(ogr.FieldDefn('tag_id', ogr.OFTString)) shp_lyr.CreateField(ogr.FieldDefn('timestamp', ogr.OFTString)) shp_row = ogr.Feature(shp_lyr.GetLayerDefn()) csv_ds = ogr.Open(csv_fn) #ogr打开scv文件 csv_lyr = csv_ds.GetLayer() #csv文件获取图层 for csv_row in csv_lyr: x = csv_row.GetFieldAsDouble('location-long') y = csv_row.GetFieldAsDouble('location-lat') shp_pt = ogr.Geometry(ogr.wkbPoint) shp_pt.AddPoint(x, y) tag_id = csv_row.GetField('individual-local-identifier') timestamp = csv_row.GetField('timestamp') shp_row.SetGeometry(shp_pt) shp_row.SetField('tag_id', tag_id) shp_row.SetField('timestamp', timestamp) shp_lyr.CreateFeature(shp_row) del csv_ds, shp_ds #有坏点需要删除 shp_ds = ogr.Open(shp_fn, True) shp_lyr = shp_ds.GetLayer() shp_lyr.SetSpatialFilterRect(-1, -1, 1, 1) #设置坏点的空间范围 for shp_row in shp_lyr: shp_lyr.DeleteFeature(shp_row.GetFID()) #获取FID删除要素 shp_lyr.SetSpatialFilter(None) #不再设置空间过滤 shp_ds.ExecuteSQL('REPACK ' + shp_lyr.GetName()) #回收空间 shp_ds.ExecuteSQL('RECOMPUTE EXTENT ON' + shp_lyr.GetName())#更新元数据空间范围 del shp_ds #投影转换，转换成兰伯特投影米制单位 #ogr2ogr -f 'ESRI Shapefile' -t_srs '+proj=lcc.....' albatross_lambert.shp yuanshi.shp #从属性列获取唯一值（从图层获取指定字段值） def get_unique(datasource, layer_name, field_name): sql = 'SELECT DISTINCT {0} FROM {1}'.format(field_name, layer_name) #搜寻图层中对应的某一字段值 lyr = datasource.ExecuteSQL(sql) values = [] #存储图层要素对应的字段值 for row in lyr: values.append(row.GetField(field_name)) #获取该鸟所有时间的字段值 datasource.ReleaseResultSet(lyr) #ReleaseResultSet(<result_layer>)是将查询结果释放，在执行下一条SQL语句之前一定要先释放。 return values #计算相邻点间距离 ds = ogr.Open(r'') lyr = ds.GetLayerByName('albatross_lambert') lyr.CreateField(ogr.FieldDefn('distance', ogr.OFTReal)) #创建举例子段 tag_ids = get_unique(ds, 'albatross_lambert', 'tag_id') #获取对应ID鸟的字段值 for tag_id in tag_ids: #遍历所有鸟种tag_id print('Processing ' + tag_id) lyr.SetAttributeFilter("tag_id='{}'".format(tag_id)) #属性查询 row = next(lyr) #获取第一个要素 previous_pt = row.geometry().Clone() #要素几何体 previous_time = row.GetField('timestamp') #该要素对应时间 for row in lyr: #遍历图层要素（单个鸟的点位） current_time = row.GetField('timestamp') #获取当前时间 if current_time < previous_time: raise Exception('Timestamps out of order') current_pt = row.geometry().Clone() distance = current_pt.Distance(previous_pt) #计算距离 row.SetField('distance', distance) #设置distance字段值 lyr.SetFeature(row) #图层添加要素 previous_pt = current_pt #更新点位信息 previous_time = current_time del ds #查找地点间最大飞行速度和所耗时间 from datetime import datetime from osgeo import ogr date_format = '%Y-%M-%d %H:%M:%S.%f' ds = ogr.Open(r'') lyr = ds.GetLayerByName('albatross_lambert') for tag_id in get_unique(ds, 'albatross', 'tag_id'): #遍历鸟种 max_speed = 0 #初始化最大速度 lyr.SetAttributeFilter("tag_id={}".format(tag_id)) #属性查询 row = next(lyr) #第一个时间点 ts = row.GetField('timestamp') previous_time = datetime.strptime(ts, date_format) for row in lyr: ts = row.GetField('timestamp') current_time = datetime.strptime(ds, date_format) #转换成可计算的时间格式 elapsed_time = current_time - previous_time hours = elapsed_time.total_seconds()/3600 distance = row.GetField('distance') speed = distance/hours #计算speed max_speed = max(max_speed, speed) #更新max_speed print('Max speed for {0}:{1}'.format(tag_id, max_speed)) #为每只鸟创建凸包多边形 ds = ogr.Open(r'') pt_lyr = ds.GetLayerByName('albatross_lambert') poly_lyr = ds.CreateLayer('albatross_ranges', pt_lyr.GetSpatialRef(), ogr.wkbPoint) #输出图层 id_field = ogr.FieldDefn('tag_id', ogr.OFTString) area_field = ogr.FieldDefn('tag_id', ogr.OFTString) area_field.SetWidth(30) area_field.SetPrecision(4) poly_lyr.CreateFields([id_field, area_field]) #新建id\area字段 poly_row = ogr.Feature(poly_lyr.GetLayerDefn()) for tag_id in get_unique(ds, 'albatross_lambert', 'tag_id'): print('Processing ' + tag_id) pt_lyr.SetAttributeFilter("tag_id = '{}'".format(tag_id)) #为lyr进行属性过滤 locations = ogr.Geometry(ogr.wkbMultiPoint) #为每只鸟(tag_id)创建复合多边形几何体 for pt_row in pt_lyr: #遍历原图层要素，把单个ID的所有点添加到一个几何对象 locations.AddGeometry(pt_row.geometry().Clone()) #把原要素几何对象添加到新几何对象 homerange = locations.ConvexHull() #所有点创建凸包多边形 poly_row.SetGeometry(homerange) #poly要素设置几何体 poly_row.SetField('tag_id', tag_id) poly_row.SetField('area', homerange.GetArea()) poly_lyr.CreateFeature(poly_row) #创建poly要素 del ds #创建地理区域分割的凸包多边形（有点难） #核心思想：不断读入要素，pt_locations用于存储连续相同位置的点，一旦位置发生改变pt_locations置0。 #此外，用poly_lyr不断添加凸多边形。 #疑惑：这样看起来如果点位是3个陆地、岛屿、陆地的情况下，pt_locations会新建，那么先前同类型的点就清空了。这样添加到最后每侧应该是若干个的凸包多边形吧？没有成为整体。这是为什么？ #自问自答：可能因为点位本身就是按照时间顺序添加的，不可能一下子在岛屿，一下又飞到了陆地，所以点位置是连续的，这样最后就是两个大凸多边形。 land_lyr = ds.GetLayer('land_lambert') land_row = next(land_lyr) land_poly = land_row.geometry().Buffer(100000) #创建陆地缓冲区 for tag_id in get_unique(ds, 'albatross_lambert', 'tag_id'): #遍历每只鸟 print('Processing ' + tag_id) pt_lyr.SetAttributeFilter("tag_id = '{}".format(tag_id)) #属性过滤id pt_locations = ogr.Geometry(ogr.wkbMultiPoint) #创建多点 last_location = None for pt_row in pt_lyr: #遍历图层的每个点 pt = pt_row.geometry().Clone() if not land_poly.Contains(pt): #如果点不在缓冲区，跳过这个点 continue if pt.GetX() < -2800000: #岛屿区域 location = 'island' else: #大陆区域 location = 'mainland' if location != last_location: #如果当前点位不同于上个点位 if pt_locations.GetGeometryCount() > 2: #除当前点以外的点，至少有3个，能组成多边形时（位置发生变化时，创建多边形） homerange = pt_locations.ConvexHull() #凸包多边形 poly_row.SetGeometry(homerange) #设置几何体和属性值 poly_row.SetField('tag_id', tag_id) poly_row.SetField('area', homerange.GetArea()) poly_row.SetField('location', last_location) poly_lyr.CreateFeature(poly_row) #添加要素，这里poly_lyr把相同的点组成的凸多边形不断地添加到里面 #点数不够3直接新建多点对象，添加点并更新点位 pt_locations = ogr.Geometry(ogr.wkbMultiPoint) #新建操作用于在点位发生变化时重新计算凸包多边形，避免不同位置点存在一起 last_location = location #点位相同，那么继续添加多点对象，直到不再相同再进行凸包操作 pt_locations.AddGeometry(pt) if pt_locations.GetGeometryCount() > 2: #由于上面只在点位发生变化的时候创建多边形，所以需要单独保存最后一个组点创建多边形 homerange = pt_locations.ConvexHull() poly_row.SetGeometry(homerange) poly_row.SetField('tag_id', tag_id) poly_row.SetField('area', homerange.GetArea()) poly_row.SetField('location', last_location) poly_lyr.CreateFeature(poly_row) #计算所有岛屿访问时公共区域所占面积的比例 ds = ogr.Open(r'') lyr = ds.GetLayerByName('albatross_ranges2') lyr.SetAttributeFilter("tag_id = '1163-1163' and location = 'island'") #属性过滤，id和island区域 row = next(lyr) #读取第一个要素 all_areas = row.geometry().Clone() #读取第一个要素，初始化all_areas\common_areas common_areas = row.geometry().Clone() for row in lyr: #遍历图层的要素，作交集和并集运算 all_areas = all_areas.Union(row.geometry()) #更新并集区域 common_areas = common_areas.Intersection(row.geometry()) #更新交集区域 percent = common_areas.GetArea()/all_areas.GetArea()*100 #面积百分比 print('Percent of all area used in every visit: {0}'.format(percent))
