Блог

У меня есть растровый файл в проекции WGS84, и я пытаюсь получить координаты случайных пикселей в области растрового геотифа слева внизу на рисунке. Сначала я вычисляю координаты центроида каждого пикселя (снова в WGS84), затем выбираю 100 случайных из них и экспортирую их в csv.

Проблема: я ожидаю, что точки будут находиться в области растра (слева внизу на рисунке), но они находятся далеко от нее. Это ошибка проекции или неправильный расчет координат? Что не так в моем коде?

Вот код

 # Get coordinates for each pixel centroid
geotiff = gdal.Open(path)
gt = geotiff.GetGeoTransform()
column_numbers, row_numbers, band_numbers = geotiff.RasterXSize, geotiff.RasterYSize, geotiff.RasterCount
minx = gt[0]
miny = gt[3]   column_numbers*gt[4]   row_numbers*gt[5]
maxx = gt[0]   column_numbers*gt[1]   row_numbers*gt[2]
maxy = gt[3]
pixelWidth = gt[1]
pixelHeight = -gt[5]
lonPxSz = (maxy - miny) / row_numbers
latPxSz = (maxx - minx) / column_numbers
total = np.array(geotiff.ReadAsArray())
res = []
for i in range(row_numbers):
    for j in range(column_numbers):
        res.append([[i,j]]   [data[i][j] for data in total])
coords = pd.DataFrame(res, columns=['Pair', 'Col1', 'Col2', 'Col3', 'Col4', 'Col5', 'Col6'])
coords[['Lat', 'Lon']] = pd.DataFrame(coords['Pair'].tolist(), index=coords.index)
coords["Lat"] = (coords["Lat"]   0.5) * 10 * latPxSz   miny
coords["Lon"] = (coords["Lon"]   0.5) * 10 * lonPxSz   minx
coords = coords.sample(n = 100)
coords[['Lat', 'Lon']].to_csv("coords.csv", sep=";")
 

Если вы хотите выбрать только 100 случайных точек на изображении:

 from osgeo import gdal
import numpy as np
import pandas as pd
import random


path = "image.tif"
geotiff = gdal.Open(path)
gt = geotiff.GetGeoTransform()
column_numbers, row_numbers, band_numbers = geotiff.RasterXSize, geotiff.RasterYSize, geotiff.RasterCount

minx = gt[0]
miny = gt[3]   column_numbers * gt[4]   row_numbers * gt[5]
maxx = gt[0]   column_numbers * gt[1]   row_numbers * gt[2]
maxy = gt[3]

pixelWidth = gt[1]
pixelHeight = -gt[5]
halfPixelWidth = pixelWidth / 2
halfPixelHeight = pixelHeight / 2

rand_point_x = random.sample([i for i in range(column_numbers)], 100)
rand_point_y = random.sample([i for i in range(row_numbers)], 100)
rand_points = np.vstack((rand_point_y, rand_point_x)).T

coords = pd.DataFrame(rand_points, columns=['Lat', 'Lon'])
coords["Lat"] = miny   (coords["Lat"] * pixelHeight)   halfPixelHeight
coords["Lon"] = minx   (coords["Lon"] * pixelWidth)   halfPixelWidth

coords.to_csv("coords.csv", sep=',')
 

Вы можете использовать координаты этих случайных точек для последующего извлечения значений пикселей.

Вы можете попробовать использовать методы обработки изображений, чтобы получить координаты растра. Например, вот как это можно сделать с помощью библиотеки cv2 (OpenCV) (назначение каждой функции, прокомментированной в коде):

 import cv2
import numpy as np

def process(img): # Function to process image for optimal contour detection
    img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    img_blur = cv2.GaussianBlur(img_gray, (5, 5), 1)
    img_canny = cv2.Canny(img_blur, 350, 150)
    kernel = np.ones((3, 3))
    img_dilate = cv2.dilate(img_canny, kernel, iterations=1)
    return cv2.erode(img_dilate, kernel, iterations=1)

def get_raster(img): # Function that uses process function to detect contour of raster
    contours, _ = cv2.findContours(process(img), cv2.RETR_TREE, cv2.CHAIN_APPROX_NONE)
    cnt = max(contours, key=cv2.contourArea)
    peri = cv2.arcLength(cnt, True)
    approx = cv2.approxPolyDP(cnt, 0.05 * peri, True)
    return cv2.boundingRect(approx)

def get_random(img, num=100): # Function that uses get_raster to get random points within raster
    x, y, w, h = get_raster(img)
    return np.vstack((np.random.randint(x, x   w, num),
                      np.random.randint(y, y   h, num))).T

img = cv2.imread("map.png") # Read in image
pts = get_random(img) # Get random points witin raster
cv2.drawContours(img, pts[:, None], -1, (0, 255, 0), 2) # Draw points onto image
cv2.imshow("Image", img)
cv2.waitKey(0)
 

Выход:

Как вы можете видеть, случайно расположенные зеленые точки были нарисованы на изображении в растровой проекции. Если вам нужны только координаты растра, вы можете просто сделать x, y, w, h = get_raster(img) .