实习背景与目的¶

随着遥感技术和数字图像处理在环境监测、资源调查、城市规划等领域的广泛应用，掌握从环境配置到图像分析的完整流程尤为重要。本次实习以“图像读写与直方图处理”为核心，旨在让学生：

熟悉基于 Conda 的 Python 环境搭建与第三方库管理； 掌握遥感影像（GeoTIFF、HDF）读取、显示与保存方法； 理解并实现影像灰度/彩色直方图统计、二值化、均衡化和匹配等常用算法； 养成编写规范代码、撰写实验报告和整理成果的良好习惯。

实习内容与步骤¶

环境搭建 1.1 安装 Miniconda 并添加国内镜像源； 1.2 创建并激活名为 env_name 的虚拟环境（Python 3.7）； 1.3 安装所需库：GDAL、NumPy、OpenCV、Matplotlib。

图像读写与显示 2.1 使用 GDAL 读取 GeoTIFF（.tif）和 HDF（.hdf）数据，输出尺寸、波段数、投影与地理变换参数； 2.2 利用 Matplotlib 显示单波段或 RGB 组合图像； 2.3 构造 100×100×3 的随机 uint8 数组，设置自定义地理变换与投影，生成 out.tif 并保存。

直方图统计与二值化 3.1 以灰度模式读取 source.tif 和 reference.tif，调用 OpenCV 的 cv2.calcHist() 计算并绘制直方图； 3.2 对 source.tif 按指定阈值执行二值化（cv2.threshold）并可视化结果；

直方图均衡与匹配 4.1 调用 cv2.equalizeHist() 对灰度图像进行全局直方图均衡，并输出均衡后直方图； 4.2 实现参考影像（reference.tif）对源影像（source.tif）的直方图匹配： a) 分别计算两幅图的累积分布函数（CDF）； b) 构建像素映射表； c) 应用映射并展示匹配结果； 4.3 将（3）和（4）的所有输出在同一窗口分子图中对比展示。

（1）图像读取

任务：读取data目录下的wuhan.tif和MODIS数据集(MYDxxx.hdf)，输出影像和数据集的基本信息(wuhan.tif：影像尺寸，波段数，地理变化参数，投影信息，数据类型；MYDxxx.hdf：子数据集数量，各个子数据集的描述，某个子数据集的尺寸、波段数)。

In [1]:
from osgeo import gdal
import numpy as np
tiff_file = r"data\wuhan.tif"
hdf_file = r"data\MYD021KM.A2017209.0520.061.2018034042549.hdf"

ds = gdal.Open(r"data\wuhan.tif")
print("=== TIFF 基本信息 ===")
print("图像尺寸 (宽 x 高)：", ds.RasterXSize, "x", ds.RasterYSize)
print("波段数：", ds.RasterCount)
print("数据类型：", gdal.GetDataTypeName(ds.GetRasterBand(1).DataType))
print("地理变换参数 (GT)：", ds.GetGeoTransform())
print("投影信息 WKT：", ds.GetProjection())
print("波段1统计信息：", ds.GetRasterBand(1).GetStatistics(False, True))
print("波段2统计信息：", ds.GetRasterBand(2).GetStatistics(False, True))
print("波段3统计信息：", ds.GetRasterBand(3).GetStatistics(False, True))
print()

ds = gdal.Open(hdf_file)
subds = ds.GetSubDatasets()
print("=== HDF 子数据集信息 ===")
print("子数据集总数：", len(subds))
for idx, (name, desc) in enumerate(subds):
    print(f"  [{idx}] 名称：{name}")
    print(f"      描述：{desc}")

if subds:
    first_name = subds[0][0]
    sub = gdal.Open(first_name)
    print("\n--- 第一个子数据集详细信息 ---")
    print("子数据集路径：", first_name)
    print("尺寸 (宽 x 高)：", sub.RasterXSize, "x", sub.RasterYSize)
    print("波段数：", sub.RasterCount)
    print("数据类型：", gdal.GetDataTypeName(sub.GetRasterBand(1).DataType))

