gdal概览(gdal官方文档)

gdal概览 GDAL1 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库

推荐整理分享gdal概览(gdal官方文档)，希望有所帮助，仅作参考，欢迎阅读内容。

文章相关热门搜索词:gdal使用,gdal_grid,gdal grid,gdal geos,gdal_grid,gdaltransform,gdal_grid,gdal_data,内容如对您有帮助，希望把文章链接给更多的朋友！

2 栅格驱动

读取栅格数据的流程： （1）获取栅格驱动(会有1个默认的，可以不创建，最好创建上，看着清晰)； （2）驱动进行注册（会有默认值，可以不创建，创建的话可以选择只读，或者读写都可以，相当于给栅格驱动设置一个权限）； （3）通过驱动获取数据集； （4）通过数据集操作栅格属性；

try: from osgeo import gdal,ogr # gdal是处理栅格，ogr处理矢量except: import gdaldriverCount = 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 gdalfilePath = '1.tif' # tif文件路径dataset = gdal.Open(filePath) # 打开tifadfGeoTransform = 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,ogrexcept: import gdalgdal.AllRegister()dataset = gdal.Open('C:\\Users\\Administrator\\Desktop\\gdal\\样例数据\\GF1\\GF1.tif')print(dir(dataset))#波段信息bandCount = dataset.RasterCountprint(bandCount)#大小weight,height = dataset.RasterXSize,dataset.RasterYSizeprint(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,ogrexcept: import gdalgdal.AllRegister() # 只读，这里前面没有指定驱动，它自动会用tifdataset = gdal.Open('C:\\Users\\Administrator\\Desktop\\gdal\\样例数据\\GF1\\GF1.tif')width = dataset.RasterXSize # 数据集中的各种属性信息height = dataset.RasterYSizebandCont = dataset.RasterCountfor 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 gdalimport numpy as npsrcFile = 'C:\\Users\\Administrator\\Desktop\\gdal\\样例数据\\GF1\\GF1.tif'dstFile = 'C:\\Users\\Administrator\\Desktop\\gdal\\样例数据\\GF1\\GF1Create.tif'dataSet = gdal.Open(srcFile)width = dataSet.RasterXSizeheight = dataSet.RasterYSizebandCount = dataSet.RasterCountdriver = 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 gdalimport numpy as npsrcFile = '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 gdalimport numpy as npgdal.AllRegister()src = 'C:\\Users\\Administrator\\Desktop\\gdal\\样例数据\\GF1\\GF1.tif'dataset = gdal.Open(src)width = dataset.RasterXSizeheigth = dataset.RasterYSizebandCount = dataset.RasterCountblock = 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 gdalimport numpy as npdef 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,transformdef 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 = Noneif __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 gdalimport numpy as npdataset = gdal.Open('C:\\Users\\Administrator\\Desktop\\gdal\\样例数据\\GF1\\GF1-NIR.tif')width = dataset.RasterXSizeheigth = dataset.RasterYSizebandCount = dataset.RasterCountif 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)#计算NDVINDVI=(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 = None4 矢量数据处理（OGR库）

4.1 读取矢量文件from osgeo import gdalfrom osgeo import ogrimport 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,ogrimport osout_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 ogrimport osout_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 ogrout_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 ogrimport osinShp = '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 osfrom 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,osrimport os,cv2import numpy as npfrom skimage import measureimagePath = '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.RasterCountif 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_nocbband = image.GetRasterBand(1)result = gdal.Polygonize(band,None,layer,shpField,[],prog_func)featCount = layer.GetFeatureCount()print(featCount)#创建一个新的shpoutShp2 = '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.RasterXSizeheight = image.RasterYSizegdalImage = 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源码剖析书籍