IT数码 购物 网址 头条 软件 日历 阅读 图书馆
TxT小说阅读器
↓语音阅读,小说下载,古典文学↓
图片批量下载器
↓批量下载图片,美女图库↓
图片自动播放器
↓图片自动播放器↓
一键清除垃圾
↓轻轻一点,清除系统垃圾↓
开发: C++知识库 Java知识库 JavaScript Python PHP知识库 人工智能 区块链 大数据 移动开发 嵌入式 开发工具 数据结构与算法 开发测试 游戏开发 网络协议 系统运维
教程: HTML教程 CSS教程 JavaScript教程 Go语言教程 JQuery教程 VUE教程 VUE3教程 Bootstrap教程 SQL数据库教程 C语言教程 C++教程 Java教程 Python教程 Python3教程 C#教程
数码: 电脑 笔记本 显卡 显示器 固态硬盘 硬盘 耳机 手机 iphone vivo oppo 小米 华为 单反 装机 图拉丁
 
   -> Python知识库 -> 使用Python和GDAL处理遥感影像数据超详细教程 -> 正文阅读

[Python知识库]使用Python和GDAL处理遥感影像数据超详细教程

提示:文章末尾有强化学习代码资源 : )

前言

? ? 在本教程中,我们将学习使用 Python 和地理空间数据抽象库 GDAL 自动处理栅格数据的基本技术。

? ? 栅格文件通常用于存储地形模型和遥感数据及其衍生产品,例如植被指数和其他环境数据集。 栅格文件往往很大(例如,想象一个以 30m x 30m 分辨率覆盖全球的栅格数据集)并且通常以较小的部分(切片)交付和处理。 因此,对此类大型数据集的高效处理和分析需要自动化。 在教程中,您将学习如何读写常见的栅格格式,以及如何使用 Python 中的 GDAL/OGR API 和 GDAL 命令行实用程序对一批文件进行基本的栅格数据处理。

? ? 在本文总结部分将一些实用的python剪裁批处理,python影像镶嵌以及一些Arcpy相关代码进行链接下载,供有需求的同学进行学习与使用。


?

一、使用 GDAL 在 Python 中打开栅格文件

? ? 使用 GDAL,可以在 Python 中读取和写入多种不同的栅格格式。 导入 GDAL 模块时,Python 会自动注册所有已知的 GDAL 驱动程序以读取支持的格式。 最常见的文件格式包括 TIFF 和 GeoTIFF、ASCII Grid 和 Erdas Imagine .img 文件。

? ? Landsat 8 波段作为单独的 GeoTIFF 文件存储在原始包中。 每个波段包含来自不同电磁波谱范围的表面反射率信息。?

? ??

from osgeo import gdal

filepath = r"LandsatData/LC81910182016153LGN00_sr_band4.tif"

# Open the file:
raster = gdal.Open(filepath)

# Check type of the variable 'raster'
type(raster)

(1)? 读取栅格文件属性?

? ? 卫星图像现在作为 GDAL 数据集对象存储在变量栅格中。 让我们仔细看看文件的属性:

# Projection
raster.GetProjection()

# Dimensions
raster.RasterXSize
raster.RasterYSize

# Number of bands
raster.RasterCount

# Metadata for the raster dataset
raster.GetMetadata()

?(2)? 获得栅格影像波段

? ? ?在我们的例子中,Landsat 8 的所有波段都存储为单独的文件。 rasterCount 为 1,因为我们只打开了一个包含 Landsat 8 波段 4 的 GeoTiff。但是,卫星图像的不同波段通常堆叠在一个栅格数据集中,在这种情况下 rasterCount 会大于 1。

# Read the raster band as separate variable
band = raster.GetRasterBand(1)

# Check type of the variable 'band'
type(band)

# Data type of the values
gdal.GetDataTypeName(band.DataType)

? ? 现在我们有一个存储在变量 band 中的 GDAL Raster Band 对象。

? ? 波段的数据类型可以在像素数据类型的 GDAL 文档的帮助下进行解释。 无符号整数始终等于或大于零,有符号整数也可以存储负值。 例如,一个无符号的 16 位整数可以存储 2^16 (=65,536) 个值,范围从 0 到 65,535。

?(3) 波段统计

? ? 接下来,让我们看一下存储在 band 中的值。 在能够打印出任何信息之前,您可能需要计算栅格的统计数据。

