首页 > 其他分享 >使用GDAL/OGR打开矢量并输出每个面外界矩形范围内的point数据

使用GDAL/OGR打开矢量并输出每个面外界矩形范围内的point数据

时间:2022-12-18 18:14:00浏览次数:65  
标签:point geotransform OGR extent print import 矩形 GDAL

原文链接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

相关文章