python 与GDAL.pdf-资料库

这段时间在学习 python 与 GDAL，网上这方面比较经典的有犹他大学开源 GIS 公开课，总共有 7 周的课时。与自己大学时老师照本念经的灌输教育，国外大学的教育似乎更偏向于学生的社会实践能力。 GDAL/OGR 官方文档的参考也比较有学习意义。我建议各位有缘人充分利用 python 的内置函数“dir”，通过该函数你可以看到该对象所有的方法与属性。我就是通过 dir 这个函数一个个去学习对象。我总是忘记 GDAL 中文件句柄的方法，我就 dir 以下就显示出来文件句柄对象所有的方法与属性。或许，你会感到惊讶，python 不就是脚本语言吗，怎么会与对象挂钩。哈哈，没错， python 是高级的面向对象语言。它就是对象为基础的高级语言。Python 功能十分强大， python 的世界也十分多姿多彩，希望你可以尽情的在这个世界刷吧。国内关于这方面中文参考资料多，但不系统。我花了大概两周的时间整理了这些资料。如果你静下心来学习，你会有以下收获： 1、最基本的矢量数据读取，输出 2、最基本的栅格数据读写、分块读写 3、 GDAL、ORG 强大的工具文档中实例代码为本人自己亲写、亲测，代码功能简单。但万事都是从简单开始。相信各位有缘人一定会比本人强，在遥感、GIS 等领域有好的建树。祝愿有缘人身体健康。没有女朋友的男票们早日找到女友。生活不易，若你感觉文档内容不错，可以给本人一个打赏。

使用开源 GIS 的原因优点：  个人或小规模公司经济实惠  开发工具非常有用，bug 快速确定修复  可以开发 windows 之外的程序  师妹会觉得你特别牛逼，后面省略无数字缺点：  没有内置的地理处理内核  用户群较少 OGR 矢量文件读写 OGR 核心类结构如下： OGR 矢量文件读取： 1. 导入 ORG 模块 2. 获取矢量文件驱动 3. 打开矢量数据，获取文件句柄 4. 基于文件句柄，获取图层，及图层范围 5. 基于图层，获取特征要素，特征要素的属性 6. 获取特征要素的几何体形状及投影参考 7. 特征要素销毁 8. 文件句柄销毁 OGRGeomeryOGRCurveOGRGeometryCollectionOGRPointOGRSurfaceOGRLineStringOGRMultiLineStringOGRMultiPointOGRMultiPolygonOGRPolygonOGRLinearRing

代码如下： import os,sys #导入 ogr 库 import ogr,osr if __name__ == '__main__': shpFile = r'F:\湖北省区县界\湖北省区县界\湖北省区县界.shp' #1:获取文件驱动，打开文件 driver = ogr.GetDriverByName('ESRI Shapefile') shpDs = driver.Open(shpFile,0) #0 表示以只读形式打开矢量文件，1 表示以可写的形式打开矢量文件 #2:矢量文件图层获取 shpLayer = shpDs.GetLayer(0) #图层获取,对于 shp 文件数据格式默认 0 为第 1 层 layerCount = shpDs.GetLayerCount()#获取矢量图层个数 layerExtent = shpLayer.GetExtent()#获取矢量图层四个角点 fieldCount = shpLayer.GetFieldCount()#获取图层属性名字 field = shpLayer.GetField(0)#获取图层属性第 0 个字段 #3:矢量图层中获取矢量要素 featureList = [] print(dir(shpLayer)) featureCount = shpLayer.GetFeatureCount() for i in range(featureCount): feature = shpLayer.GetFeature(i) #图层要素获取 geometryRef =feature.GetGeometryRef()#图层要素形状获取 featureList.append(feature) feature.Destroy()#图层要素销毁 #4:获取图层的空间参考 spatialRef = shpLayer.GetSpatialRef() #矢量文件句柄销毁 shpDs.Destroy()

OGR 矢量文件写入： 1. 导入 ORG 模块 2. 获取矢量文件驱动 3. 创建矢量数据，获取文件句柄 4. 基于文件句柄，创建图层及图层属性 5. 基于图层，获取特征要素 6. 创建图层要素，要素属性 7. 创建要素几何体形状，投影 8. 特征要素销毁 9. 文件句柄销毁代码如下：课程中的实例仅仅是面向教学，下面实例经过修改可以直接读取一个矢量文件，并将该矢量文件输出。如果你感兴趣，可以增加其他功能 import os,sys #导入 ogr 库 import ogr,osr if __name__ == '__main__': shpFile = r'F:\湖北省区县界\湖北省区县界\湖北省区县界.shp' #1:获取文件驱动，打开文件，获取文件驱动 driver = ogr.GetDriverByName('ESRI Shapefile') shpDs = driver.Open(shpFile,0) #0 表示以只读形式打开矢量文件，1 表示以可写的形式打开矢量文件 shpLayer = shpDs.GetLayer(0) #图层获取,对于 shp 文件数据格式默认 0 为第 1 层 orgSpatialRef = shpLayer.GetSpatialRef()#获取图层参考系 featureCount = shpLayer.GetFeatureCount()#获取图层要素个数 outShp = r'F:\11.shp' outDs = driver.CreateDataSource(outShp) #通过图层复制的方式将源图层字段定义及图层要素定义复制到目标文件中 outLayer = outDs.CopyLayer(shpLayer,'11') featureDefn = outLayer.GetLayerDefn() #将目标图层的参考系赋值

