topografoi.com

Συζητήσεις => Καταγγελίες => Μαθήματα και Εξετάσεις => Topic started by: chiossif on 02 Δεκ 2019, 20:27

Title: Ψηφιακή Τηλεπισκόπηση Σ.Α.Τ.Μ. Ε.Μ.Π. Χειμερινό 2019 'Ασκηση 7
Post by: chiossif on 02 Δεκ 2019, 20:27
Κυρίες, Δεσποσύνες και Κύριοι,

καλησπέρα σας :-)

Σήμερα 2/12 ολοκληρώθηκε η 7η Άσκηση (http://mycourses.ntua.gr/courses/SURVEY1011/projects/DRS_LAB_07.pdf) του τρέχοντος εξαμήνου στην ΨΤ. Έτσι, για να βοηθήσουμε την διαδικασία επίλυσης αποριών, ανοίγουμε αυτό το θέμα στο οποίο μπορείτε να υποβάλετε τις ερωτήσεις σας και να διαβάζετε τις ερωτήσεις και τις απαντήσεις των άλλων.

Μην ξεχνάτε ΠΡΙΝ ρωτήσετε να διαβάζετε κατά σειρά: την γενική περιγραφή για το μάθημα (http://mycourses.ntua.gr/courses/SURVEY1011/document/DRS_Intro_2019.pdf) και τις οδηγίες για τις ασκήσεις (http://mycourses.ntua.gr/courses/SURVEY1011/document/DRS_LABintro_2019.pdf) του, τις σημειώσεις σας και ότι έχει ήδη ερωτηθεί εδώ.

Καλή και γόνιμη μελέτη :-)

Λέφτερα,
Ch Iossif
Title: Απ: Ψηφιακή Τηλεπισκόπηση Σ.Α.Τ.Μ. Ε.Μ.Π. Χειμερινό 2019 'Ασκηση 7
Post by: chiossif on 03 Δεκ 2019, 10:40
Γεια και χαρά σε όλες και όλους :-)

Το πρόγραμμα RasterKMeans_multiband.py εκτελεί μη επιβλεπόμενη ταξινόμηση:

