python-3.x - 如何获得真正的 Landsat 图像conrers
问题描述
如何获得 Landsat 图像角的实际坐标(请参阅图像以了解)?从元数据文件(..._MTL.txt)我可以得到红角的坐标,但我需要得到绿角的坐标。
我使用 GDAL 处理 GeoTIFF 文件。
我需要获得正确的绿点纬度和经度。
感谢帮助
元数据文件
组 = 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 结尾
解决方案
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)
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).
推荐阅读
- c# - Asp.net core如何在viewmodel实体框架中加载相关数据
- matlab - 将温度时间序列插值到压力网格上
- android - 如何修复“错误:链接引用失败”更新 Admob 所需的 firebase 依赖项后出现此错误
- javascript - 为什么我从 vue.js 组件脚本中的简单 for 循环中得到“我未定义”错误?
- javascript - 在 while 循环中设置不允许的内容的问题
- javascript - '!!user' 是做什么用的?这条线的结果是什么?
- haskell - 了解 Haskell 中的键值构造函数
- javascript - 观察现有对象的对象属性变化
- angular - 如何转换此 DateTime 使其在 Angular Material 中看起来不错
- python - 如何用 mock.patch 覆盖 Django 视图进行测试?