# Compute statistics if needed
if band.GetMinimum() is None or band.GetMaximum()is None:
    band.ComputeStatistics(0)
    print("Statistics computed.")
    
# Fetch metadata for the band
band.GetMetadata()
    
# Print only selected metadata:
print ("[ NO DATA VALUE ] = ", band.GetNoDataValue()) # none
print ("[ MIN ] = ", band.GetMinimum())
print ("[ MAX ] = ", band.GetMaximum())

?

二、将栅格文件作为数值数组打开

? ? 在访问地理空间栅格数据方面,GDAL 是一个功能强大的库,但它没有提供很多计算功能。 对于更高级的计算,我们将以数值数组的形式读入栅格数据,以便使用 NumPy 库中的功能。

? ? 如果您想将现有的 Gdal 数据集或 Band 转换为 numpy 数组,您可以使用 ReadAsArray 方法将其转换:

In [ ]:
# read raster data as numeric array from Gdal Dataset
rasterArray = raster.ReadAsArray()

? ? 那有什么区别呢??

In [ ]:
#Check the datatype of variables
type(rasterArray)
type(raster)

? ?如您所见,我们现在已将相同的栅格数据存储到两种不同类型的变量中:

? -raster 是一个 Gdal 数据集

? -rasterArray 是一个 numpy 数组

? ? GDAL 库还附带一个模块 gdal_array,用作 NumPy 和 GDAL 之间的接口。 Gdal_array 直接从文件读取和写入栅格文件到 NumPy 数组:

from osgeo import gdal_array
#from osgeo import gdalnumeric

# Read raster data as numeric array from file
rasterArray = gdal_array.LoadFile(filepath)

(1)注意NODATA数据

# What is the minimum value of the array?
rasterArray.min() #-9999

? ? 如您所见,numpy 数组仍然包含原始无数据值。 如果不注意这些,任何计算都将是错误的。

import numpy as np

# Get nodata value from the GDAL band object
nodata = band.GetNoDataValue()

#Create a masked array for making calculations without nodata values
rasterArray = np.ma.masked_equal(rasterArray, nodata)
type(rasterArray)

# Check again array statistics
rasterArray.min()

? ? 现在您有一个二维数组,可以进行进一步的计算。 但是我们将使用一个非常简单的命令行工具 gdal_calc.py。?

?(2)关闭栅格数据?

? ? 如果您想释放资源并从内存中删除不必要的变量,则在代码中间关闭现有的 GDAL 对象可能很有用。

raster = None
band = None

三、GDAL 命令行实用程序

? ? 我们现在已经测试了 Python GDAL/OGR API 中用于读取和检查栅格文件的一些基本函数。 然而,GDAL 还包括其他强大的数据转换和处理功能,这些功能没有直接在库中实现。 我们将仔细研究几个这样的函数:

  • gdalwarp 用于裁剪、镶嵌、重投影和其他过程
  • gdal_merge.py 用于拼接/堆叠图像
  • gdal_calc.py 用于栅格计算

? ? 这些工具需要从终端/命令提示符或作为 Python 中的子进程运行。 我们现在将在终端窗口中快速测试这些工具。

(1) 使用 gdalwarp 裁剪图像

? ? 在其他技巧中,gdalwarp 是一个非常方便的快速裁剪图像的工具。 我们现在将练习如何基于边界框裁剪卫星图像波段。使用选项 -te 指定输出文件的所需范围:

gdalwarp -te xmin ymin xmax ymax inputfile.tif outputfile.tif

? ? 接下来,让我们同时对所有其余波段重复剪裁。 为此,我们将使用 Python 为场景中的每个光谱带生成命令。?

import glob
import os

# List filepaths for all bands in the scence
FileList = glob.glob(os.path.join(r'/home/geo/LandsatData','*band*.tif'))

# Define clipping extent
xmin, ymin, xmax,ymax = (0, 0, 0, 0) # INSERT HERE THE CORRECT COORDINATES

# Generate gdalwarp command for each band
command = ""

for fp in FileList:
    inputfile = fp
    outputfile = inputfile[:-4] + "_clip.tif"
    
    command += "gdalwarp -te %s %s %s %s %s %s \n" % (xmin, ymin, xmax, ymax, inputfile, outputfile)
    
# Write the commands to an .sh file
cmd_file = "ClipTurkufromLandsat.sh"
f = open(os.path.join(cmd_file), 'w')

