gdal概览

大家好，欢迎来到IT知识分享网。

1 gdal库

2 栅格驱动

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方法读取坐标，简单介绍一下该方法，该方法会返回以下六个参数

3.1.2 tif文件的地理坐标（两种情况）

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 读取栅格像元

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 创建栅格影像

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 创建点要素

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 < ' .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 一些处理工具

6 总结（简单，半天就学会了）

1. 读取栅格，看3.3节
2. 写出栅格，看3.4节
3. 计算栅格，看3.4.4，就是numpy的运算

参考资料