网站建设安全标准,新昌网站制作,校园文化创意产品设计,天津模板建站代理从GDAL提供的实用程序来看#xff0c;很多程序的后缀都是 .py #xff0c;这充分地说明了Python语言在GDAL的开发中得到了广泛的应用。
1. 打开已有的GeoTIF文件
下面我们试着读取一个GeoTiff文件的信息。第一步就是打开一个数据集。 from osgeo import gdal…从GDAL提供的实用程序来看很多程序的后缀都是 .py 这充分地说明了Python语言在GDAL的开发中得到了广泛的应用。
1. 打开已有的GeoTIF文件
下面我们试着读取一个GeoTiff文件的信息。第一步就是打开一个数据集。 from osgeo import gdaldataset gdal.Open(/gdata/geotiff_file.tif)既然已经将一个GeoTIFF文件打开为一个GDAL可操作的对象 下面来看一下都能对其进行怎样的操作。
Python提供了 dir() 内省函数 可以快速查看一下当前对象可用的操作 dir() 函数可能是 Python 自省机制中最著名的部分了。它可以返回传递给它的任何对象的属性名称经过排序的列表。 如果不指定对象则 dir() 返回当前作用域中的名称。 dir(dataset)[:3] [... ...] dir(dataset)[-3:][AbortSQL,AddBand,AddFieldDomain,... ...,__weakref__,this,thisown]下面看一下如何获取文件的一些基本信息需要用到下面的一些函数与属性。 dataset.GetDescription() 获得栅格的描述信息 dataset.RasterCount 获得栅格数据集的波段数 dataset.RasterXSize 栅格数据的宽度(X方向上的像素个数) dataset.RasterYSize 栅格数据的高度(Y方向上的像素个数) dataset.GetGeoTransform() 栅格数据的六参数。 GetProjection() 栅格数据的投影
通过下面小节我们逐一解释一下。
2. 读取影像的元数据
从元数据的作用来看它更多地是为工程服务的。 客观地说GDAL对元数据的支持并不好 它并没有直接的元数据处理接口 更没有实现元数据的相关标准。 但是GDAL已经提供了足够方便的函数可以读取影像的一些信息 从而方便对数据进行处理。 GDAL一般是以字典的形式对元数据进行组织的 但是对于不同的栅格数据类型元数据的类型与键值可能都不一样。 目前 国际上对空间元数据标准内容进行研究的组织主要有三个分别是欧洲标准化委员会(CEN/TC 287)、 美国联邦地理数据委员会(FGDC)和国际标准化组织地理信息/地球信息技术委员会(ISO/TC 211)。 如果要进行元数据处理可以考虑将元数据信息写入到XML文件中。 这个问题扩展开来就不是本书关心的内容的在此就不再多说了。
我们先来看一下最常用的GeoTIFF文件的元数据信息。 GDAL可以作为数据集级别的元数据来处理下面的基本的TIFF标志。 TIFFTAG_DOCUMENTNAME TIFFTAG_IMAGEDESCRIPTION TIFFTAG_SOFTWARE TIFFTAG_DATETIME TIFFTAG_ARTIST TIFFTAG_HOSTCOMPUTER TIFFTAG_COPYRIGHT TIFFTAG_XRESOLUTION TIFFTAG_YRESOLUTION TIFFTAG_RESOLUTIONUNIT TIFFTAG_MINSAMPLEVALUE (read only) TIFFTAG_MAXSAMPLEVALUE (read only)
使用Python 打开数据来读取一下数据的信息 from osgeo import gdaldataset gdal.Open(/gdata/geotiff_file.tif)dataset.GetMetadata(){AREA_OR_POINT: Area, PyramidResamplingType: NEAREST}dir(dataset)[:10][AbortSQL,AddBand,AddFieldDomain,AddRelationship,AdviseRead,BeginAsyncReader,BuildOverviews,ClearStatistics,CommitTransaction,CopyLayer]GetMetadata() 方法可以访问数据的元数据信息元数据信息对于每个数据都是不一样的。 比如再打开另外一个文件 ds gdal.Open(/gdata/lu75c.tif)ds.GetMetadata(){AREA_OR_POINT: Area,TIFFTAG_XRESOLUTION: 1,TIFFTAG_YRESOLUTION: 1}这个文件只对两个TIFF标志进行了定义还有一个并不是TIFF标志定义的。
GetDescription() 获得栅格的描述信息 dataset.GetDescription()/gdata/geotiff_file.tif看来这里的图像描述是图像的路径名 但是这是和各种不同数据集相关的 不同数据集有不同的描述。
获取栅格数目
栅格数据集是由多个数据构成的在GDAL中每一个波段都是一个数据集 不仅如此栅格数据集还可能包含有子数据集每子数据集又可能包含有波段。
先来看一下刚才打开的数据的RasterCount dataset.RasterCount3这是一个由7个波段构成的Landsat遥感影像。
再看一个MODIS L1B数据 mds gdal.Open(/gdata/MOD09A1.A2009193.h28v06.005.2009203125525.hdf)mds.RasterCount0运行结果居然是 0 。这意味着当前的数据集 dataset 中的栅格数目是 0 。 实际上MODISL1B的数据格式是HDF格式的它的数据是以子数据集组织的要获取其相关的数据的信息需要继续访问其子数据集。
影像大小
栅格数据的大小指出了影像以像元为单位的宽度与高度。 dataset.RasterXSize,dataset.RasterYSize(1500, 900)可以看出我们的影像大小。
获得空间参考
下面看一下如果从栅格数据集中获取其投影与空间参考信息。 更多的关于投影与空间参考的讨论会在后面章节介绍。
GetGeoTransform() 地理仿射变换参数。
对于遥感影像来说它需要在地理空间中进行定位。 在GDAL中这有两种方式其中一种是使用六个参数坐标转换模型。 这个模型的具体实现在不同的软件中是不一样的。 在GDAL中这六个参数包括 左上角坐标 像元X、Y方向大小 旋转 等信息。 要注意 Y 方向的像元大小为负值。 ds.GetGeoTransform()(1852951.7603168152, 30.0, 0.0, 5309350.360150607, 0.0, -30.0)获得投影信息
使用 GetProjection() 函数可以比较容易地获取数据集的投影信息。 但是对于什么是地图投影以及如何在GDAL中实现就不是这么容易了。 此处我们只是简单地运行一下更详细的解释会在后面章节中展开。 ds.GetProjection()PROJCS[unnamed,GEOGCS[unknown,DATUM[unnamed,SPHEROID[unretrievable - using WGS84,6378137,298.257223563]],PRIMEM[Greenwich,0],UNIT[degree,0.0174532925199433,AUTHORITY[EPSG,9122]]],PROJECTION[Albers_Conic_Equal_Area],PARAMETER[latitude_of_center,0],PARAMETER[longitude_of_center,105],PARAMETER[standard_parallel_1,25],PARAMETER[standard_parallel_2,47],PARAMETER[false_easting,0],PARAMETER[false_northing,0],UNIT[metre,1,AUTHORITY[EPSG,9001]],AXIS[Easting,EAST],AXIS[Northing,NORTH]]3. 使用GDAL获取栅格数据波段信息
上面我们介绍了针对数据集操作的主要函数。 但是如果需要了解栅格数据的更多信息我们就需要看一下遥感图像处理中更常用到的波段操作的函数。
获取数据集的波段
GetRasterBand() 函数可以获得栅格数据集的波段。这是函数的参数使用波段的索引值。 from osgeo import gdaldataset gdal.Open(/gdata/lu75c.tif)dataset.RasterCount1band dataset.GetRasterBand(1)这里我们获取了第一个波段(红色值组成的表)。
注意这里的波段获取和一般编程语言中数组下标的索引方式不一样开始是 1 不是 0 。
查看波段的基本信息
下面我们现在来看看刚才读取出来的 band 有些什么东西可以供我们操作。
与操作数据集一样GDAL同样提供了查看波段基本信息的函数。 band dataset.GetRasterBand(1)dir(band)[:6] [... ...] dir(band)[-3:][AdviseRead,AsMDArray,Checksum,ComputeBandStats,ComputeRasterMinMax,ComputeStatistics,... ...,__weakref__,this,thisown]看一下常用的操作。这些也是用来获取波段的属性信息。
获取波段大小 band.XSize, band.YSize, band.DataType(6122, 4669, 3)执行以上代码得到了波段图像的宽和高像元为单位。 对于我们所使用的影像 这个与 dataset 中使用 RasterXSize() 与 RasterYSize() 获取的值一致。 DataType 是图像中实际数值的数据类型表示 8 位无符整型。
获取波段数据的属性 band.GetNoDataValue()band.GetMaximum()print(band.GetMinimum())Noneband.ComputeRasterMinMax()(-1.0, 66.0)Maximum 是表示在本波段数值中最大的值当然 Minimum 就是表示本波段中最小的值。 通过运行结果我们可以看到在一开始RasterXSize()和RasterYSize()都没有值。 因为对于文件格式不会有固有的最大最小值。 所以我们可以通过函数 ComputeRasterMinMax() 计算得到。 注意这里的最大最小值不包括“无意义值” 也就是上面显示的 NoDataValue 。
4. 其他数据格式格式
使用GDAL读取ENVI数据格式
ENVI栅格文件格式
ENVI使用的是通用栅格数据格式包含一个简单的二进制文件a simple flat binary和一个相关的ASCII文本的头文件。这也保证了单个ENVI栅格文件没有大小上限。
1头文件
ENVI头文件包含用于读取图像数据文件的信息它通常创建于一个数据文件第一次被ENVI读取时。单独的ENVI头文本文件提供关于图像尺寸、 嵌入的头文件若存在、数据格式及其它相关信息。所需信息通过交互式输入或自动地用“文件吸取”创建并且以后可以编辑修改。 使用者可以在ENVI之外使用一个文本编辑器生成一个ENVI头文件。
2数据文件
通用栅格数据都会存储为二进制的字节流通常它将以BSQband sequential按波段顺序储存、BIPband interleaved by pixel按波段像元交叉储存或者BILband interleaved by line按波段行交叉储存的方式进行存储。
ENVI栅格文件必须包含着两个文件其中头文件的后缀名为.hdr数据文件的后缀随意甚至可以不带后缀名。 这两个文件是通过文件名来关联即数据文件和头文件名称一致。
GDAL读取HDF数据格式
由于modis卫星数据跟我们经常遇到的geotif数据组织方式不一样读取的时候一定要特别注意。geotif数据一般是一个文件包含了多个波段的数据 而modis呢一个文件包含了多各SUBDATASETSGDAL 每个SUBDATASETS又包含多个波段数据。 另外默认编译的GDAL并不包含对MODIS数据支持 需要单独下载针对HDF4HDF5的源码再修改下make.opt文件 这时再编译GDAL就支持modis数据的读写了。
