gdal统计指定区域内的栅格值（gdal概览）
导读：1 gdal库 2 栅格驱动...
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 = None4 矢量数据处理（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源码剖析书籍
声明：本站所有文章 ，如无特殊说明或标注 ，均为本站原创发布 。任何个人或组织 ，在未征得本站同意时 ，禁止复制 、盗用 、采集 、发布本站内容到任何网站 、书籍等各类媒体平台 。如若本站内容侵犯了原著者的合法权益 ，可联系我们进行处理 。
创心域SEO版权声明：以上内容作者已申请原创保护，未经允许不得转载，侵权必究！授权事宜、对本内容有异议或投诉，敬请联系网站管理员，我们将尽快回复您，谢谢合作！