GDALpython教程（1）——⽤OGR读写⽮量数据本教程的讲义和源码都是取⾃Utah State University的openGIS课程相关资料，包括讲义、源码、数据样例，请从此处下载本⼈只是做点翻译，写写学习体会⽽已，版权属于原作者。
欢迎转载，不过别忘了上⾯这段话。
==================================================为什么⽤open source？优点1. 免费，适合个⼈和⼩公司2. 强⼤的开发⼯具，找bug更容易3. 跨平台，windows和linux都能⽤4. 拉风！缺点1. 没有内嵌地理处理器2. ⽤的⼈少Open source RS/GIS模块1. OGR⽮量库：简单的⽮量数据读写，是GDAL的⼀部分2. GDAL地理空间数据抽象库：a) 读写栅格数据b) ArcGIS也是基于GDAL开发的c) C++库，但是可以⽤python调⽤相关模块1. Numeric：⾼速的数组处理，对栅格数据尤其重要2. NumPy：下⼀代的Numeric3. 更强⼤的gis库导⼊库：import ogr或者：from osgeo import ogr万能的⽅法是：try:from osgeo import ogrexcept:import ogr要读取某种类型的数据，必须要先载⼊数据驱动，也就是初始化⼀个对象，让它“知道”某种数据结构。
import ogrdriver = ogr.GetDriverByName(‘ESRI Shapefile’)数据驱动driver的open()⽅法返回⼀个数据源对象open(<filename>, <update>)其中update为0是只读，为1是可写例如：from osgeo import ogrdriver = ogr.GetDriverByName('ESRI Shapefile')filename = 'C:/Users/gongwei/Documents/My eBooks/python_and_sage/GDAL python/test/ospy_data1/sites.shp' dataSource = driver.Open(filename,0)if dataSource is None:print 'could not open'sys.exit(1)print 'done!'注意filename⼀定要写绝对路径！因为⼀定要⽤绝对路径，为了简化代码，经常会使⽤到os.chdir()读取数据层layer = dataSource.GetLayer(0)⼀般ESRI的shapefile都是填0的，如果不填的话默认也是0.再看看这个数据层⾥⾯有⼏个点呢？n = layer.GetFeatureCount()print 'feature count:', n读出上下左右边界extent = layer.GetExtent()print 'extent:', extentprint 'ul:', extent[0], extent[3]print 'lr:', extent[1], extent[2]读取某⼀要素feature（总算切⼊正题了），这⾥读取的是⼀个点feat = layer.GetFeature(41)fid = feat.GetField('id')print fidfeat = layer.GetFeature(0)fid = feat.GetField('id') #should be a different idprint fid另外还有按顺序读取feature，循环遍历所有的featurefeat = layer.GetNextFeature() #读取下⼀个while feat:feat = layer.GetNextFeature()later.ResetReading() #复位提取feature的⼏何形状geom = feat.GetGeometryRef()geom.GetX()geom.GetY()print geom.释放内存feature.Destroy()关闭数据源，相当于⽂件系统操作中的关闭⽂件dataSource.Destroy()读完了再说怎么写创建新⽂件driver.CreateDataSource(<filename>)但是这个⽂件不能已经存在了，否则会出错创建新的layerdataSource.CreateLayer(<name>,CreateLayer(<name>, geom_type=<OGRwkbGeometryType>, [srs])举个例⼦：ds2 = driver.CreateDataSource('test.shp')layer2 = ds2.CreateLayer('test', geom_type=ogr.wkbPoint)要删除⼀个shp⽂件driver.DeleteDataSource('test.shp')要添加⼀个新字段，只能在layer⾥⾯加，⽽且还不能有数据添加的字段如果是字符串，还要设定宽度fieldDefn = ogr.FieldDefn('id', ogr.OFTString)fieldDefn.SetWidth(4)layer.CreateField(fieldDefn)添加⼀个新的feature，⾸先得完成上⼀步，把字段field都添加齐了然后从layer中读取相应的feature类型，并创建featurefeatureDefn = layer.GetLayerDefn()feature = ogr.Feature(featureDefn)设定⼏何形状feature.SetGeometry(point)设定某字段的数值feature.SetField('id', 23)将feature写⼊layer layer.CreateFeature(feature)。
GDAL的基本操作 上⼀节简单介绍了GDAL，这⼀节将介绍⼀些GDAL的基本操作，如影像读写、波段提取、波段合成等。
代码均⽤python编写。
1.遥感影像的读写 众所周知，遥感影像是以栅格形式存储的，GDAL中使⽤dataset表⽰⼀个栅格数据（使⽤抽象类表⽰），⼀个dataset包含了对于栅格数据的波段，空间参考以及元数据等信息。
⼀张GeoTIFF遥感影像，⼀张DEM影像，或者⼀张⼟地利⽤图，在GDAL中都是⼀个GDALDataset。
坐标系统（使⽤OGC WKT格式表⽰的空间坐标系统或者投影系统）地理放射变换（使⽤放射变换表⽰图上坐标和地理坐标的关系）GCPs（⼤地控制点记录了图上点及其⼤地坐标的关系，通过多个⼤地控制点可以重建图上坐标和地理坐标的关系）元数据（键值对的集合，⽤于记录和影像相关的元数据信息）栅格波段（使⽤GDALRasterBand类表⽰，真正⽤于存储影像栅格值，⼀个栅格数据可以有多个波段）颜⾊表（Color Table⽤于图像显⽰） GDAL使⽤Open（）函数读取影像数据，函数返回Dataset对象。
GDAL⽬前⽀持约100种格式的栅格数据读取，包括ERDAS Imagine、ENVI、GRASS、GeoTIFF、HDF4、HDF5、TIFF、JPEG、JPEG2000、PNG、GIF、BMP等。
import numpy as npfrom osgeo import gdalin_ds = gdal.Open(r"D:\data\test.tif") # 打开样本⽂件xsize = in_ds.RasterXSize # 获取⾏列数ysize = in_ds.RasterYSizebands = in_ds.RasterCount # 获取波段数geotransform = in_ds.GetGeoTransform() # 获取投影信息projection = in_ds.GetProjectionRef()block_data = in_ds.ReadAsArray(0,0,xsize,ysize).astype(np.float32)# 获取影像信息B = block_data[0, :, :] #获取第⼀波段数据G = block_data[1,:, :] #获取第⼆波段数据R = block_data[2,:, :] #获取第三波段数据 Dataset对象的RasterYSize、RasterXSize和RasterCount属性分别返回栅格数据的⾏数、列数和波段数。
很好的GDAL初学者资料1、编译加使用：/gisblog/archive/2008/10/06/3021933.aspx如何在CSharp中使用GDAL？收藏GDAL库的简洁、高效深受开发人员的喜爱，很多开源的GIS软件甚至是商业GIS软件都使用了这个库。
GDAL使用C++，在Visual studio环境下开发，对C，C++的支持当然不会有什么问题。
但对于C#、Python、Java来说，过程就相对繁琐一些，对初学者造成了不少的麻烦。
本文以GDAL1.50为例，详细说明在CSharpe语言中使用GDAL的环境配置问题，至于Python的环境配置问题，在另一篇文章“如何在CSharpe中使用GDAL”中另行讨论。
希望对初用GDAL的人员有所帮助。
GDAL从1.4.1开始添加了swig，以支持.NET环境下的C#语言。
在swig 文件夹中包含了C#源代码。
以下是主要的处理步骤：1、准备一个C++编译器。
建议使用Visual studio2005。
2、下载源代码。
地址：/gdal/gdal150.zip。
3、把gdal150.zip解压至“C:\gdal-1.5.0”，当然放在其它文件夹也可以。
4、修改源代码中的错误。
源文件“C:\gdal-1.5.0\frmts\leveller\levellerdataset.cpp”的第171行的内容如下：{ "?, kPI / 180.0, UNITLABEL_DEGREE },在?号后面加上双引呈"，改成如下形式：{ "?", kPI / 180.0, UNITLABEL_DEGREE },你可以使用任何一种文本文件的编辑工具来进行修改。
笔者尝试过下载GDAL1.52的版本，存在同样的问题，真想不到GDAL库也会犯如此低级的错误。
5、创建安装的目标文件夹“C:\gdal-runtime”，稍后的编译配置中，我们将配置GDAL把编译后将会安装到这个文件夹内。
pythongdal读取栅格数据1.gdal包简介gdal是空间数据处理的开源包，其⽀持超过100种栅格数据类型，涵盖所有主流GIS与RS数据格式，包括Arc/Info ASCII Grid(asc)，GeoTiff (tiff)，Erdas Imagine Images(img)，ASCII DEM(dem) 等格式。
2.安装gdal包（1）通过此链接查找并下载gdal包：（根据⾃⼰的系统及python版本选择对应的gdal包）（2）在dos命令⾏下切换到包所在的⽂件夹，并运⾏pip install 包名进⾏安装。
3.读取栅格数据#导⼊gdal包from osgeo import gdal#导⼊numpy包（⽀持⾼维数组和矩阵运算，也提供了许多数组和矩阵运算的函数）import numpy as np#打开⽂件dataset=gdal.Open("fdem.tif")#栅格矩阵的列数im_width = dataset.RasterXSize#栅格矩阵的⾏数im_height = dataset.RasterYSize#波段数im_bands = dataset.RasterCount#仿射矩阵，左上⾓像素的⼤地坐标和像素分辨率。
#共有六个参数，分表代表左上⾓x坐标；东西⽅向上图像的分辨率；如果北边朝上，地图的旋转⾓度，0表⽰图像的⾏与x轴平⾏；左上⾓y 坐标；#如果北边朝上，地图的旋转⾓度，0表⽰图像的列与y轴平⾏；南北⽅向上地图的分辨率。
im_geotrans = dataset.GetGeoTransform()#地图投影信息im_proj = dataset.GetProjection()#读取某⼀像素点的值#（1）读取⼀个波段，其参数为波段的索引号，波段索引号从1开始(我打开的这幅图像只有⼀个波段)band=dataset.GetRasterBand(1)#（2）⽤ReadAsArray(<xoff>, <yoff>, <xsize>, <ysize>)，读出从(xoff,yoff)开始，⼤⼩为(xsize,ysize)的矩阵。
