【ChatGPT+Python】Landsat卫星图像黑边去云及旋转校正

引言

下图是一张Landsat图像的示例(右图)。我们可以明显地看到四周的黑边和倾斜的角度。这是由于卫星传感器成像导致的。一般情况下,我们是不需要去除黑边和选择的,因为这样做之后投影信息和位置信息就不正确了。但对于做深度学习图像处理任务的同学来说,可能不需要考虑此类约束。因此,出于节省存储和内存的考虑,借助ChatGPT直接帮助我们生成了黑边去云及旋转校正的Python代码,最终效果如左图。

在这里插入图片描述

代码

Prompt: 现有一RGB彩色图像,图像内仅有一个旋转的矩形物体,其周围的像素值全是0,要求裁剪出其中的矩形物体,并将其旋转校正,注意矩形物体的边界也可能有值为0的像素,输出结果保存为彩色图像。使用Python实现。

  • 安装依赖
pip install opencv-python numpy
  • 算法流程
  1. 读取图像
  2. 转换为灰度图像,然后通过阈值化操作将所有非零像素转换为白色(255),所有零像素保持为黑色(0)。
  3. 使用OpenCV的findContours方法找到最大的轮廓,这应该是你的矩形物体。
  4. 使用minAreaRect获取轮廓的最小外接矩形(包括旋转角度)。
  5. 通过旋转原始图像将矩形校正。
  6. 裁剪并保存结果图像。
  • 黑边去除及旋转校正的代码如下:
import cv2
import numpy as np

# 1. 读取图像
img = cv2.imread('input.jpg')

# 2. 转换为灰度图像
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

# 3. 找到最大的轮廓
_, thresh = cv2.threshold(gray, 1, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
contours, _ = cv2.findContours(thresh, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
contour = max(contours, key=cv2.contourArea)

# 4. 获取最小外接矩形
rect = cv2.minAreaRect(contour)
box = cv2.boxPoints(rect)
box = np.int0(box)

# 5. 旋转原始图像以校正矩形
width = int(rect[1][0])
height = int(rect[1][1])
src_pts = box.astype("float32")
dst_pts = np.array([[0, height-1],
                    [0, 0],
                    [width-1, 0],
                    [width-1, height-1]], dtype="float32")
M = cv2.getPerspectiveTransform(src_pts, dst_pts)
warped = cv2.warpPerspective(img, M, (width, height))

# 6. 显示中间过程(用于可解释性研究)
cv2.drawContours(img, [box], 0, (0, 0, 255), 2)
cv2.imshow('edge', img)
cv2.waitKey(0)
cv2.destroyAllWindows()

# 7. 裁剪并保存结果图像
output = warped[0:height, 0:width]
cv2.imwrite('output.jpg', output)

以上代码首先读取图像,并通过二值化将图像分割成背景和物体。然后找到物体的外接矩形,使用仿射变换将物体旋转到水平位置,最后对旋转后的图像进行裁剪。

但请注意,对于非常复杂或噪声很大的图像,这种方法可能不适用。你可能需要对图像进行更多的预处理,比如滤波、形态学操作等等。

实验

表1:边缘(轮廓)检测方法采用不同函数的对比结果

在这里插入图片描述

拓展

上述代码不仅对于Landsat卫星图像有效,对于普通的图像(周围是黑色填充,内接一个矩形物体)也同样有效。下面给出一个演示结果。

在这里插入图片描述

此外,上述代码中寻找最大轮廓除使用 cv2.findContours() 函数外,也可以使用如下方式替代,只不过代码的执行时间会变长,但是鲁棒性会更好,即抗噪声等干扰的能力强。

import cv2
import numpy as np

# 1. 读取图像
img = cv2.imread('input.jpg')

# 2. 转换为灰度图像
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

# 3. 找到最大的轮廓
# 使用Canny边缘检测算法检测图像中的边缘
edges = cv2.Canny(gray, 50, 150, apertureSize=3)

# 使用霍夫变换检测图像中的直线
lines = cv2.HoughLines(edges, 1, np.pi/180, 200)

# 获取检测到的直线中最长的一条
longest_line = None
max_length = 0

for line in lines:
    for rho, theta in line:
        a = np.cos(theta)
        b = np.sin(theta)
        x0 = a*rho
        y0 = b*rho
        x1 = int(x0 + 1000*(-b))
        y1 = int(y0 + 1000*(a))
        x2 = int(x0 - 1000*(-b))
        y2 = int(y0 - 1000*(a))
        length = np.sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1))
        if length > max_length:
            max_length = length
            longest_line = line

# 获取直线的端点坐标
for rho, theta in longest_line:
    a = np.cos(theta)
    b = np.sin(theta)
    x0 = a*rho
    y0 = b*rho
    x1 = int(x0 + 1000*(-b))
    y1 = int(y0 + 1000*(a))
    x2 = int(x0 - 1000*(-b))
    y2 = int(y0 - 1000*(a))

# 计算直线的角度
angle = np.arctan2(y2-y1, x2-x1) * 180 / np.pi