for i in range(featureCount): feature = shpLayer.GetFeature(i) #获取文件属性要素 outFeature = ogr.Feature(featureDefn) #获取输出要素特征定义 outFeature.SetGeometry(feature.GetGeometryRef()) #特征要素几何参考定义 outLayer.CreateFeature(outFeature)#创建图层要素 feature.Destroy() #输入特征销毁 outFeature.Destroy()#输出特征销毁 shpDs.Destroy() #源文件句柄销毁 outDs.Destroy() #目标文件句柄销毁 OGR 矢量文件几何体创建 1. point 几何体创建 2. line 几何体创建 3. polygon 几何体创建 4. linering 几何体创建 5. multipoint 几何体创建 6. multipolygone 几何体创建 OGR 矢量文件分析 OGR 矢量文件要素选择： 1. 根据属性选择 SetAttributeFilter() 2. 根据空间范围选择 SetSpatialFilter() SetSpatialFilterRect(, , , ) 3. 根据 SQL 语句选择 ExecuteSQL() OGR 矢量文件要素判断： 1. 面状要素求交

poly2.Intersect(poly1) 2. 面状几何要素相离 poly2.Disjoint(poly1) 3. 线面几何要素相切 poly2.Touches(poly1) 4. 线面几何要素相交 poly2.Crosses(poly1) 5. 几何要素在里面 poly3.Within(poly2) 6. 几何要素包含 poly3.Contains(poly2) 7. 几何要素重叠 poly2.Overlaps(poly1) OGR 矢量文件要素操作： 1. 几何要素交集 poly3.Intersection(poly2) 2. 几何要素并集 poly3.Union(poly2) 3. 几何要素差集 poly3.Difference(poly2) 4. 几何要素非重叠部分合并 poly3.SymmetricDifference(poly2) 5. 几何要素缓冲区建立 .Buffer() 6. 几何要素相等判断 .Equal() 7. 几何要素最短距离分析 .Distance() 8. 几何要素外包范围 .GetEnvelope()

GDAL 栅格文件读取 GDAL 核心类如下图所示： GDAL 文件读取： 1. 导入 GDAL 模块 2. 获取栅格文件驱动 3. 打开栅格文件，获取文件句柄 4. 基于文件句柄，获取文件波段 5. 基于波段读取影像数据 6. 基于数据集获取影像地理六参数 7. 基于投影获取影像波段个数 8. 波段销毁 9. 文件句柄销毁代码如下： from osgeo import gdal #导入 gdal 库 gdal.AllRegister() #获取文件类型驱动 Driver = gdal.GetDriverByName('GTiff') #打开 tiff 数据 fh = gdal.Open(r'path') #获取文件属性 bands = fh.RasterCount #文件波段数 xsize = fh.RasterXSize #文件行数 ysize = ysize = fh.RasterYSize #文件列数 band_type = fh.GetRasterBand(1).DataType #数据类型 projection = fh.GetProjection() #投影参数信息 geotransform = fh.GetGeoTransform() #文件几何信息 GDALMajorObjectGDALDatasetGDALDriverGDALDriverManagerGDALRasterBand

ulx = geotransform[0] #左上 x 坐标 resolutionx = geotransform[1] #x 方向分辨率 rotationx = geotransform[2] #x 方向旋转 uly = geotransform[3] #左上 y 坐标 rotationy = geotransform[4] #y 方向旋转 resolutiony = geotransform[5] #y 方向分辨率 band = fh.GetRasterBand(1) #获取波段 ct = band.GetRasterColorTable() #数据颜色表 data = band.ReadAsArray() #获取影像数据 band = None #波段数据清除 fh = None #解除程序对文件的占用影像分块读取方式： 1. 影像逐行读取 for i in range(rows): data = band.ReadAsArray(0, i, cols, 1) 确保按 Y 方向能够读取完数据，若数据超出则报错 2. 影像按多行读取 Bsize = n for i in range(0, rows, bSize): if i + bSize < rows: size = bSize else: size = rows -i data = band.ReadAsArray(0, i, cols, size) 3. 影像按行*列块逐块读取 xBSize = m yBSize = n for i in range(0, rows, yBSize): if i + yBSize < rows: numRows = yBSize

资料库

python 与GDAL.pdf

相关推荐

开发技术

热门标签

最新资料