Code: [Select]
#
# 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
"""
)

ΠΡΟΣΟΧΗ
Το όνομα της εικόνας το αλλάζετε στην γραμμή :
Code: [Select]
imagefilename='Cropped_A' ### Cropped_A.tif multiband image file name
ενώ στην γραμμή
Code: [Select]
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

Title: Απ: Ψηφιακή Τηλεπισκόπηση Σ.Α.Τ.Μ. Ε.Μ.Π. Χειμερινό 2019 'Ασκηση 7
Post by: chiossif on 03 Δεκ 2019, 16:28
Γεια και χαρά σε όλες και όλους :-)

Μια και ανοίξαμε τον χορό με κώδικα σε πύθωνα λάβετε και το ConfusionMatrix.py (http://users.ntua.gr/chiossif/Free_As_Freedom_Software/ConfusionMatrix.py) .

Οδηγίες:

1. Ανοίγουμε τον κώδικα με έναν διορθωτή κειμένου (https://notepad-plus-plus.org) και βάζουμε στις γραμμές:

Code: [Select]
### Enter your classified image full path filename here
class_filename="Cropped_A_k-means_12_clusters.tif"
### Enter your truth image full path filename here
truth_filename="Cropped_A_k-means_8_clusters.tif"
### Enter your confusion matrix full path filename here
conf_filename="Cropped_A_k-means_8_12.csv"

το όνομα των αρχείων με την ταξινομημένη, την αναφοράς και το αρχείο csv για τον πίνακα σύγχυσης.

2. Αποθηκεύουμε το αρχείο ConfusionMatrix.py στον φάκελο C:\DRS\Ex7 όπου έχουμε και τις εικόνες

3. Ανοίγουμε το QGIS και τρέχουμε με την σειρά τις εντολές στην επιλογή Plugins->Open python console :
3.1. import os
3.2. os.chdir('C:\DRS\Ex7')
3.3. exec(open('ConfusionMatrix.py').read())

4. Διαβάζουμε προς τα πίσω τα αποτελέσματα του προγράμματος με κύλιση στην οθόνη και ανοίγουμε το αρχείο csv σε ένα πρόγραμμα διαχείρισης λογιστικών φύλλων (https://www.libreoffice.org/discover/calc/) ώστε να κάνουμε εκεί ότι θέλουμε.

The sky's the limit.

Λέφτερα,
Ch Iossif

υγ.
ΠΡΟΣΟΧΗ!
Στο παράδειγμα συγκρίνω δύο μη επιβλεπόμενες ταξινομήσεις και με μεταβλητό πλήθος κατηγοριών οπότε οι υπολογισμοί ακριβειών κτλ δεν ισχύουν ως τέτοιοι.
Το πρόγραμμα αυτό λειτουργεί ΑΠΟΚΛΕΙΣΤΙΚΑ με δεδομένα raster (άρα ΔΕΝ κάνει για την άσκηση 6 άμεσα).
Title: Απ: Ψηφιακή Τηλεπισκόπηση Σ.Α.Τ.Μ. Ε.Μ.Π. Χειμερινό 2019 'Ασκηση 7
Post by: chiossif on 04 Ιαν 2020, 13:19
Γεια και χαρά σε όλες και όλους :-)

Καλή Χρονιά και πάντα με Υγεία :-)

Ο υπολογιστής μας έχει ήδη το saga εγκατεστημένο και μερικές φορές η εκτέλεση από το QGIS είναι πιο δύσκολη. Εδώ θα γράψω τα βήματα για μια εύκολη και γρήγορη εκτέλεση του αλγόριθμου K-means απευθείας από το saga :-)


Ελπίζω να βοήθησα :-)

Καλή και Γόνιμη Χρονιά :-)

Λέφτερα,
Ch Iossif
Title: Απ: Ψηφιακή Τηλεπισκόπηση Σ.Α.Τ.Μ. Ε.Μ.Π. Χειμερινό 2019 'Ασκηση 7
Post by: chiossif on 08 Ιαν 2020, 19:27
Κυρίες, Δεσποσύνες και Κύριοι,

καλησπέρα σας και Καλή Χρονιά :-)

Ερώτηση:
«για την μη επιβλεπόμενη ταξινόμηση με αλγόριθμο Κ means clustering κατέβασα το orfeo toolbox από το λινκ που μας δίνεται στην εκφώνηση της άσκησης. Ύστερα έκανα unzip το αρχείο και το αντέγραψα στο c: . Μετά από το qgis Πήγα στο processing --> options --> providers και δεν βρήκα το orfeo ώστε να κάνω activate.. τι μπορεί να πήγε λάθος;»
Απάντηση:
Τίποτε στο οποίο να φταις εσύ. Το πρόγραμμα orfeo έπαψε να είναι εντελώς συμβατό με το ΚιουΤζιΑιΕς. Προτείνω να χρησιμοποιήσεις τις άλλες λύσεις που επιδείχθηκαν στην άσκηση με GrassGIS (https://grasswiki.osgeo.org/wiki/Image_classification για περισσότερα , εύκολο εντός του ΚιουΤζιΑιΕς, δύσκολο εκτός)  ή SAGA (βλ. και απάντηση εδώ (http://www.topografoi.com/forum/index.php?topic=655.0) εύκολο εντός και εκτός του ΚιουΤζιΑιΕς ;-) ).

Ελπίζω να βοήθησα :-)

Χρόνια Πολλά !

Λέφτερα,
Ch Iossif
Title: Απ: Ψηφιακή Τηλεπισκόπηση Σ.Α.Τ.Μ. Ε.Μ.Π. Χειμερινό 2019 'Ασκηση 7
Post by: chiossif on 09 Ιαν 2020, 01:49
Κυρίες, Δεσποσύνες και Κύριοι,

Καλή Χρονιά και πάντα με Υγεία :-)

Το SNAP (https://step.esa.int/main/toolboxes/snap/) είναι ένα σύγχρονο και εξαιρετικό πρόγραμμα διαχείρισης δεδομένων παρατήρησης της γης. Και ως τέτοιο κάνει και μη επιβλεπόμενη ταξινόμηση με φιλικό τρόπο παρότι ΔΕΝ συνεργάζεται με το Κιου Τζι Άι Ές. Πάμε να δούμε πως:


Με την ολοκλήρωση της εκτέλεσης η ταξινομημένη εικόνα kommeni_kmeans_8_clusters.tif αποθηκεύεται στον δίσκο στον ίδιο φάκελο με την αρχική εικόνα. Έτσι συνεχίζουμε την εκπόνηση της άσκησης στο Κιου Τζι Άι Ές όπως μας ζητά η εκφώνηση.

Καλή και Γόνιμη Χρονιά :-)

Λέφτερα,
Ch Iossif