3_读取遥感影像信息

最新推荐文章于 2023-11-01 14:55:19 发布

Hexk_

最新推荐文章于 2023-11-01 14:55:19 发布

阅读量1.8k

点赞数

分类专栏： GDAL 文章标签： python 开发语言

本文链接：https://blog.csdn.net/Asphy/article/details/127709048

版权

GDAL 专栏收录该内容

15 篇文章 3 订阅

订阅专栏

读取遥感影像信息

本文中的所有方法可以读取遥感影像信息，但是仅限于概括地获取元数据信息，无法具体到对像元进行处理。

打开已有的GeoTiff

GDAL不能通过with打开影像，会报错

GDAL打开影像的方式很简单，直接gdal.open(...)就可以打开影像。

不需要声明上一章所用的driver。

l1b1 = gdal.Open('I:...\\N11_RPRO_AVH_L1B_1P_19920417T060234_19920417T060633_018353\\image.l1b')
tif1 = gdal.Open('I:\\MasterLife\\资源与环境遥感\\实验作业一\\caijian2.tif')

通过dir(l1b1)可以查看打开影像的数据集下所有的方法。

查看影像描述信息

l1b1.GetDescription()方法可以查看当前影像的描述信息，一般是返回影像的调用路径。

l1b1.GetDescription()
#  D:\Anaconda\envs\GDAL_Learning\python.exe I:\PythonLearning\Learning\GDAL_Learning\demo3_打开GeoTiff.py 
#  I:\MasterLife\遥感机理\结课作业\N11_RPRO_AVH_L1B_1P_19920417T060234_19920417T060633_018353_v0100\N11_RPRO_AVH_L1B_1P_19920417T060234_19920417T060633_018353\N11_RPRO_AVH_L1B_1P_19920417T060234_19920417T060633_018353\image.l1b

查看影像元数据

l1b1.GetMetadata()可以获取影像的元数据，不同影像的元数据结果有所差异。

l1b1.GetMetadata()
# {'DATASET_NAME': 'BRN.HRPT.NH.D92108.S0607.E0603.B1835353.MN', 'DATA_TYPE': 'AVHRR HRPT', 'LOCATION': 'Ascending', 'PROCESSING_CENTER': 'Unknown processing center', 'REVOLUTION': '18353', 'SATELLITE': 'NOAA-11(H)', 'SOURCE': 'Unknown receiving station', 'START': 'year: 1992, day: 108, millisecond: 21754153', 'STOP': 'year: 1992, day: 108, millisecond: 21993986'}

查看影像波段信息和像元数目

属性：

RasterCount获取影像集的波段数目

RasterXSize获取影像集的X方向像元数目

RasterYSize获取影像集的Y方向像元数目

print(l1b1.RasterCount)  # 5
print(l1b1.RasterXSize)  # 2048
print(l1b1.RasterYSize)  # 1440

获取空间参考信息

获取影像的地理坐标仿射变换参数。

print(tif1.GetGeoTransform()) 
# (435435.0, 30.0, 0.0, 4467435.0, 0.0, -30.0)

获取影像投影信息

print(tif1.GetProjection())
# PROJCS["WGS 84 / UTM zone 50N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32650"]]

获取栅格集的波段

# 获取栅格集的波段
band1 = l1b1.GetRasterBand(1)
print(band1)
# <osgeo.gdal.Band; proxy of <Swig Object of type 'GDALRasterBandShadow *' at 0x00000253EF7174E0> >

获取波段的大小

print(band1.XSize, band1.YSize)  # 2048 1440
print(band1.DataType)  # 2

获取波段数据的属性

Maximum和Minimum分别表示band中的最大和最小值，但是直接print(band1.GetMaximum())却并没有结果。

根据解释：因为对于文件格式不会有固有的最大最小值。

所以通过：ComputeRasterMinMax()计算最大最小值。

band1.GetNoDataValue()
print(band1.GetMaximum())  # None
print(band1.ComputeRasterMinMax())  # (66.0, 695.0)
print(band1.GetMaximum())  # None

存在疑问：GetNoDataValue()和GetMaximum()这种函数有什么用处？用在什么地方？

读取ENVI文件

ENVI格式的文件包含HDR头文件和一个dat文件，所以在gdal里面，确定目标文件是ENVI格式的后，直接读取，不用加后缀即可。

envi_file = gdal.Open("I:\\MasterLife\\资源与环境遥感\\Finial\\b19")
print(envi_file.RasterCount)  # 1

读取HDF数据文件

由于modis卫星数据跟我们经常遇到的geotif数据组织方式不一样，读取的时候一定要特别注意。

geotif数据，一般是一个文件，包含了多个波段的数据

而modis，一个文件包含了多各SUBDATASETSGDAL，每个SUBDATASETS又包含多个波段数据。

另外默认编译的GDAL并不包含对MODIS数据支持，需要单独下载针对HDF4，HDF5的源码，再修改下make.opt文件，这时再编译GDAL，就支持modis数据的读写了。

这里我没有尝试