如何使用GDAL/OGR打开矢量并输出图层数据范围和每个ploygon要素范围
本文主要目的:我们有的时候需要获取矢量数据的外接矩形范围,但是一个图层数据有好几个面要素,如果获取整体的外接矩形进行处理的话,还是会造成数据量较大;那么我们就逐个面要素进行判断其外接矩形;然后按外接矩形范围生成point数据。
1)如何使用GDAL/OGR打开矢量并输出图层数据范围和每个ploygon要素范围
2)按矩形范围生成以及像元的大小为间隔,生成point数据
0、构想
来源:
数据形式,有一个面图层数据,
一个栅格数据:
我需要对矢量面数据内,提取出每个面图层内的point数据,point数据的间隔大小是 栅格数据的像元大小; 也即是上图;需要提取的point数据为:
放大查看即是这种point数据:
1、arcmap查看数据属性信息
首先用arcmap打开的我的数据,查看数据的范围,以及数据的属性信息:
2、输出每个数据的extent外接矩形范围
获取图层数据范围使用GetExtent()函数;
获取要素范围,则需要循环每个要素,然后使用GetEnvelope()函数即可;
#如何使用GDAL/OGR打开矢量文件、读取属性表字段,并将数据范围和每个ploygon要素的范围
import ogr,sys,os
import numpy as np
os.chdir('D:/20200309/shujushouji/python+jupyternoterbook/zuoye1(1)/6栅格数据操作/shiyan')
#设置driver,并打开矢量文件
driver = ogr.GetDriverByName('ESRI Shapefile')
ds = driver.Open('polygon.shp', 0)
if ds is None:
print("Could not open", 'sites.shp')
sys.exit(1)
#获取图册
layer = ds.GetLayer()
#要素数量
numFeatures = layer.GetFeatureCount()
print("Feature count: "+str(numFeatures))
#获取范围
extent = layer.GetExtent()
print("Extent:", extent)
print("UL:", extent[0],extent[3])
print("LR:", extent[1],extent[2])
#获取要素
#feature = layer.GetNextFeature()
#循环每个要素属性
for i in range(numFeatures):
feature = layer.GetNextFeature()
#获取字段“id”的属性
id = feature.GetField('type')
#获取空间属性
print(id)
geometry = feature.GetGeometryRef()
#x = geometry.GetX()
polygonextent=geometry.GetEnvelope()
print(geometry.GetEnvelope())
#print(y)
#y = geometry.GetY()
print("UL:", polygonextent[0],polygonextent[3])
print("LR:", polygonextent[1],polygonextent[2])
3、按外界矩形范围生成point数据
下面代码是根据polygon的extent的外接矩形进行采样生成point生成;
其中涉及到知识有gdal读取栅格数据,获取其影像的左上角起始坐标和像元大小,以及按找矢量获取到的外界矩形范围的空间坐标,把空间坐标转换成 影像的行列号;然后循环生成矢量point数据,并写到矢量文件中。
# 先获取polygon.shp外界矩形 的范围,以及坐标点; 然后将坐标点转换为 栅格的行列号,只对sample_region范围内的栅格区域,进行栅格转point
%matplotlib inline
import matplotlib.pyplot as plt
from shapely.geometry import Point
import geopandas as gpd
import ogr,sys,os
from pyproj import Transformer
import numpy as np
from osgeo import gdal
os.chdir('D:/20200309/shujushouji/python+jupyternoterbook/zuoye1(1)/shiyan/polygon')
#设置driver,并打开矢量文件
driver = ogr.GetDriverByName('ESRI Shapefile')
dsshp = driver.Open('sample_regions.shp', 0)
if dsshp is None:
print("Could not open", 'sites.shp')
sys.exit(1)
#获取图册
layer = dsshp.GetLayer()
#要素数量
numFeatures = layer.GetFeatureCount()
print("Feature count: "+str(numFeatures))
def world2Pixel(geotransform, x, y):
originX = geotransform[0]
originY = geotransform[3]
pixelWidth = geotransform[1]
pixelHeight = geotransform[5]
line = int((y-originY)/pixelHeight)+1
column = int((x-originX)/pixelWidth)+1
return (line,column)
def Pixel2world(geotransform, line, column):
originX = geotransform[0]
originY = geotransform[3]
pixelWidth = geotransform[1]
pixelHeight = geotransform[5]
x = column*pixelWidth + originX - pixelWidth/2
y = line*pixelHeight + originY - pixelHeight/2
return(x,y)
ds = gdal.Open( "D:/20200309/shujushouji/python+jupyternoterbook/zuoye1(1)/6栅格数据操作/landsat8_20180523/20180523.img" )
geotransform = ds.GetGeoTransform()
# 获取 空间坐标 范围
extent = layer.GetExtent()
print("Extent:", extent)
print("UL:", extent[0],extent[3])
print("LR:", extent[1],extent[2])
# 将获取研究区polygon.shp的空间坐标,转换成 栅格数据的行列数
linex1,columny1 = world2Pixel(geotransform, extent[0],extent[3])
print(f"linex1:{linex1}")
print(f"columny1:{columny1}")
linex2,columny2 = world2Pixel(geotransform,extent[1],extent[2])
print(f"linex2:{linex2}")
print(f"columny2:{columny2}")
points = []
for i in range(linex1,linex2):
for j in range(columny1,columny2):
x,y = Pixel2world(geotransform,i+1,j+1)
point = Point((x,y))
points.append(point)
gdf = gpd.GeoDataFrame()
gdf["geometry"] = points
#将栅格数据的坐标系 赋予 感兴趣区的point数据
print(ds.GetProjection())
gdf.crs=ds.GetProjection()
gdf.to_file("D:/20200309/shujushouji/python+jupyternoterbook/zuoye1(1)/6栅格数据操作/test/d_points1.shp")
gdf.plot()