0

如何获得 Landsat 图像角的实际坐标（请参阅图像以了解）？从元数据文件（..._MTL.txt）我可以得到红角的坐标，但我需要得到绿角的坐标。

我使用 GDAL 处理 GeoTIFF 文件。
我需要获得正确的绿点纬度和经度。

我可以使用python3吗？ 在此处输入图像描述

感谢帮助

元数据文件

    组 = L1_METADATA_FILE
      组 = METADATA_FILE_INFO
        ORIGIN =“图片由美国地质调查局提供”
        REQUEST_ID = "9991103150002_00325"
        PRODUCT_CREATION_TIME = 2011-03-16T20:14:24Z
        STATION_ID = "EDC"
        LANDSAT5_XBAND = "1"
        GROUND_STATION = "IKR"
        LPS_PROCESSOR_NUMBER = 0
        DATEHOUR_CONTACT_PERIOD = "1016604"
        SUBINTERVAL_NUMBER = "01"
      END_GROUP = METADATA_FILE_INFO
      组 = PRODUCT_METADATA
        PRODUCT_TYPE = "L1T"
        ELEVATION_SOURCE = "GLS2000"
        PROCESSING_SOFTWARE = "LPGS_11.3.0"
        EPHEMERIS_TYPE = "确定"
        SPACECRAFT_ID = "Landsat5"
        SENSOR_ID = "TM"
        SENSOR_MODE = "保险杠"
        收购日期 = 2010-06-15
        SCENE_CENTER_SCAN_TIME = 04:57:44.2830500Z
        WRS_PATH = 145
        STARTING_ROW = 26
        ENDING_ROW = 26
        BAND_COMBINATION = "1234567"
        PRODUCT_UL_CORNER_LAT = 49.8314223
        PRODUCT_UL_CORNER_LON = 84.0018859
        PRODUCT_UR_CORNER_LAT = 49.8694055
        PRODUCT_UR_CORNER_LON = 87.4313889
        PRODUCT_LL_CORNER_LAT = 47.8261840
        PRODUCT_LL_CORNER_LON = 84.1192898
        PRODUCT_LR_CORNER_LAT = 47.8615913
        PRODUCT_LR_CORNER_LON = 87.4144676
        PRODUCT_UL_CORNER_MAPX = 284400.000
        PRODUCT_UL_CORNER_MAPY = 5524200.000
        PRODUCT_UR_CORNER_MAPX = 531000.000
        PRODUCT_UR_CORNER_MAPY = 5524200.000
        PRODUCT_LL_CORNER_MAPX = 284400.000
        PRODUCT_LL_CORNER_MAPY = 5301000.000
        PRODUCT_LR_CORNER_MAPX = 531000.000
        PRODUCT_LR_CORNER_MAPY = 5301000.000
        PRODUCT_SAMPLES_REF = 8221
        PRODUCT_LINES_REF = 7441
        PRODUCT_SAMPLES_THM = 4111
        PRODUCT_LINES_THM = 3721
        BAND1_FILE_NAME = "L5145026_02620100615_B10.TIF"
        BAND2_FILE_NAME = "L5145026_02620100615_B20.TIF"
        BAND3_FILE_NAME = "L5145026_02620100615_B30.TIF"
        BAND4_FILE_NAME = "L5145026_02620100615_B40.TIF"
        BAND5_FILE_NAME = "L5145026_02620100615_B50.TIF"
        BAND6_FILE_NAME = "L5145026_02620100615_B60.TIF"
        BAND7_FILE_NAME = "L5145026_02620100615_B70.TIF"
        GCP_FILE_NAME = "L5145026_02620100615_GCP.txt"
        METADATA_L1_FILE_NAME = "L5145026_02620100615_MTL.txt"
        CPF_FILE_NAME = "L5CPF20100401_20100630_09"
      END_GROUP = PRODUCT_METADATA
      组 = MIN_MAX_RADIANCE
        LMAX_BAND1 = 193.000
        LMIN_BAND1 = -1.520
        LMAX_BAND2 = 365.000
        LMIN_BAND2 = -2.840
        LMAX_BAND3 = 264.000
        LMIN_BAND3 = -1.170
        LMAX_BAND4 = 221.000
        LMIN_BAND4 = -1.510
        LMAX_BAND5 = 30.200
        LMIN_BAND5 = -0.370
        LMAX_BAND6 = 15.303
        LMIN_BAND6 = 1.238
        LMAX_BAND7 = 16.500
        LMIN_BAND7 = -0.150
      END_GROUP = MIN_MAX_RADIANCE
      组 = MIN_MAX_PIXEL_VALUE
        QCALMAX_BAND1 = 255.0
        QCALMIN_BAND1 = 1.0
        QCALMAX_BAND2 = 255.0
        QCALMIN_BAND2 = 1.0
        QCALMAX_BAND3 = 255.0
        QCALMIN_BAND3 = 1.0
        QCALMAX_BAND4 = 255.0
        QCALMIN_BAND4 = 1.0
        QCALMAX_BAND5 = 255.0
        QCALMIN_BAND5 = 1.0
        QCALMAX_BAND6 = 255.0
        QCALMIN_BAND6 = 1.0
        QCALMAX_BAND7 = 255.0
        QCALMIN_BAND7 = 1.0
      END_GROUP = MIN_MAX_PIXEL_VALUE
      GROUP = PRODUCT_PARAMETERS
        CORRECTION_METHOD_GAIN_BAND1 = "CPF"
        CORRECTION_METHOD_GAIN_BAND2 = "CPF"
        CORRECTION_METHOD_GAIN_BAND3 = "CPF"
        CORRECTION_METHOD_GAIN_BAND4 = "CPF"
        CORRECTION_METHOD_GAIN_BAND5 = "CPF"
        CORRECTION_METHOD_GAIN_BAND6 = "IC"
        CORRECTION_METHOD_GAIN_BAND7 = "CPF"
        CORRECTION_METHOD_BIAS = "IC"
        太阳方位角 = 141.2669762
        SUN_ELEVATION = 59.9909680
        OUTPUT_FORMAT = "GEOTIFF"
      END_GROUP = PRODUCT_PARAMETERS
      GROUP = CORRECTIONS_APPLIED
        STRIPING_BAND1 = "无"
        STRIPING_BAND2 = "无"
        STRIPING_BAND3 = "无"
        STRIPING_BAND4 = "无"
        STRIPING_BAND5 = "无"
        STRIPING_BAND6 = "无"
        STRIPING_BAND7 = "无"
        条带=“N”
        COHERENT_NOISE = "N"
        MEMORY_EFFECT = "是"
        SCAN_CORRELATED_SHIFT = "是"
        INOPERABLE_DETECTORS = "N"
        DROPPED_LINES = "N"
      END_GROUP = CORRECTIONS_APPLIED
      组 = PROJECTION_PARAMETERS
        REFERENCE_DATUM = "WGS84"
        REFERENCE_ELLIPSOID = "WGS84"
        GRID_CELL_SIZE_THM = 60.000
        GRID_CELL_SIZE_REF = 30.000
        方向=“NUP”
        RESAMPLING_OPTION = "抄送"
        MAP_PROJECTION = "UTM"
      END_GROUP = PROJECTION_PARAMETERS
      组 = UTM_PARAMETERS
        ZONE_NUMBER = 45
      END_GROUP = UTM_PARAMETERS
    END_GROUP = L1_METADATA_FILE
    结尾

