GDAL
- 1 gdal库
- 2 栅格驱动
- 3 栅格数据集（就是包含各种栅格属性的一个类）
- 3.1 坐标（6个参数）
- 3.1.2 tif文件的地理坐标（两种情况）
- 3.2 波段数、大小、投影等信息
- 3.3 读取栅格像元
- 3.4 创建栅格影像
- 3.4.1 直接用数组创建数据集
- 3.4.2 用CreateCopy直接复制现有的数据集
- 3.4.3 分块读取（解决大文件读取慢的问题）
- 3.4.4 随机裁剪栅格（制作深度学习样本数据）
- 3.4.5 计算NDVI波段
- 4 矢量数据处理（OGR库）
- 4.1 读取矢量文件
- 4.2 创建点要素
- 4.3 创建线要素（和点要素步骤一样）
- 4.4 创建面要素（和点要素步骤一样）
- 4.5 选择要素
- 4.5.1 按属性信息选择要素
- 4.5.2 按空间位置或sql语句选择要素
- 4.6 栅格矢量化，并将线段进行平滑操作
- 5 一些处理工具
- 6 总结（简单，半天就学会了）
- 参考资料
1 gdal库
2 栅格驱动
读取栅格数据的流程：
（1）获取栅格驱动(会有1个默认的，可以不创建，最好创建上，看着清晰)；
（2）驱动进行注册（会有默认值，可以不创建，创建的话可以选择只读，或者读写都可以，相当于给栅格驱动设置一个权限）；
（3）通过驱动获取数据集；
（4）通过数据集操作栅格属性；
try:
from osgeo import gdal,ogr # gdal是处理栅格，ogr处理矢量
except:
import gdal
driverCount = gdal.GetDriverCount()
print(driverCount)
#gdal.AllRegister()
driver = gdal.GetDriverByName('GTiff')
driver.Register()
print(driver.ShortName)
print(driver.LongName)
print(dir(driver))
3 栅格数据集（就是包含各种栅格属性的一个类）
3.1 坐标（6个参数）
# 读取地理坐标
from osgeo import gdal
filePath = '1.tif' # tif文件路径
dataset = gdal.Open(filePath) # 打开tif
adfGeoTransform = dataset.GetGeoTransform() # 读取地理信息
# 左上角地理坐标
print(adfGeoTransform[0])
print(adfGeoTransform[3])
nXSize = dataset.RasterXSize # 列数
nYSize = dataset.RasterYSize # 行数
print(nXSize, nYSize)
arrSlope = [] # 用于存储每个像素的（X，Y）坐标
for i in range(nYSize):
row = []
for j in range(nXSize):
px = adfGeoTransform[0] + i * adfGeoTransform[1] + j * adfGeoTransform[2]
py = adfGeoTransform[3] + i * adfGeoTransform[4] + j * adfGeoTransform[5]
col = [px, py] # 每个像素的经纬度
row.append(col)
print(col)
arrSlope.append(row)
上面的代码其实已经实现获取tif中经纬度，如果大家仔细研究一下会发现，其实我们使用的就是gdal里面的GetGeoTransform方法读取坐标，简单介绍一下该方法，该方法会返回以下六个参数
GT(0) 左上像素左上角的x坐标。
GT(1) w-e像素分辨率/像素宽度。
GT(2) 行旋转（通常为零）。
GT(3) 左上像素左上角的y坐标。
GT(4) 列旋转（通常为零）。
GT(5) n-s像素分辨率/像素高度（北上图像为负值）
3.1.2 tif文件的地理坐标（两种情况）
一是只有tif文件，地理坐标存储在这个tif文件里；
二是有附带的tfw文件，这个文件里存储的地理坐标的6个参数；
arcgis读取tif文件的时候首先看tif文件里有没有地理坐标，有的话直接用，没有的话才会去找tfw文件；
3.2 波段数、大小、投影等信息
try:
from osgeo import gdal,ogr
except:
import gdal
gdal.AllRegister()
dataset = gdal.Open('C:\\Users\\Administrator\\Desktop\\gdal\\样例数据\\GF1\\GF1.tif')
print(dir(dataset))
#波段信息
bandCount = dataset.RasterCount
print(bandCount)
#大小
weight,height = dataset.RasterXSize,dataset.RasterYSize
print(weight,height)
#空间信息
transform = dataset.GetGeoTransform()
print(transform)
#投影信息
proj = dataset.GetProjection()
#描述信息
descrip = dataset.GetDescription()
print(descrip)
3.3 读取栅格像元
读取栅格数据的流程：
（1）获取栅格驱动(会有1个默认的，可以不创建)；
（2）驱动进行注册（会有默认值，可以不创建，创建的话可以选择只读，或者读写都可以，相当于给栅格驱动设置一个权限）；
（3）通过驱动获取数据集；
（4）通过数据集操作栅格属性；
try:
from osgeo import gdal,ogr
except:
import gdal
gdal.AllRegister() # 只读，这里前面没有指定驱动，它自动会用tif
dataset = gdal.Open('C:\\Users\\Administrator\\Desktop\\gdal\\样例数据\\GF1\\GF1.tif')
width = dataset.RasterXSize # 数据集中的各种属性信息
height = dataset.RasterYSize
bandCont = dataset.RasterCount
for i in range(bandCont):
band = dataset.GetRasterBand(i+1) # 获取波段，波段从1开始，不是从0开始
bandinform = band.ReadRaster(0,0,3,3) # 获取栅格数值0-3
print(bandinform)
#获得波段的数值类型
dataType = band.DataType
print(dataType)
#获得影像的nodata
nodata = band.GetNoDataValue()
print(nodata)
#获得最大值最小值，这个波段中的最大最小值
maxmin = band.ComputeRasterMinMax()
print(maxmin)
imageArray = band.ReadAsArray(0,0,3,3,10,10) # 读取栅格像元并输出为numpy的array形式
print(imageArray)
3.4 创建栅格影像
创建栅格影像流程：
（1）创建驱动；
（2）定义数据集名字，数据集影像宽和高；
（3）创建数据集；
（4）添加栅格像元值；
上面两种读写方式在下面的例子中都有实际运用：
（1）WriteRaster: 看3.4.4 随机裁剪栅格；
（2）WriteArray: 看3.4.5 计算NDVI波段。
3.4.1 直接用数组创建数据集
像元值存储的类型，不同数值类型得到的像元值范围不一样。
from osgeo import gdal
import numpy as np
srcFile = 'C:\\Users\\Administrator\\Desktop\\gdal\\样例数据\\GF1\\GF1.tif'
dstFile = 'C:\\Users\\Administrator\\Desktop\\gdal\\样例数据\\GF1\\GF1Create.tif'
dataSet = gdal.Open(srcFile)
width = dataSet.RasterXSize
height = dataSet.RasterYSize
bandCount = dataSet.RasterCount
driver = gdal.GetDriverByName('GTiff')
print(dir(driver))
metadata = driver.GetMetadata()
if gdal.DCAP_CREATE in metadata and metadata['DCAP_CREATE'] == 'YES':
print('支持create方法')
else:
print('不支持create方法')
if gdal.DCAP_CREATECOPY in metadata and metadata['DCAP_CREATECOPY'] == 'YES':
print('支持createCopy方法')
else:
print('不支持createCopy方法')
dataSetOut = driver.Create(dstFile,width,height,bandCount,gdal.GDT_Byte) # 8位无符号整数
#空间参考
srcProj = dataSet.GetProjection()
srcTranform = dataSet.GetGeoTransform()
dataSetOut.SetProjection(srcProj)
dataSetOut.SetGeoTransform(srcTranform)
datas= []
for i in range(bandCount):
band = dataSet.GetRasterBand(i+1)
dataArray = band.ReadAsArray(0,0,width,height)
#图像处理
datas.append(dataArray)
#bandOut = dataSetOut.GetRasterBand(i+1) # 1这种是分波段写入
#bandOut.WriteArray(dataArray)
datas = np.concatenate(datas)
dataSetOut.WriteRaster(0,0,width,height,datas.tobytes()) # 2这种是直接写入数据集
gdal.SetConfigOption('USE_RRD','YES')
dataSetOut.BuildOverviews(overviewlist = [2,4,8,16,32,64,128]) # 生成金字塔，不同分辨率的显示
3.4.2 用CreateCopy直接复制现有的数据集
from osgeo import gdal
import numpy as np
srcFile = 'C:\\Users\\Administrator\\Desktop\\gdal\\样例数据\\GF1\\GF1.tif'
dataSet = gdal.Open(srcFile)
dstFile = 'C:\\Users\\Administrator\\Desktop\\gdal\\样例数据\\GF1\\GF1Copy.tif'
driver = gdal.GetDriverByName('GTiff')
proc = gdal.TermProgress_nocb # 这个是显示进度条，可以不加这句话
datasetCopy = driver.CreateCopy(dstFile,dataSet,0,callback = proc)
print()
3.4.3 分块读取（解决大文件读取慢的问题）
# 分块读取
from osgeo import gdal
import numpy as np
gdal.AllRegister()
src = 'C:\\Users\\Administrator\\Desktop\\gdal\\样例数据\\GF1\\GF1.tif'
dataset = gdal.Open(src)
width = dataset.RasterXSize
heigth = dataset.RasterYSize
bandCount = dataset.RasterCount
block = 64
# 每次读取1块儿内容
for i in range(0,width,block):
if i+block < width:
numCols = block
else:
numCols = width-i
for j in range(0,heigth,block):
if j+block < heigth:
numRows = block
else:
numRows = heigth - j
for m in range(bandCount):
band = dataset.GetRasterBand(m+1)
imageBlock = band.ReadAsArray(i,j,numCols,numRows)
#之后进行波段数据的处理
3.4.4 随机裁剪栅格（制作深度学习样本数据）
from osgeo import gdal
import numpy as np
def ReadImage(src):
dataset = gdal.Open(src)
if dataset is None:
print('数据集无法打开！！')
width = dataset.RasterXSize # 宽
heigth = dataset.RasterYSize # 高
bandCount = dataset.RasterCount # 波段数
proj = dataset.GetProjection() # 投影信息
transform = dataset.GetGeoTransform() # 空间参考6个参数
return dataset,width,heigth,bandCount,proj,transform
def WriteImage(clipImageData,outName,imageSize,bandCount,proj,transform):
driver = gdal.GetDriverByName('GTiff')
dataset = driver.Create(outName,imageSize,imageSize,bandCount,gdal.GDT_Byte)
dataset.SetProjection(proj) # 设置投影
dataset.SetGeoTransform(transform) # 设置空间参考6要素
dataset.WriteRaster(0,0,imageSize,imageSize,clipImageData.tobyte()) # 写入栅格值
dataset.FlushCache() # 刷新到硬盘
dataset = None
if __name__ == '__main__':
src_image = 'C:\\Users\\Administrator\\Desktop\\gdal\\样例数据\\GF1\\GF1.tif'
outPath = ''
src_label = ''
dataset,width,heigth,bandCount,proj, transform = ReadImage(src_image)
datasetLabel, width, heigth, bandCount, proj, transform = ReadImage(src_label)
#裁剪的大小和数量
sampleCount = 50
imageSize = 512
count = 0
while count < sampleCount:
#生成随机的角点
randomWidth = np.random.randint(0,width-imageSize-1)
randomHeight = np.random.randint(0, heigth - imageSize - 1)
#裁剪
clipImageData = dataset.ReadAsArray(randomWidth,randomHeight,imageSize,imageSize)
clipLabelData = datasetLabel.ReadAsArray(randomWidth,randomHeight,imageSize,imageSize)
count+=1
#保存影像
outName = outPath+'\\'+'%s.tif'%(count)
WriteImage(clipImageData,outName,imageSize,bandCount,proj,transform)
WriteImage(clipLabelData, outName, imageSize, bandCount, proj, transform)
3.4.5 计算NDVI波段
from osgeo import gdal
import numpy as np
dataset = gdal.Open('C:\\Users\\Administrator\\Desktop\\gdal\\样例数据\\GF1\\GF1-NIR.tif')
width = dataset.RasterXSize
heigth = dataset.RasterYSize
bandCount = dataset.RasterCount
if bandCount < 4:
print('请确认是否存在近红外波段')
#获得相应的波段
bandNir = dataset.GetRasterBand(1)
bandNirArray = bandNir.ReadAsArray(0,0,width,heigth).astype(np.float16)
bandRed = dataset.GetRasterBand(3)
bandRedArray = bandRed.ReadAsArray(0,0,width,heigth).astype(np.float16)
#计算NDVI
NDVI=(bandNirArray - bandRedArray)/(bandNirArray+bandRedArray)
#排除分母为0情况
mask = np.greater(bandNirArray + bandRedArray,0)
NDVI = np.choose(mask,(-99,NDVI))
#ndvi写出去
outName = 'C:\\Users\\Administrator\\Desktop\\gdal\\样例数据\\GF1\\NDVI_FF.tif'
driver = gdal.GetDriverByName('GTiff')
datasetOut =driver.Create(outName,width,heigth,1,gdal.GDT_Float32)
datasetOut.GetRasterBand(1).WriteArray(NDVI)
datasetOut.FlushCache()
datasetOut = None
4 矢量数据处理（OGR库）
4.1 读取矢量文件
from osgeo import gdal
from osgeo import ogr
import os
##打开矢量数据集
os.chdir('C:\\Users\\Administrator\\Desktop\\gdal\\样例数据\\shp')
driver = ogr.GetDriverByName('ESRI Shapefile')
dataset = driver.Open('prov_capital.shp',0) # 0代表只读
if dataset is None:
print('矢量数据打开出错！！！')
#打开图层
layer = dataset.GetLayer(0) # 获取第1个图层，当然shp格式只有1个图层
#获取图层的属性信息
#图层里有共多少要素、图层的空间范围、属性表的结构信息
featureCount = layer.GetFeatureCount()
print(featureCount)
#西东南北及空间参考信息
extent = layer.GetExtent()
print(extent)
proj = layer.GetSpatialRef() # 获取空间投影
print(proj)
#获得属性表的信息，属性的字段，就是每列的列名等信息
layerDef = layer.GetLayerDefn() # 获取属性表
layerDefCount = layerDef.GetFieldCount() # 获取字段数量
for i in range(layerDefCount):
defn = layerDef.GetFieldDefn(i) # 获取这个字段
# 输出字段名字，字段类型，精度等
print(defn.GetName(),defn.GetWidth(),defn.GetType(),defn.GetPrecision())
#获取要素及要素信息
for i in range(featureCount):
feature = layer.GetNextFeature() # 获取要素
geo = feature.GetGeometryRef() # 获取几何对象
#获取指定字段的内容
name = feature.GetField('name')
#获得点的坐标
X = geo.GetX()
Y = geo.GetY()
dataset.Destroy() #释放内存
4.2 创建点要素
分为一下六步：
（1）获取驱动；
（2）创建数据集；
（3）创建图层；
（4）创建属性表；
（5）创建要素并添加到图层中(包括集合信息和属性信息）；
（6）保存到磁盘。
from osgeo import gdal,osr,ogr
import os
out_shp = 'C:\\Users\\Administrator\\Desktop\\gdal\\样例数据\\4.2\\test.shp'
driver = ogr.GetDriverByName('ESRI Shapefile')
#创建数据集
ds = driver.CreateDataSource(out_shp)
#判断文件夹下是否已经存在同名文件
if os.path.exists(out_shp):
driver.DeletDataSource(out_shp)
#创建图层
# 方法的3个参数，图层名，空间参考，图层类型
layer = ds.CreateLayer('point',None ,geom_type = ogr.wkbPoint)
#定义属性结构，字段一共有3个
fieldID = ogr.FieldDefn('id',ogr.OFTString)
fieldID.SetWidth(4)
layer.CreateField(fieldID)
fieldName = ogr.FieldDefn('x',ogr.OFTString)
fieldName.SetWidth(10)#设置宽度
layer.CreateField(fieldName)
fieldName = ogr.FieldDefn('y',ogr.OFTString)
fieldName.SetWidth(10)#设置宽度
layer.CreateField(fieldName)
#设置地理位置
pointList = [(100,100),(200,200),(300,300),(400,400)]
#point = ogr.Geometry(ogr.wkbPoint)
# 创建点要素
# 从layer中获得要素的结构
featureDefn = layer.GetLayerDefn()
feature = ogr.Feature(featureDefn)
for points in pointList:
ss = str(pointList.index(points))
# 设置要素的属性信息和几何信息
feature.SetField('id', ss)
feature.SetField('x', points[0])
feature.SetField('y', points[1])
#point.AddPoint(points[0],points[1])
#feature.SetGeometry(point)
#使用wkt创建要素
# wkt = 'Point(%f %f)'%(points[0],points[1])
# geo = ogr.CreateGeometryFromWkt(wkt)
# feature.SetGeometry(geo)
# 创建，将要素添加到图层中
layer.CreateFeature(feature)
#释放内存，刷新到磁盘中
ds.Destroy()
4.3 创建线要素（和点要素步骤一样）
from osgeo import ogr
import os
out_shp = 'C:\\Users\\Administrator\\Desktop\\gdal\\样例数据\\4.2\\line.shp'
driver = ogr.GetDriverByName('ESRI Shapefile')
ds = driver.CreateDataSource(out_shp)
#判断文件夹下是否已经存在同名文件
if os.path.exists(out_shp):
driver.DeletDataSource(out_shp)
layer = ds.CreateLayer('line',None,ogr.wkbLineString)
#定义属性表结构
fieldID = ogr.FieldDefn('id',ogr.OFTString)
fieldID.SetWidth(4)
layer.CreateField(fieldID)
#创建要素
featureDefn = layer.GetLayerDefn()
feature = ogr.Feature(featureDefn)
pointList = [(100,100),(100,500)]#每两个点连成一个线
#属性信息
feature.SetField('id',0)
#几何信息
line = ogr.Geometry(ogr.wkbLineString)
#wkt = '' # 用wkt字符串也可以创建几何，是第2种创建几何的方式
for idx,point in enumerate(pointList):
line.AddPoint(point[0],point[1])
# if idx == len(pointList)-1:
# wkt = wkt + '%f %f' % (point[0], point[1])
# else:
# wkt = wkt + '%f %f,'%(point[0],point[1])
# wkt = 'Linestring(%s)'% wkt
# geo = ogr.CreateGeometryFromWkt(wkt)
# feature.SetGeometry(geo)
# 下面这个代码的作用是使得线闭合，不加这个代码，得到的是不闭合的线，首尾不相连
line.CloseRings()
feature.SetGeometry(line)
layer.CreateFeature(feature)
ds.Destroy()
4.4 创建面要素（和点要素步骤一样）
from osgeo import ogr
out_shp = 'C:\\Users\\Administrator\\Desktop\\gdal\\样例数据\\4.2\\polygon.shp'
driver = ogr.GetDriverByName('ESRI Shapefile')
ds = driver.CreateDataSource(out_shp)
layer = ds.CreateLayer('polygon',None,ogr.wkbPolygon)
#属性结构
fieldID = ogr.FieldDefn('id',ogr.OFTString)
fieldID.SetWidth(4)
layer.CreateField(fieldID)
#创建要素
featureDefn = layer.GetLayerDefn()
feature = ogr.Feature(featureDefn)#创建一个对象
#设置feature对象的属性
feature.SetField('id','1')
#设置几何属性
#1、封闭的线
pointList = [(100,100),(100,500),(500,500),(500,100),(100,100)]
#wkt = 'Polygon((100 100,100 500,500 500,500 100,100 100))' # wkt字符串创建几何的方式
lineClosing = ogr.Geometry(ogr.wkbLinearRing)
for point in pointList:
lineClosing.AddPoint(point[0],point[1])
lineClosing.CloseRings()
polygon = ogr.Geometry(ogr.wkbPolygon)
polygon.AddGeometry(lineClosing)
#polygon = ogr.CreateGeometryFromWkt(wkt)
feature.SetGeometry(polygon)
layer.CreateFeature(feature)
ds.Destroy()
4.5 选择要素
4.5.1 按属性信息选择要素
from osgeo import ogr
import os
inShp = 'C:\\Users\\Administrator\\Desktop\\gdal\\样例数据\\shp\\county_popu.shp'
ds = ogr.Open(inShp)
layer = ds.GetLayer()
featureCout = layer.GetFeatureCount()
print(featureCout)
#按照属性条件选择
layer.SetAttributeFilter('Area<200')
layerSlectCount = layer.GetFeatureCount()
print(layerSlectCount)
#把选中的要素输出
outShp = 'C:\\Users\\Administrator\\Desktop\\gdal\\样例数据\\select\\select.shp'
driver = ogr.GetDriverByName('ESRI Shapefile')
if os.access(outShp,os.F_OK):
driver.DeleteDataSource(outShp)
dsOut = driver.CreateDataSource(outShp)
layerOut = dsOut.CreateLayer('test',None,ogr.wkbPolygon)
#获得要素
feature = layer.GetNextFeature() # GetNextFeature每次只获取1个要素
while feature is not None:
layerOut.CreateFeature(feature)#只能拷贝几何位置，属性信息还没有复制了
feature = layer.GetNextFeature()
#选择会对后续操作产生影响
layerSlectCount = layer.GetFeatureCount()
print(layerSlectCount) # 上面已经选择了要素，这里只显示选择了的要素
ds = ogr.Open(inShp) # 如果要获取全部要素，不被前面的选择影响，从新再打开1次shp文件，或者可以直接清空上面的选择就行
# layer.SetSpatialFilter(None) #清空上边的选择也可以
layer2 = ds.GetLayer()
featureCout2 = layer2.GetFeatureCount()
print(featureCout2)
dsOut.Destroy()
4.5.2 按空间位置或sql语句选择要素
import os
from osgeo import ogr
#打开已有矢量数据
inshp = 'C:\\Users\\Administrator\\Desktop\\gdal\\样例数据\\select\\GSHHS_c.shp'
covershp = 'C:\\Users\\Administrator\\Desktop\\gdal\\样例数据\\select\\cover.shp'
ds_in = ogr.Open(inshp)
layer_in = ds_in.GetLayer()
fcount= layer_in.GetFeatureCount()
print(fcount)
ds_cover = ogr.Open(covershp)
layer_cover = ds_cover.GetLayer()
feature_cover = layer_cover.GetNextFeature()
cover = feature_cover.GetGeometryRef()
#cover = feature_cover.geometry()
#根据空间位置选择
layer_in.SetSpatialFilter(cover)
fcount = layer_in.GetFeatureCount()
print(fcount)
#清空上边的选择
layer_in.SetSpatialFilter(None)
#按照矩形坐标*(minx,miny,maxx,maxy)
layer_in.SetSpatialFilterRect()
fcount = layer_in.GetFeatureCount()
print(fcount)
#SQL查询
layer_in.SetSpatialFilterRect(None) # 清除上面矩形选择
fcount = layer_in.GetFeatureCount()
print(fcount)
layername = layer_in.GetName()
result = ds_in.ExecteSQL('select * from {layer_in} where area < 800000'
.format(layer_in=layername))
fcount = result.GetFeatureCount()
print(fcount)
ds_in.ReleaseResultSet(result)
4.6 栅格矢量化，并将线段进行平滑操作
from osgeo import gdal,ogr,osr
import os,cv2
import numpy as np
from skimage import measure
imagePath = 'C:\\Users\\DELL\\Desktop\\Code\\data\\raster2shp.tif'
outShp = 'C:\\Users\\DELL\\Desktop\\Code\\data\\raster2shp1.shp'
image = gdal.Open(imagePath)
if image is None:
print('数据不能正常打开！！！')
bandCount = image.RasterCount
if bandCount != 1:
print('输入数据必须是单波段数据！！！')
#获得影像的投影信息
imgproj = image.GetProjection()
#定义矢量数据
driver = ogr.GetDriverByName('ESRI Shapefile')
shp = driver.CreateDataSource(outShp)
proj = osr.SpatialReference()
proj.ImportFromWkt(imgproj)
layer = shp.CreateLayer('shp',proj,ogr.wkbPolygon)
field = ogr.FieldDefn('id',ogr.OFTString)
shpField = layer.CreateField(field)
#栅格矢量化
prog_func = gdal.TermProgress_nocb
band = image.GetRasterBand(1)
result = gdal.Polygonize(band,None,layer,shpField,[],prog_func)
featCount = layer.GetFeatureCount()
print(featCount)
#创建一个新的shp
outShp2 = 'C:\\Users\\DELL\\Desktop\\Code\\data\\raster2shp2.shp'
shp2 = driver.CreateDataSource(outShp2)
layer2 = shp2.CreateLayer('name',None,ogr.wkbPolygon)
featDefn2 = layer2.GetLayerDefn()
feature2 = ogr.Feature(featDefn2)
# 下面使用3种方法将矢量化后的线段圆滑化
# 第一种方法直接简化
# feature = layer.GetNextFeature()
# while feature is not None:
# geo = feature.GetGeometryRef()
# geom = geo.SimplifyPreserveTopology(12)
# feature2.SetGeometry(geom)
# layer2.CreateFeature(feature2)
# feature = layer.GetNextFeature()
# shp.Destroy() # 写到磁盘里
#第二种方法借助skimage
# width = image.RasterXSize
# height = image.RasterYSize
# gdalImage = image.ReadAsArray(0,0,width,height)
# imageGeoTransform = image.GetGeoTransform()
# contours, hierarchy = cv2.findContours(
# gdalImage, cv2.RETR_TREE, cv2.CHAIN_APPROX_TC89_KCOS)
#
# polygons = ogr.Geometry(ogr.wkbMultiPolygon)
# for i in range(len(contours)):
# data = np.squeeze(contours[i])
# print(len(data.shape))
# if len(data.shape) == 1:
# continue
# contours2 = measure.subdivide_polygon(data)
# lineRing = ogr.Geometry(ogr.wkbLinearRing)
# for contour in contours2:
# x_col = float(contour[1])
# y_row = float(contour[0])
# # 把图像行列坐标转化成地理坐标
# # x_geo = 左上角x + xpixel*x分辨率 + ypixel *旋转角
# row_geo = imageGeoTransform[0] + imageGeoTransform[1] * y_row + x_col * imageGeoTransform[2]
# col_geo = imageGeoTransform[3] + imageGeoTransform[4] * y_row + x_col * imageGeoTransform[5]
# lineRing.AddPoint(row_geo, col_geo)
#
# lineRing.CloseRings()
# polygon = ogr.Geometry(ogr.wkbPolygon)
# polygon.AddGeometry(lineRing)
# polygons.AddGeometry(polygon)
# feature2.SetGeometry(polygons)
# layer2.CreateFeature(feature2)
# shp2.Destroy()
#第三种方法，借助OPencv
#用Opencv中的方法寻找角点并进行简化，然后把角点用gdal构建矢量面
width = image.RasterXSize
height = image.RasterYSize
gdalImage = image.ReadAsArray(0,0,width,height)
imageGeoTransform = image.GetGeoTransform()
contours, hierarchy = cv2.findContours(
gdalImage, cv2.RETR_TREE, cv2.CHAIN_APPROX_TC89_KCOS)
polygons = ogr.Geometry(ogr.wkbMultiPolygon)
for contour in contours:
lineRing = ogr.Geometry(ogr.wkbLinearRing)
for point in contour:
x_col = float(point[0, 1])
y_row = float(point[0, 0])
# 把图像行列坐标转化成地理坐标
# x_geo = 左上角x + xpixel*x分辨率 + ypixel *旋转角
row_geo = imageGeoTransform[0] + imageGeoTransform[1] * y_row + x_col * imageGeoTransform[2]
col_geo = imageGeoTransform[3] + imageGeoTransform[4] * y_row + x_col * imageGeoTransform[5]
lineRing.AddPoint(row_geo, col_geo)
lineRing.CloseRings()
polygon = ogr.Geometry(ogr.wkbPolygon)
polygon.AddGeometry(lineRing)
polygons.AddGeometry(polygon)
feature2.SetGeometry(polygons)
layer2.CreateFeature(feature2)
shp2.Destroy()
5 一些处理工具
GDAL已经将好多常用的工具封装好了，主要包括两类:
（1）.exe可执行程序，我们用python代码直接调用就行；
（2）.python文件，这里有非常多的库和功能，我们直接用，可以直接放到我们的代码中。
Gdal官方文档中对这些工具都有详细介绍说明。
.exe文件在安装包的下面目录中
.py文件在安装包的下面目录中
6 总结（简单，半天就学会了）
gdal非常简答，就包括3部分，org矢量处理,gdal栅格处理，常用的1些工具集。
记住处理的都是一些类，栅格类，数量数据集类，我们主要用这些类的方法操作它的属性。
一共就分为三部分：
第一部分是gdal处理栅格影像：
- 读取栅格，看3.3节
- 写出栅格，看3.4节
- 计算栅格，看3.4.4，就是numpy的运算
第二部分就是ORG处理矢量数据：
1.读取shp文件，就直接创建驱动读取，就可以读取数据集，图层，要素的属性信息；
2.创建点、线、面文件，看上面4.2，4.3，4.4节的代码，创建文件有两种方式，1是直接创建几何要素，2是使用wkt字符串创建几何要素。
3.对要素的操作，选择要素，栅格矢量化，矢量线段圆滑化等。
第三部分就是一些工具集：
主要就是.exe和.py工具集。
参考资料
[1] python与GDAL-空间数据处理入门教程_Python视频-51CTO学堂
51cto官网中的python与GDAL-空间数据处理入门教程的课程代码和资料，个人仓库
https://github.com/GisXiaoMa/gdal_51cto
[2] Gdal官方说明文档: https://gdal.org/
[3] 犹他州立大学的教程:
https://www.osgeo.cn/python gdal utah tutorial/index.html
[4] OSGEO中国官网:https://www.osgeo.cn/
[5] 卜坤的python与开源GIS
[6] 李民录 GDAL源码剖析书籍
文章出处登录后可见！
已经登录？立即刷新
