在栅格区域内选择随机像素的质心 - Python+gdal

Question

我在 WGS84 投影中有一个光栅文件，我正在尝试获取图片左侧下方光栅 GeoTIFF 区域内随机像素的坐标。首先，我计算每个像素的质心坐标（再次在 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=";")

Answer 1

您可以尝试使用图像处理技术来获取光栅的坐标。例如，这里是如何使用 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).

Answer 2

如果您只想在图像上随机选取 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=',')

之后您可以使用这些随机点的坐标来检索像素值。

在栅格区域内选择随机像素的质心 - Python+gdal

Pick random pixels' centroids within raster area - Python+gdal

python

gis

raster

gdal

python-3.x