Przetwarzanie obrazu przy użyciu Pythona, GDAL i Scikit-Image

Mam problem z przetwarzaniem i mam nadzieję, że uda mi się tutaj rozwiązać.

Pracuję z teledetekcją stosowaną w leśnictwie, szczególnie z danymi LiDAR. Pomysł polega na użyciu obrazu Scikit do wykrywania wierzchołków drzew. Ponieważ jestem nowy w Pythonie, uznałem wielki osobisty triumf, aby wykonać następujące czynności:

1. Zaimportuj CHM (z matplotlib);
2. Uruchom filtr gaussowski (z pakietem scikit-image);
3. Uruchom filtr maksima (z pakietem scikit-image);
4. Uruchom peak_local_max (z pakietem scikit-image);
5. Pokaż CHM z lokalnymi maksimami (z matplotlib);

Teraz mój problem. Podczas importowania za pomocą programu matplot obraz traci współrzędne geograficzne. Współrzędne, które mam, są tylko podstawowymi współrzędnymi obrazu (tj. 250 312). Potrzebuję uzyskać wartość piksela poniżej lokalnych maksymalnych punktów na obrazie (czerwone kropki na obrazie). Tutaj na forum widziałem jednego faceta pytającego o to samo (uzyskanie wartości piksela rastra GDAL pod punktem OGR bez NumPy? ), Ale on już miał te punkty w pliku kształtu. W moim przypadku punkty zostały obliczone za pomocą scikit-image (jest to tablica ze współrzędnymi każdego wierzchołka drzewa). Więc nie mam pliku kształtu.

Podsumowując, na końcu chcę to plik txt ze współrzędnymi każdej lokalnej maksima we współrzędnych geograficznych, na przykład:

525412 62980123 1150 ...

Po pierwsze, witamy na stronie!

Tablice numpy nie mają wbudowanych w układ koncepcji współrzędnych. W przypadku rastra 2D są one indeksowane według kolumn i wierszy.

Uwaga: Zakładam, że czytasz format rastrowy obsługiwany przez GDAL .

W Pythonie najlepszym sposobem importowania przestrzennych danych rastrowych jest rasteriopakiet. Surowe dane importowane przez rasterio nadal są tablicą liczb liczbowych bez dostępu do układów współrzędnych, ale rasterio daje również dostęp do metody afinicznej w tablicy źródłowej, której można użyć do przekształcenia kolumn i wierszy rastrowych w rzutowane współrzędne. Na przykład:

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)

I stamtąd możesz zapisać nasze wyniki do pliku tekstowego, jak chcesz (sugeruję na przykład zajrzeć do wbudowanego csvmodułu ).

Z szybkiego spojrzenia na matplotlib powiedziałbym, że po imporcie musisz zmienić skale osi .

Spróbuj użyć następującego fragmentu kodu. Można to wykorzystać do odczytu danych obrazu z rastra i zapisu przetworzonych danych do rastra (plik .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()