原文链接:https://blog.csdn.net/weixin_40625478/article/details/106851352
本文主要目的:我们有的时候需要获取矢量数据的外接矩形范围,但是一个图层数据有好几个面要素,如果获取整体的外接矩形进行处理的话,还是会造成数据量较大;那么我们就逐个面要素进行判断其外接矩形;然后按外接矩形范围生成point数据。
1)如何使用GDAL/OGR打开矢量并输出图层数据范围和每个ploygon要素范围
2)按矩形范围生成以及像元的大小为间隔,生成point数据
1、输出每个数据的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(i) #获取字段“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])
2、按外界矩形范围生成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()
标签:point,geotransform,OGR,extent,print,import,矩形,GDAL From: https://www.cnblogs.com/suoyike1001/p/16990709.html