在此处输入图像描述

4

1 回答 1

0

You might first find the contour with the biggest area. Then try some algorithm to find the points you want. It seems that the satellite picture in the image is not a perfect rectangle, so you can't fit a rectangle on it using OpenCV's built-in methods.

You should try something like that:

import cv2
import numpy as np

img = cv2.imread('z_edited.jpg')
imgray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
blurred = cv2.GaussianBlur(imgray, (11, 11), 0)
ret, thresh = cv2.threshold(blurred, 27, 255, 0)
cnts, _ = cv2.findContours(thresh, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
max_area = 0
max_area_index = 0
for i, cnt in enumerate(cnts):
    area = cv2.contourArea(cnt)
    if area > max_area:
        max_area = area
        max_area_index = i

x_min = np.min(cnts[max_area_index][:, 0, 0])
x_max = np.max(cnts[max_area_index][:, 0, 0])
y_min = np.min(cnts[max_area_index][:, 0, 1])
y_max = np.max(cnts[max_area_index][:, 0, 1])

(x_left, y_left) = (x_min, cnts[max_area_index][np.max(np.where(cnts[max_area_index][:, 0, 0] == x_min)), 0, 1])
(x_right, y_right) = (x_max, cnts[max_area_index][np.max(np.where(cnts[max_area_index][:, 0, 0] == x_max)), 0, 1])
(x_down, y_down) = (cnts[max_area_index][np.max(np.where(cnts[max_area_index][:, 0, 1] == y_max)), 0, 0], y_max)
(x_top, y_top) = (cnts[max_area_index][np.max(np.where(cnts[max_area_index][:, 0, 1] == y_min)), 0, 0], y_min)


cv2.circle(img, (x_left, y_left), 10, (0, 0, 255), thickness=8)
cv2.circle(img, (x_right, y_right), 10, (0, 0, 255), thickness=8)
cv2.circle(img, (x_down, y_down), 10, (0, 0, 255), thickness=8)
cv2.circle(img, (x_top, y_top), 10, (0, 0, 255), thickness=8)
# cv2.drawContours(img, cnts, max_area_index, (0, 255, 0), 2)
cv2.namedWindow('s', cv2.WINDOW_NORMAL)
cv2.imshow('s', img)
cv2.waitKey(0)

And the result looks like: enter image description here

Using this code you can find the coordinates of the corners of the satellite picture inside the image(red points).

Also need to say I have assumed that your satellite picture background is completely black(the image you have uploaded, has a thin gray strip around the whole image).

于 2020-12-01T15:16:56.563 回答