=== TIFF 基本信息 ===
图像尺寸 (宽 x 高)： 1501 x 1501
波段数： 3
数据类型： UInt16
地理变换参数 (GT)： (771345.0, 30.0, 0.0, 3375045.0, 0.0, -30.0)
投影信息 WKT： PROJCS["WGS 84 / UTM zone 49N",GEOGCS["WGS 84",DATUM["World Geodetic System 1984",SPHEROID["WGS 84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",111],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]]
波段1统计信息： [6003.0, 24519.0, 8078.3388063299, 1035.6633992489]
波段2统计信息： [6657.0, 22312.0, 8364.6720365415, 749.76280038079]
波段3统计信息： [7799.0, 21258.0, 8836.6612584726, 622.01119009827]

=== HDF 子数据集信息 ===
子数据集总数： 19
  [0] 名称：HDF4_EOS:EOS_SWATH:"data\MYD021KM.A2017209.0520.061.2018034042549.hdf":MODIS_SWATH_Type_L1B:EV_1KM_RefSB
      描述：[15x2030x1354] EV_1KM_RefSB MODIS_SWATH_Type_L1B (16-bit unsigned integer)
  [1] 名称：HDF4_EOS:EOS_SWATH:"data\MYD021KM.A2017209.0520.061.2018034042549.hdf":MODIS_SWATH_Type_L1B:EV_1KM_RefSB_Uncert_Indexes
      描述：[15x2030x1354] EV_1KM_RefSB_Uncert_Indexes MODIS_SWATH_Type_L1B (8-bit unsigned integer)
  [2] 名称：HDF4_EOS:EOS_SWATH:"data\MYD021KM.A2017209.0520.061.2018034042549.hdf":MODIS_SWATH_Type_L1B:EV_1KM_Emissive
      描述：[16x2030x1354] EV_1KM_Emissive MODIS_SWATH_Type_L1B (16-bit unsigned integer)
  [3] 名称：HDF4_EOS:EOS_SWATH:"data\MYD021KM.A2017209.0520.061.2018034042549.hdf":MODIS_SWATH_Type_L1B:EV_1KM_Emissive_Uncert_Indexes
      描述：[16x2030x1354] EV_1KM_Emissive_Uncert_Indexes MODIS_SWATH_Type_L1B (8-bit unsigned integer)
  [4] 名称：HDF4_EOS:EOS_SWATH:"data\MYD021KM.A2017209.0520.061.2018034042549.hdf":MODIS_SWATH_Type_L1B:EV_250_Aggr1km_RefSB
      描述：[2x2030x1354] EV_250_Aggr1km_RefSB MODIS_SWATH_Type_L1B (16-bit unsigned integer)
  [5] 名称：HDF4_EOS:EOS_SWATH:"data\MYD021KM.A2017209.0520.061.2018034042549.hdf":MODIS_SWATH_Type_L1B:EV_250_Aggr1km_RefSB_Uncert_Indexes
      描述：[2x2030x1354] EV_250_Aggr1km_RefSB_Uncert_Indexes MODIS_SWATH_Type_L1B (8-bit unsigned integer)
  [6] 名称：HDF4_EOS:EOS_SWATH:"data\MYD021KM.A2017209.0520.061.2018034042549.hdf":MODIS_SWATH_Type_L1B:EV_250_Aggr1km_RefSB_Samples_Used
      描述：[2x2030x1354] EV_250_Aggr1km_RefSB_Samples_Used MODIS_SWATH_Type_L1B (8-bit integer)
  [7] 名称：HDF4_EOS:EOS_SWATH:"data\MYD021KM.A2017209.0520.061.2018034042549.hdf":MODIS_SWATH_Type_L1B:EV_500_Aggr1km_RefSB
      描述：[5x2030x1354] EV_500_Aggr1km_RefSB MODIS_SWATH_Type_L1B (16-bit unsigned integer)
  [8] 名称：HDF4_EOS:EOS_SWATH:"data\MYD021KM.A2017209.0520.061.2018034042549.hdf":MODIS_SWATH_Type_L1B:EV_500_Aggr1km_RefSB_Uncert_Indexes
      描述：[5x2030x1354] EV_500_Aggr1km_RefSB_Uncert_Indexes MODIS_SWATH_Type_L1B (8-bit unsigned integer)
  [9] 名称：HDF4_EOS:EOS_SWATH:"data\MYD021KM.A2017209.0520.061.2018034042549.hdf":MODIS_SWATH_Type_L1B:EV_500_Aggr1km_RefSB_Samples_Used
      描述：[5x2030x1354] EV_500_Aggr1km_RefSB_Samples_Used MODIS_SWATH_Type_L1B (8-bit integer)
  [10] 名称：HDF4_EOS:EOS_SWATH:"data\MYD021KM.A2017209.0520.061.2018034042549.hdf":MODIS_SWATH_Type_L1B:Height
      描述：[406x271] Height MODIS_SWATH_Type_L1B (16-bit integer)
  [11] 名称：HDF4_EOS:EOS_SWATH:"data\MYD021KM.A2017209.0520.061.2018034042549.hdf":MODIS_SWATH_Type_L1B:SensorZenith
      描述：[406x271] SensorZenith MODIS_SWATH_Type_L1B (16-bit integer)
  [12] 名称：HDF4_EOS:EOS_SWATH:"data\MYD021KM.A2017209.0520.061.2018034042549.hdf":MODIS_SWATH_Type_L1B:SensorAzimuth
      描述：[406x271] SensorAzimuth MODIS_SWATH_Type_L1B (16-bit integer)
  [13] 名称：HDF4_EOS:EOS_SWATH:"data\MYD021KM.A2017209.0520.061.2018034042549.hdf":MODIS_SWATH_Type_L1B:Range
      描述：[406x271] Range MODIS_SWATH_Type_L1B (16-bit unsigned integer)
  [14] 名称：HDF4_EOS:EOS_SWATH:"data\MYD021KM.A2017209.0520.061.2018034042549.hdf":MODIS_SWATH_Type_L1B:SolarZenith
      描述：[406x271] SolarZenith MODIS_SWATH_Type_L1B (16-bit integer)
  [15] 名称：HDF4_EOS:EOS_SWATH:"data\MYD021KM.A2017209.0520.061.2018034042549.hdf":MODIS_SWATH_Type_L1B:SolarAzimuth
      描述：[406x271] SolarAzimuth MODIS_SWATH_Type_L1B (16-bit integer)
  [16] 名称：HDF4_EOS:EOS_SWATH:"data\MYD021KM.A2017209.0520.061.2018034042549.hdf":MODIS_SWATH_Type_L1B:gflags
      描述：[406x271] gflags MODIS_SWATH_Type_L1B (8-bit unsigned integer)
  [17] 名称：HDF4_EOS:EOS_SWATH:"data\MYD021KM.A2017209.0520.061.2018034042549.hdf":MODIS_SWATH_Type_L1B:EV_Band26
      描述：[2030x1354] EV_Band26 MODIS_SWATH_Type_L1B (16-bit unsigned integer)
  [18] 名称：HDF4_EOS:EOS_SWATH:"data\MYD021KM.A2017209.0520.061.2018034042549.hdf":MODIS_SWATH_Type_L1B:EV_Band26_Uncert_Indexes
      描述：[2030x1354] EV_Band26_Uncert_Indexes MODIS_SWATH_Type_L1B (8-bit unsigned integer)

--- 第一个子数据集详细信息 ---
子数据集路径： HDF4_EOS:EOS_SWATH:"data\MYD021KM.A2017209.0520.061.2018034042549.hdf":MODIS_SWATH_Type_L1B:EV_1KM_RefSB
尺寸 (宽 x 高)： 1354 x 2030
波段数： 15
数据类型： UInt16

d:\projects\遥感数字图像处理\code\.venv\Lib\site-packages\osgeo\gdal.py:330: FutureWarning: Neither gdal.UseExceptions() nor gdal.DontUseExceptions() has been explicitly called. In GDAL 4.0, exceptions will be enabled by default.
  warnings.warn(

（2）图像显示

任务：读取data目录下的wuhan.tif影像，使用matplotlib库显示该影像。

In [2]:
import rasterio
from rasterio.plot import show
import matplotlib.pyplot as plt

with rasterio.open('data/wuhan.tif') as src:
    r = src.read(1).astype(float)
    g = src.read(2).astype(float)
    b = src.read(3).astype(float)


def norm(arr):
    arr = arr - arr.min()
    return arr / arr.max()

rgb = np.dstack([norm(r), norm(g), norm(b)])

ds = rasterio.open('data\MYD021KM.A2017209.0520.061.2018034042549.hdf')
subdataset = ds.subdatasets[0]
with rasterio.open(subdataset) as subsrc:
    band1 = subsrc.read(1)

plt.figure(figsize=(10, 4))

plt.subplot(1, 2, 1)
plt.title("Wuhan RGB Composite")
plt.imshow(rgb)
plt.axis('off')
plt.subplot(1, 2, 2)
plt.title("MODIS HDF Subdataset Band 1")
plt.imshow(band1, cmap='gray')
plt.axis('off')
plt.tight_layout()
plt.show()

<>:17: SyntaxWarning: "\M" is an invalid escape sequence. Such sequences will not work in the future. Did you mean "\\M"? A raw string is also an option.
<>:17: SyntaxWarning: "\M" is an invalid escape sequence. Such sequences will not work in the future. Did you mean "\\M"? A raw string is also an option.
C:\Users\Administrator\AppData\Local\Temp\ipykernel_15944\2022487637.py:17: SyntaxWarning: "\M" is an invalid escape sequence. Such sequences will not work in the future. Did you mean "\\M"? A raw string is also an option.
  ds = rasterio.open('data\MYD021KM.A2017209.0520.061.2018034042549.hdf')
d:\projects\遥感数字图像处理\code\.venv\Lib\site-packages\rasterio\__init__.py:350: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.
  dataset = DatasetReader(path, driver=driver, sharing=sharing, **kwargs)
No description has been provided for this image

（3）图像保存

任务：使用0~255之间的uint8格式随机数构建一个尺寸大小为100*100的三波段数据，自定义地理变换参数与投影信息，保存到data/out.tif

In [3]:
from osgeo import osr
bands = 3
width = 100
height = 100
data = np.random.randint(256, size=(bands, height, width), dtype=np.uint8)

driver = gdal.GetDriverByName('GTiff')
ds = driver.Create("data/out.tif", width, height, bands, gdal.GDT_Byte)
wuhan_ds = gdal.Open(tiff_file)

ds.SetGeoTransform(wuhan_ds.GetGeoTransform())
ds.SetProjection(wuhan_ds.GetProjection())

for i in range(bands):
    band = ds.GetRasterBand(i+1)
    band.WriteArray(data[i])

ds.FlushCache()

# 显示新创建的图像
with rasterio.open('data/out.tif') as src:
    r = src.read(1).astype(float)
    g = src.read(2).astype(float)
    b = src.read(3).astype(float)
rgb = np.dstack([norm(r), norm(g), norm(b)])
plt.figure(figsize=(8,8))
plt.imshow(rgb)
plt.axis('off')
plt.title('Random RGB')
plt.show()
No description has been provided for this image

（4）直方图统计

任务：以灰度模式分别读取data目录下的source.tif和reference.tif影像，使用opencv库计算图像的直方图。

In [4]:
import cv2
import matplotlib.pyplot as plt


src_path = 'data/source.tif'
ref_path = 'data/reference.tif'
img_src = cv2.imread(src_path, cv2.IMREAD_GRAYSCALE)
img_ref = cv2.imread(ref_path, cv2.IMREAD_GRAYSCALE)


hist_src = cv2.calcHist([img_src], [0], None, [256], [0, 256])
hist_ref = cv2.calcHist([img_ref], [0], None, [256], [0, 256])

cv2.normalize(hist_src, hist_src, 0, 1, cv2.NORM_MINMAX)
cv2.normalize(hist_ref, hist_ref, 0, 1, cv2.NORM_MINMAX)


plt.figure(figsize=(10, 4))

plt.subplot(1, 2, 1)
plt.title('Source Image Histogram')
plt.plot(hist_src, color='gray')
plt.xlim([0, 255])
plt.xlabel('gray level')
plt.ylabel('Normalized Frequency')

plt.subplot(1, 2, 2)
plt.title('Reference Image Histogram')
plt.plot(hist_ref, color='gray')
plt.xlim([0, 255])
plt.xlabel('gray level')
plt.ylabel('Normalized Frequency')

plt.tight_layout()
plt.show()
No description has been provided for this image

（5）二值化

任务：以灰度模式读取data目录下的source.tif影像，自定义分割阈值，使用opencv库将图像二值化。

In [5]:
img_gray = cv2.imread("data/source.tif", cv2.IMREAD_GRAYSCALE)

ret, img_bin = cv2.threshold(img_gray, 128, 255, cv2.THRESH_BINARY)
print(f"使用的阈值: {ret}")

cv2.imwrite("data/source_binary.tif", img_bin)
print(f"二值化图像已保存到: data/source_binary.tif")

plt.figure(figsize=(10, 4))
plt.subplot(1, 2, 1)
plt.title("Binary Image")
plt.imshow(img_bin, cmap="gray")
plt.axis("off")
plt.subplot(1, 2, 2)
hist_bin, bins = np.histogram(img_bin.flatten(), bins=256, range=[0,256])
plt.plot(hist_bin, color="black")
plt.title("Binary Image Histogram")
plt.xlim([0, 255])
plt.xlabel("gray level")
plt.ylabel("pixel count")
plt.tight_layout()
plt.show()

使用的阈值: 128.0
二值化图像已保存到: data/source_binary.tif
No description has been provided for this image

（6）直方图均衡

任务：以灰度模式读取data目录下的source.tif影像，使用opencv库将影像进行直方图均衡化，输出均衡化后的直方图。

In [6]:
import cv2
import numpy as np
from matplotlib import pyplot as plt


src_path = "data/source.tif"
img = cv2.imread(src_path, cv2.IMREAD_GRAYSCALE)


img_eq = cv2.equalizeHist(img)

hist_eq, bins = np.histogram(img_eq.flatten(), bins=256, range=[0,256])


cv2.imwrite("data/source_hist_eq.tif", img_eq)
print(f"均衡化后图像已保存到: data/source_hist_eq.tif")


plt.figure(figsize=(6, 5))

plt.subplot(2, 2, 1)
plt.imshow(img, cmap="gray")
plt.title("source image")
plt.axis("off")

plt.subplot(2, 2, 2)
plt.imshow(img_eq, cmap="gray")
plt.title("after histogram equalization")
plt.axis("off")

plt.subplot(2, 1, 2)
plt.plot(hist_eq, color="black")
plt.title("Histogram after Equalization")
plt.xlabel("gray level")
plt.ylabel("pixel count")
plt.xlim([0, 255])

plt.tight_layout()
plt.show()

均衡化后图像已保存到: data/source_hist_eq.tif
No description has been provided for this image

（7）直方图匹配

任务：以灰度模式分别读取data目录下的source.tif和reference.tif影像，以reference.tif为参考影像，对source.tif影像进行直方图匹配。

In [7]:
import numpy as np
import rasterio

def match_histogram(source, reference, n_bins=256):

    src_values, src_counts = np.unique(source.ravel(), return_counts=True)
    ref_values, ref_counts = np.unique(reference.ravel(), return_counts=True)
    src_cdf = np.cumsum(src_counts).astype(np.float64)
    src_cdf /= src_cdf[-1]
    ref_cdf = np.cumsum(ref_counts).astype(np.float64)
    ref_cdf /= ref_cdf[-1]

    interp_ref_values = np.interp(src_cdf, ref_cdf, ref_values)

    matched = np.interp(source.ravel(), src_values, interp_ref_values)
    return matched.reshape(source.shape)


with rasterio.open("data/source.tif") as src_ds:
    src = src_ds.read(1).astype(np.float32)
    profile = src_ds.profile
with rasterio.open("data/reference.tif") as ref_ds:
    ref = ref_ds.read(1).astype(np.float32)


matched = match_histogram(src, ref)


dtype = profile["dtype"]
if np.issubdtype(np.dtype(dtype), np.integer):
    info = np.iinfo(dtype)
    matched = np.clip(matched, info.min, info.max)
matched = matched.astype(dtype)


profile.update(dtype=dtype, count=1)
with rasterio.open("data/matched.tif", "w", **profile) as dst:
    dst.write(matched, 1)

print("直方图匹配完成：data/matched.tif")


with rasterio.open("data/matched.tif") as src:
    matched_img = src.read(1).astype(float)
plt.figure(figsize=(10, 4))
plt.subplot(1, 2, 1)
plt.title("Matched Image")
plt.imshow(matched_img, cmap="gray")
plt.axis("off")
plt.subplot(1, 2, 2)
hist_matched, bins = np.histogram(matched_img.flatten(), bins=256, range=[0,256])
plt.plot(hist_matched, color="black")
plt.title("Matched Image Histogram")
plt.xlim([0, 255])
plt.xlabel("gray level")
plt.ylabel("pixel count")
plt.tight_layout()
plt.show()

直方图匹配完成：data/matched.tif

d:\projects\遥感数字图像处理\code\.venv\Lib\site-packages\rasterio\__init__.py:360: NotGeoreferencedWarning: The given matrix is equal to Affine.identity or its flipped counterpart. GDAL may ignore this matrix and save no geotransform without raising an error. This behavior is somewhat driver-specific.
  dataset = writer(
No description has been provided for this image

（7）图片整合

任务：将data目录下原始影像（source.tif）、参考影像（reference.tif）、二值化影像（source_binary.tif）、直方图均衡化影像（source_hist_eq.tif）和直方图匹配影像（matched.tif），以及他们的灰度直方图整合到一张图中显示。

In [8]:
import numpy as np
import matplotlib.pyplot as plt
import rasterio


def read_tif(path):
    with rasterio.open(path) as src:
        arr = src.read(1)
    return arr

paths = [
    'data/source.tif',
    'data/reference.tif',
    'data/source_binary.tif',
    'data/source_hist_eq.tif',
    'data/matched.tif'
]
titles = [
    'source.tif',
    'reference.tif',
    'source_binary.tif',
    'source_hist_eq.tif',
    'matched.tif'
]

imgs = [read_tif(p) for p in paths]

fig, axes = plt.subplots(2, 5, figsize=(20, 8))
fig.subplots_adjust(hspace=0.3, wspace=0.2)

for i, (img, title) in enumerate(zip(imgs, titles)):
    ax_img = axes[0, i]
    ax_img.imshow(img, cmap='gray')
    ax_img.set_title(title, fontsize=12)
    ax_img.axis('off')
    
    ax_h = axes[1, i]
    data = img.flatten()
    ax_h.hist(data, bins=256, range=(np.min(data), np.max(data)),
              color='black', alpha=0.7)
    ax_h.set_xlim(np.min(data), np.max(data))
    ax_h.set_xlabel('grey level')
    ax_h.set_ylabel('pixel count')
    ax_h.grid(True, linestyle='--', linewidth=0.5, alpha=0.5)

plt.tight_layout()
plt.savefig('data/combined_figure.png', dpi=300)
plt.show()
No description has been provided for this image

实习收获与心得¶

全流程实践加深了对遥感影像数据格式和处理方法的理解； 掌握了 Conda 环境管理、GDAL 与 OpenCV 的基本用法； 通过直方图分析深入认识了图像亮度分布及其均衡／匹配的原理与实现； 强化了代码规范编写和实验报告撰写能力，为后续更复杂的数字图像处理研究打下坚实基础。