f.write(command)
f.close()

注意:如果您在 Windows 环境中工作,请将 .sh 扩展名更改为 .bat。

运行上述脚本后,您的工作目录中应该有一个文件“ClipTurkufromLandsat.sh”。 打开文件(使用文本编辑器)并检查命令是否已正确写入文件。

接下来,在终端窗口中运行该文件:

bash ClipTurkufromLandsat.sh

(2) 使用 gdal_merge.py 生成layer stack

? ? 裁剪图像后,可以堆叠波段 3(绿色)、4(红色)和 5(近红外)以假彩色合成。 使用 gdal_merge.py 合并图层并使用 -separate 选项指示您希望将输入保存为输出文件中的单独波段。

让我们尝试在 python 中将命令作为子进程运行:

import os

# Define input and output files
inputfiles =  "band3_clip.tif band4_clip.tif band5_clip.tif"
outputfile =  "Landsat8_GreenRedNir.tif"

# Generate the command
command = "gdal_merge.py -separate %s -o %s" % (inputfiles, outputfile)

# Run the command. os.system() returns value zero if the command was executed succesfully
os.system(command)

(3)?Gdal_calc.py?

? ?Gdal_calc.py 是一个命令行栅格计算器,可用于栅格数据的简单重复计算。 打开终端窗口并输入:

Gdal_calc.py

? ? 然后按回车。 您应该看到有关该工具的用法和选项的说明。gdal_calc.py 的基本语法如下:

gdal_calc.py -A input1.tif - B input2.tif [other_options] --outfile=outputfile.tif

?

? ? 从其他选项中,至少要注意用于指定计算语法的参数 --calc 和用于控制输出文件大小的 --creation-option(或 --co)参数:

  • ? ?在两个输入文件的情况下 --calc="A+B" 会将文件 A 和 B 添加到一起。
  • ? ?默认情况下,输出文件往往很大,这将很快导致磁盘大小和内存出现问题。 使用 gdal_calc.py 您可以添加参数 --co="COMPRESS=LZW" 以减小输出文件的大小。

总结

? ? 本博客介绍了Python GDAL对遥感数据的处理。下面附三个代码链接供有需要的同学进行强化学习。

(1)基于Geopandas的矢量文件裁剪矢量文件方法:矢量文件剪裁矢量文件(Python)-Python文档类资源-CSDN下载

(2)

基于Python-GDAL的遥感影像镶嵌脚本:

基于Python-GDAL的遥感影像镶嵌脚本_pythongdal库镶嵌-Python文档类资源-CSDN下载

(3)基于Arcpy以及GDAL的批量分类栅格图转shp矢量文件程序:

遥感影像分类结果栅格转矢量(Python)_python栅格转矢量,分类结果转矢量-Python文档类资源-CSDN下载?

  Python知识库 最新文章
Python中String模块
【Python】 14-CVS文件操作
python的panda库读写文件
使用Nordic的nrf52840实现蓝牙DFU过程
【Python学习记录】numpy数组用法整理
Python学习笔记
python字符串和列表
python如何从txt文件中解析出有效的数据
Python编程从入门到实践自学/3.1-3.2
python变量
上一篇文章      下一篇文章      查看所有文章
加:2022-12-25 11:04:55  更:2022-12-25 11:05:29 
 
开发: C++知识库 Java知识库 JavaScript Python PHP知识库 人工智能 区块链 大数据 移动开发 嵌入式 开发工具 数据结构与算法 开发测试 游戏开发 网络协议 系统运维
教程: HTML教程 CSS教程 JavaScript教程 Go语言教程 JQuery教程 VUE教程 VUE3教程 Bootstrap教程 SQL数据库教程 C语言教程 C++教程 Java教程 Python教程 Python3教程 C#教程
数码: 电脑 笔记本 显卡 显示器 固态硬盘 硬盘 耳机 手机 iphone vivo oppo 小米 华为 单反 装机 图拉丁

360图书馆 购物 三丰科技 阅读网 日历 万年历 2024年12日历 -2024/12/26 3:08:00-

图片自动播放器
↓图片自动播放器↓
TxT小说阅读器
↓语音阅读,小说下载,古典文学↓
一键清除垃圾
↓轻轻一点,清除系统垃圾↓
图片批量下载器
↓批量下载图片,美女图库↓
  网站联系: qq:121756557 email:121756557@qq.com  IT数码
数据统计