Обработка изображений с использованием Python, GDAL и Scikit-Image

Я борюсь с обработкой и, надеюсь, я смогу решить здесь.

Я работаю с дистанционным зондированием, примененным к лесному хозяйству, особенно работаю с данными LiDAR. Идея состоит в том, чтобы использовать Scikit-изображение для определения вершины дерева. Так как я новичок в Python, я считаю большим личным триумфом сделать следующее:

1. Импортировать CHM (с matplotlib);
2. Запустите гауссов фильтр (с пакетом scikit-image);
3. Запустите фильтр максимума (с пакетом scikit-image);
4. Запустите peak_local_max (с пакетом scikit-image);
5. Показать CHM с локальными максимумами (с matplotlib);

Теперь моя проблема. Когда я импортирую с помощью matplot, изображение теряет свои географические координаты. Таким образом, координаты, которые у меня есть, являются просто базовыми координатами изображения (то есть 250 312). Что мне нужно, это получить значение пикселя под локальной точкой максимума на изображении (красные точки на изображении). Здесь, на форуме, я видел одного парня, спрашивающего то же самое ( Получить значение пикселя растра GDAL под точкой OGR без NumPy? ), Но у него уже были точки в шейп-файле. В моем случае точки были вычислены с помощью scikit-image (это массив с координатами каждой вершины дерева). Так что у меня нет шейп файла.

В заключение, в конце я хочу получить текстовый файл с координатами каждого локального максимума в географических координатах, например:

525412 62980123 1150 ...

Во-первых, добро пожаловать на сайт!

В массивах Numpy отсутствует концепция систем координат, встроенных в массив. Для 2D-растра они индексируются по столбцам и строкам.

Заметьте, я предполагаю, что вы читаете растровый формат, который поддерживается GDAL .

В Python лучший способ импортировать пространственные растровые данные - это rasterioпакет. Необработанные данные, импортированные растерио, по-прежнему являются массивом без доступа к системам координат, но растерио также предоставляет вам доступ к аффинному методу в исходном массиве, который можно использовать для преобразования растровых столбцов и строк в проекционные координаты. Например:

import rasterio

# The best way to open a raster with rasterio is through the context manager
# so that it closes automatically

with rasterio.open(path_to_raster) as source:

    data = source.read(1) # Read raster band 1 as a numpy array
    affine = source.affine

# ... do some work with scikit-image and get an array of local maxima locations
# e.g.
# maxima = numpy.array([[0, 0], [1, 1], [2, 2]])
# Also note that convention in a numy array for a 2d array is rows (y), columns (x)

for point in maxima: #Loop over each pair of coordinates
    column = point[1]
    row = point[0]
    x, y = affine * (column, row)
    print x, y

# Or you can do it all at once:

columns = maxima[:, 1]
rows = maxima[:, 0]

xs, ys = affine * (columns, rows)

И оттуда вы можете записать свои результаты в текстовый файл так, как вам нравится (я бы посоветовал взглянуть на встроенный csvмодуль, например).

Беглый взгляд на matplotlib, я бы сказал, что вам нужно изменить масштаб оси после импорта.

Пожалуйста, попробуйте следующий фрагмент кода. Это может быть использовано для считывания данных изображения из растра и записи обработанных данных в растр (файл .geotiff).

from PIL import Image,ImageOps
import numpy as np
from osgeo import gdal
#from osgeo import gdal_array
#from osgeo import osr
#from osgeo.gdalconst import *
#import matplotlib.pylab as plt

#from PIL import Image, ImageOps
#import gdal
#from PIL import Image
gdal.AllRegister()

################## Read Raster #################
inRaster='C:\python\Results\Database\Risat1CRS\CRS_LEVEL2_GEOTIFF\scene_HH\imagery_HH.tif'

inDS=gdal.Open(inRaster,1)
geoTransform = inDS.GetGeoTransform()
band=inDS.GetRasterBand(1)
datatype=band.DataType
proj = inDS.GetProjection()
rows = inDS.RasterYSize
cols=inDS.RasterXSize
data=band.ReadAsArray(0,0,cols,rows)#extraction of data to be processed#
############write raster##########
driver=inDS.GetDriver()
outRaster='C:\\python\\Results\\Database\\Temporary data base\\clipped_26July2017\\demo11.tif'
outDS = driver.Create(outRaster, cols,rows, 1,datatype)
geoTransform = inDS.GetGeoTransform()
outDS.SetGeoTransform(geoTransform)
proj = inDS.GetProjection()
outDS.SetProjection(proj)
outBand = outDS.GetRasterBand(1)
outBand.WriteArray(data1,0,0)
#data is the output array to written in tiff file
outDS=None 
im2=Image.open('C:\\python\\Results\\Database\\Temporary data base\\clipped_26July2017\\demo11.tif');
im2.show()