Γεια και χαρά σε όλες και όλους :-)
Το πρόγραμμα RasterKMeans_multiband.py εκτελεί μη επιβλεπόμενη ταξινόμηση:
#
# DRSEx. 7 - Remote Sensing Lab. - SRSE NTUA
#
# Image unsupervised classification (multiband)
#
# License: https://www.gnu.org/licenses/gpl.txt
#
import numpy, sys, gdal, gdalconst
from scipy.cluster.vq import *
# register all of the GDAL drivers
gdal.AllRegister()
# open the image
imagefilename='Cropped_A' ### Cropped_A.tif multiband image file name
inDs = gdal.Open(imagefilename+'.tif', gdalconst.GA_ReadOnly)
if inDs is None:
print('Could not open '+imagefilename+'.tif')
sys.exit(1)
print('Image file name = '+imagefilename+'.tif')
# get image size
rows = inDs.RasterYSize
cols = inDs.RasterXSize
bands = inDs.RasterCount
print('Rows = '+str(rows)+' Cols = '+str(cols)+' Bands = ',str(bands))
# read image into Numpy Array
raster = inDs.ReadAsArray().astype(numpy.float32)
print('Shape of input image : '+str(raster.shape))
# flatten image to get line of values [ http://users.ntua.gr/chiossif/Free_As_Freedom_Software/BIL_BIP_BSQ.pdf ]
flatraster = numpy.zeros((rows*cols, bands))
for r in range(rows):
for c in range(cols):
for b in range(bands):
flatraster[r*cols+c,b]=raster[b,r,c]
print('Shape of K-Means algorithm input : '+str(flatraster.shape))
# run k-means algorithm with multiple clusters
for i in range(8,17,4):
kmeansfilename=imagefilename+'_k-means_'+str(i)+'_clusters.tif' # Output image file name
driver = inDs.GetDriver()
outDs = driver.Create(kmeansfilename, cols, rows, 1, gdalconst.GDT_Byte)
if outDs is None:
print('Could not create '+ kmeansfilename)
sys.exit(1)
outBand = outDs.GetRasterBand(1)
print("Calculating k-means with ", i, " clusters.")
# This scipy code classifies with k-mean unsupervised classifier in two steps:
# centroid estimation and image classification.
# centroids have cluster size and
# code has same length as flattened image
# and defines which class the value corresponds to
centroids, variance = kmeans(flatraster, i)
print('Shape of K-Means centroid array : '+str(centroids.shape))
code, distance = vq(flatraster, centroids)
print('Shape of K-Means classified array : '+str(code.shape))
#Since code contains the classified values, reshape back to image dimensions
codeim = code.reshape(raster.shape[1], raster.shape[2])
print('Shape of K-Means classified image : '+str(codeim.shape))
# write data to disk
outBand.WriteArray(codeim)
outBand.FlushCache()
# georeference the image and set the projection
outDs.SetGeoTransform(inDs.GetGeoTransform())
outDs.SetProjection(inDs.GetProjection())
outDs = None
print('Classified (KMeans) image file name = '+kmeansfilename+' DONE !')
inDs = None
# Display QGIS message
print(
"""
# Run inside QGIS at Plugins->Open python console >>> prompt
import os # make os functions avaliable here
os.chdir('C:\DRS\Ex7') # change to working folder
execfile('RasterKMeans_multiband.py') # run your code Python2 QGIS 2.x OR
exec(open('RasterKMeans_multiband.py').read()) # run your code Python3 QGIS 3.x
"""
)
ΠΡΟΣΟΧΗ
Το όνομα της εικόνας το αλλάζετε στην γραμμή :
imagefilename='Cropped_A' ### Cropped_A.tif multiband image file name
ενώ στην γραμμή
for i in range(8,17,4):ρυθμίζετε πόσες κλάσεις θέλετε. Εδώ από 8 μέχρι 17 με βήμα 4 δηλαδή 8 12 και 16.
Οδηγίες:
1. Φτιάχνουμε τον φάκελο της άσκησης έστω C:\DRS\Ex7 και βάζουμε μέσα το πρόγραμμα RasterKMeans_multiband.py και την εικόνα μας έστω Cropped_A.tif
2. Ανοίγουμε το QGIS και τρέχουμε με την σειρά τις εντολές στην επιλογή Plugins->Open python console :
2.1. import os
2.2. os.chdir('C:\DRS\Ex7')
2.3. exec(open('RasterKMeans_multiband.py').read())
3. Το πρόγραμμα τρέχει (αργεί αρκετά και ζεσταίνει το πισί) αλλά φτιάχνει 3 μη επιβλεπόμενες ταξινομήσεις με 8 12 και 16 τάξεις με ονόματα Cropped_A_k-means_8_clusters.tif, Cropped_A_k-means_12_clusters.tif και Cropped_A_k-means_12_clusters.tif αντίστοιχα.
4. ΑΦΟΥ τελειώσει το πρόγραμμα (το QGIS ΔΕΝ θα αντιδρά νωρίτερα) ανοίγουμε τις τρεις ταξινομήσεις.
5. Για να συνεχίσουμε την άσκηση με τον χρωματισμό των κατηγοριών και την ονοματοδοσία τους προχωράμε κανονικά: δεξί κλικ στο Layers, Properties, Symbology, Επιλέγουμε Palleted /Unique values στο Render type, πατάμε Classify, βάζουμε κόκκινο το 1ο και τα λοιπά όπως δείξαμε...
Ελπίζω να βοήθησα και με Πάιθον :-)
Λέφτερα,
Ch Iossif
