How to get real Landsat image conrers

213 Views Asked by At

How I can get actual coordinates of Landsat image corners (see image to understand) ? From metadata file (..._MTL.txt) I can get coordinates of red corners, but I need to get coordinates of green corners.

I work with GeoTIFF files using GDAL.
I need to get correct latitude and longitude of green points.

Can I do it using python3? enter image description here

Thanks for help

Metadata file


    GROUP = L1_METADATA_FILE
      GROUP = METADATA_FILE_INFO
        ORIGIN = "Image courtesy of the U.S. Geological Survey"
        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
      GROUP = PRODUCT_METADATA
        PRODUCT_TYPE = "L1T"
        ELEVATION_SOURCE = "GLS2000"
        PROCESSING_SOFTWARE = "LPGS_11.3.0"
        EPHEMERIS_TYPE = "DEFINITIVE"
        SPACECRAFT_ID = "Landsat5"
        SENSOR_ID = "TM"
        SENSOR_MODE = "BUMPER"
        ACQUISITION_DATE = 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
      GROUP = 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
      GROUP = 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"
        SUN_AZIMUTH = 141.2669762
        SUN_ELEVATION = 59.9909680
        OUTPUT_FORMAT = "GEOTIFF"
      END_GROUP = PRODUCT_PARAMETERS
      GROUP = CORRECTIONS_APPLIED
        STRIPING_BAND1 = "NONE"
        STRIPING_BAND2 = "NONE"
        STRIPING_BAND3 = "NONE"
        STRIPING_BAND4 = "NONE"
        STRIPING_BAND5 = "NONE"
        STRIPING_BAND6 = "NONE"
        STRIPING_BAND7 = "NONE"
        BANDING = "N"
        COHERENT_NOISE = "N"
        MEMORY_EFFECT = "Y"
        SCAN_CORRELATED_SHIFT = "Y"
        INOPERABLE_DETECTORS = "N"
        DROPPED_LINES = "N"
      END_GROUP = CORRECTIONS_APPLIED
      GROUP = PROJECTION_PARAMETERS
        REFERENCE_DATUM = "WGS84"
        REFERENCE_ELLIPSOID = "WGS84"
        GRID_CELL_SIZE_THM = 60.000
        GRID_CELL_SIZE_REF = 30.000
        ORIENTATION = "NUP"
        RESAMPLING_OPTION = "CC"
        MAP_PROJECTION = "UTM"
      END_GROUP = PROJECTION_PARAMETERS
      GROUP = UTM_PARAMETERS
        ZONE_NUMBER = 45
      END_GROUP = UTM_PARAMETERS
    END_GROUP = L1_METADATA_FILE
    END

enter image description here

1

There are 1 best solutions below

0
On

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).