1-概要

提示：本教程介绍如何使用 Python 中的 GDAL 库将栅格数据读取为数组并将数组另存为GeoTiff 文件

视频地址：B站对应教程

目标：
（1）读写GeoTiff影像；
（2）正确显示单波段影像；
（3）学会制作掩膜影像用于提取影像中有用部分

2.代码实现

实现主要功能如下

1. Geotiff文件读取并显示
2. 数据操作-制作掩膜；
3. 输出影像

（1）代码：Geotiff文件读取并显示

# -*- coding:utf-8 -*-
from osgeo import gdal
import numpy as np
import matplotlib.pyplot as plt
# 1.加载并显示影像 
ds = gdal.Open("dem.tif")# gdal.Dataset
gt = ds.GetGeoTransform()#地理变换：图像到地理坐标的变换矩阵
proj = ds.GetProjection()#投影信息band = ds.GetRasterBand(1)#查询波段数：ds.RasterCount
array = band.ReadAsArray()plt.figure()
plt.imshow(array)

（2）代码：数据操作-制作掩膜

# -*- coding:utf-8 -*-
from osgeo import gdal
import numpy as np
import matplotlib.pyplot as plt
# 1.加载并显示影像 
ds = gdal.Open("dem.tif")# gdal.Dataset
gt = ds.GetGeoTransform()#地理变换：图像到地理坐标的变换矩阵
proj = ds.GetProjection()#投影信息band = ds.GetRasterBand(1)#查询波段数：ds.RasterCount
array = band.ReadAsArray()# 2.数据操作：过滤高于数据平均值的值，
#高于的赋值为1，否则为0，组成了掩膜图像binmask
binmask = np.where((array  >= np.mean(array)),1,0)#制作二进制掩膜，用于提取影像
plt.figure()
plt.imshow(binmask)

（3）代码：输出影像

# -*- coding:utf-8 -*-
from osgeo import gdal
import numpy as np
import matplotlib.pyplot as plt
# 1.加载并显示影像 
ds = gdal.Open("dem.tif")# gdal.Dataset
gt = ds.GetGeoTransform()#地理变换：图像到地理坐标的变换矩阵
proj = ds.GetProjection()#投影信息band = ds.GetRasterBand(1)#查询波段数：ds.RasterCount
array = band.ReadAsArray()plt.figure()
plt.imshow(array)# 2.数据操作：过滤高于数据平均值的值，高于赋值为1，否则为0
binmask = np.where((array  >= np.mean(array)),1,0)#制作二进制掩膜，用于提取影像
plt.figure()
plt.imshow(binmask)#3.输出影像
driver = gdal.GetDriverByName("GTiff")
driver.Register()
outds = driver.Create("binmask.tif", xsize = binmask.shape[1],ysize = binmask.shape[0], bands = 1, eType = gdal.GDT_Int16)
outds.SetGeoTransform(gt)
outds.SetProjection(proj)
outband = outds.GetRasterBand(1)
outband.WriteArray(binmask)
outband.SetNoDataValue(np.nan)
outband.FlushCache()# 4.关闭数据和波段
outband = None
outds = No