# 4. 获取最小外接矩形
h, w = img.shape[:2]
rect = cv2.minAreaRect(np.array([(x, y) for x in range(w) for y in range(h) if edges[y, x] > 0]))
box = cv2.boxPoints(rect)
box = np.int0(box)

# 5. 旋转原始图像以校正矩形
width = int(rect[1][0])
height = int(rect[1][1])
src_pts = box.astype("float32")
dst_pts = np.array([[0, height-1],
                    [0, 0],
                    [width-1, 0],
                    [width-1, height-1]], dtype="float32")
M = cv2.getPerspectiveTransform(src_pts, dst_pts)
warped = cv2.warpPerspective(img, M, (width, height))

# 6. 显示中间过程(用于可解释性研究)
cv2.drawContours(img, [box], 0, (0, 0, 255), 2)
cv2.imshow('edge', img)
cv2.waitKey(0)
cv2.destroyAllWindows()

# 7. 裁剪并保存结果图像
output = warped[0:height, 0:width]
cv2.imwrite('output.jpg', output)

注意

Landsat的遥感影像四个角有黑色区域,这是正常的。表明那个黑色区的地方没有数据,在实际应用中也不需要去掉四个角的黑色区域。一般我们使用shp矢量量面来对遥感影像进行剪裁,提取出我们所需要研究或者显示的区域就可以啦。但是如果你真的想去掉黑色区域的话,你可以使用重分类,把黑色的区域变成白色,这样和背景就一致了,在发布服务的时候设为白色透明就可以了。还一个问题,Landsat8遥感影像是具有准确的投影坐标系的。倾斜是因为在卫星扫描的时候就是这个样子,拿到的影像是由准确的地理坐标的,不能够把它给旋转正了。旋转正了的话,它的投影信息和位置信息就不正确了,如果感觉不好看,可以直接剪裁出来正方形,或者把相邻的影像进行拼接,然后再剪裁。

参考

https://www.zhihu.com/question/497775240/answer/2216988184

  • 2
    点赞
  • 7
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 2
    评论
Python实现Landsat大气校正可以利用遥感图像处理库,如GDAL和Rasterio等。以下是一个基本的Python代码示例,用于实现Landsat 8的大气校正: 1. 导入所需的库和模块: ```python import numpy as np import rasterio from rasterio import plot from rasterio.warp import calculate_default_transform, reproject, Resampling ``` 2. 定义校正函数: ```python def atmospheric_correction(input_path, output_path): with rasterio.open(input_path) as src: # 读取波段数据 red = src.read(4) green = src.read(3) blue = src.read(2) nir = src.read(5) # 获取元数据 transform, width, height = calculate_default_transform(src.crs, 'EPSG:4326', src.width, src.height, *src.bounds) # 创建新的校正图像 dst_profile = src.profile dst_profile.update({'crs': 'EPSG:4326', 'transform': transform, 'width': width, 'height': height}) # 大气校正公式 A = np.array([[0.3561, 0.3972, 0.3904, 0.6966], [-0.3344, -0.3544, -0.4556, 0.6966], [0.2626, 0.2141, 0.0926, 0.6830]]) # 校正计算 red_corr = (red*A[0, 0] + green*A[0, 1] + blue*A[0, 2] + nir*A[0, 3]) green_corr = (red*A[1, 0] + green*A[1, 1] + blue*A[1, 2] + nir*A[1, 3]) blue_corr = (red*A[2, 0] + green*A[2, 1] + blue*A[2, 2] + nir*A[2, 3]) # 保存校正图像 with rasterio.open(output_path, 'w', **dst_profile) as dst: dst.write(red_corr, 1) dst.write(green_corr, 2) dst.write(blue_corr, 3) ``` 3. 调用校正函数进行图像校正: ```python input_path = 'input_image.tif' output_path = 'output_image.tif' atmospheric_correction(input_path, output_path) ``` 在上述代码中,我们首先使用[`rasterio.open()`](https://rasterio.readthedocs.io/en/latest/api/rasterio.open.html)打开输入图像,然后使用[`read()`](https://rasterio.readthedocs.io/en/latest/api/rasterio.io.html#rasterio.io.DatasetReader.read)读取红、绿、蓝和近红外波段的数据。接下来,使用[`calculate_default_transform()`](https://rasterio.readthedocs.io/en/latest/api/rasterio.warp.html#rasterio.warp.calculate_default_transform)函数计算出输出图像的投影和变换信息。 然后,我们创建一个新的输出图像,并将大气校正公式应用于每个波段。最后,使用[`rasterio.open()`](https://rasterio.readthedocs.io/en/latest/api/rasterio.open.html)并结合[`write()`](https://rasterio.readthedocs.io/en/latest/api/rasterio.io.html#rasterio.io.DatasetWriter.write)方法将校正后的波段写入输出图像。 以上就是利用Python实现Landsat大气校正的基本代码示例。需要注意的是,校正公式中的参数值可以根据具体情况进行调整。此外,还可以通过其他方法来进一步改进校正效果,如大气遥感模型、气溶胶光谱等。
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

Xavier Jiezou

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值