如何使用Python中的GDAL库对遥感影像进行读取和存储 |
您所在的位置：网站首页 › python保存jpg文件 › 如何使用Python中的GDAL库对遥感影像进行读取和存储 |
|
与自然影像稍有不同，遥感影像一般有多个波段(波段数大于4)，数据格式多为tiff，比如国产高分一号/二号遥感影像，多光谱影像有4个波段，分别为蓝(1)、绿(2)、红(3)以及近红外(4)波段，真彩色是321显示，假彩色是432显示。 利用遥感或GIS处理软件，如ENVI、ERDAS、ArcGIS等，可以很方便快捷地实现对图像进行读取、处理等操作。 然而，很多时候我们需要对大量的图像数据进行批量处理操作，这时候就需要我们自己编写代码来实现。这篇博客就简单地介绍，如何在Python中实现遥感影像的读取和存储操作。Python中，可以读取和存储遥感影像的库有一些，这里我仅介绍比较常用且经典的库--GDAL。 GDAL的安装百度百科：GDAL(Geospatial Data Abstraction Library)是一个在X/MIT许可协议下的开源栅格空间数据转换库。它利用抽象数据模型来表达所支持的各种文件格式。它还有一系列命令行工具来进行数据转换和处理。GDAL提供对多种栅格数据的支持，包括Arc/Info ASCII Grid(asc)，GeoTiff (tiff)，Erdas Imagine Images(img)，ASCII DEM(dem) 等格式。 GDAL的安装有几种方式，可以通过下载对应Python版本的whl文件来安装，也可以使用pip命令来安装。建议使用pip或conda命令安装，比较简单方便，注意要加上版本号。 conda/pip install gdal=版本号虽然我们安装的时候直接写的是gdal，但调用时，需要从osgeo里面调用，这是由于GDAL成为了OSGeo的子项目后（1.6版本以后），GDAL被包含在osgeo这个名字空间里了。 from osgeo import gdal 遥感影像读取废话不多说，直接上代码 # -*- coding: utf-8 -*- from osgeo import gdal #读图像文件 def read_img(filename): dataset = gdal.Open(filename) #打开文件 im_width = dataset.RasterXSize #栅格矩阵的列数 im_height = dataset.RasterYSize #栅格矩阵的行数 im_bands = dataset.RasterCount #波段数 im_geotrans = dataset.GetGeoTransform() #仿射矩阵，左上角像素的大地坐标和像素分辨率 im_proj = dataset.GetProjection() #地图投影信息，字符串表示 im_data = dataset.ReadAsArray(0,0,im_width,im_height) del dataset return im_width, im_height, im_bands, im_proj, im_geotrans, im_data可以看到，使用GDAL可以很方便的对多波段数据进行读取，还可以获得其投影参数和转换参数信息，这对影像是十分重要的。 遥感影像存储 # -*- coding: utf-8 -*- from osgeo import gdal #写GeoTiff文件 def write_img(filename, im_proj, im_geotrans, im_data): #判断栅格数据的数据类型 if 'int8' in im_data.dtype.name: datatype = gdal.GDT_Byte elif 'int16' in im_data.dtype.name: datatype = gdal.GDT_UInt16 else: datatype = gdal.GDT_Float32 #判读数组维数 if len(im_data.shape) == 3: im_bands, im_height, im_width = im_data.shape else: im_bands, (im_height, im_width) = 1, im_data.shape #创建文件 driver = gdal.GetDriverByName("GTiff") dataset = driver.Create(filename, im_width, im_height, im_bands, datatype) dataset.SetGeoTransform(im_geotrans) #写入仿射变换参数 dataset.SetProjection(im_proj) #写入投影 if im_bands == 1: dataset.GetRasterBand(1).WriteArray(im_data) #写入数组数据 else: for i in range(im_bands): dataset.GetRasterBand(i+1).WriteArray(im_data[i]) del dataset这里，需要注意几点： 1、GDAL的数据类型，一般包括以下几种： #gdal数据类型 #gdal.GDT_Byte, #gdal .GDT_UInt16, gdal.GDT_Int16, gdal.GDT_UInt32, gdal.GDT_Int32, #gdal.GDT_Float32, gdal.GDT_Float642、注意数据的存储波段顺序：im_bands, im_height, im_width。 3、创建文件时 driver = gdal.GetDriverByName("GTiff")，数据类型必须要指定，因为要计算需要多大内存空间。 总结这里，我们可以影像读取和存储做成一个类文件，这样后续就可以方便的进行调用啦！ # -*- coding: utf-8 -*- from osgeo import gdal class IMAGE: #读图像文件 def read_img(self, filename): dataset = gdal.Open(filename) #打开文件 im_width = dataset.RasterXSize #栅格矩阵的列数 im_height = dataset.RasterYSize #栅格矩阵的行数 im_bands = dataset.RasterCount #波段数 im_geotrans = dataset.GetGeoTransform() #仿射矩阵，左上角像素的大地坐标和像素分辨率 im_proj = dataset.GetProjection() #地图投影信息，字符串表示 im_data = dataset.ReadAsArray(0,0,im_width,im_height) del dataset return im_width, im_height, im_bands, im_proj, im_geotrans, im_data #写GeoTiff文件 def write_img(self, filename, im_proj, im_geotrans, im_data): #判断栅格数据的数据类型 if 'int8' in im_data.dtype.name: datatype = gdal.GDT_Byte elif 'int16' in im_data.dtype.name: datatype = gdal.GDT_UInt16 else: datatype = gdal.GDT_Float32 #判读数组维数 if len(im_data.shape) == 3: im_bands, im_height, im_width = im_data.shape else: im_bands, (im_height, im_width) = 1, im_data.shape #创建文件 driver = gdal.GetDriverByName("GTiff") dataset = driver.Create(filename, im_width, im_height, im_bands, datatype) dataset.SetGeoTransform(im_geotrans) #写入仿射变换参数 dataset.SetProjection(im_proj) #写入投影 if im_bands == 1: dataset.GetRasterBand(1).WriteArray(im_data) #写入数组数据 else: for i in range(im_bands): dataset.GetRasterBand(i+1).WriteArray(im_data[i]) del dataset
|
今日新闻 |
推荐新闻 |
| CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |