topografoi.com



Author Topic: Φωτοερμηνεία Τηλεπισκόπηση Εαρινό 2019 - Άσκηση 3  (Read 8863 times)

chiossif

  • Posts: 334
Κυρίες, Δεσποσύνες και Κύριοι,

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

Σήμερα 27/2 ολοκληρώθηκε η 3η Άσκηση του τρέχοντος εξαμήνου στην ΦΤ. Έτσι, για να βοηθήσουμε την διαδικασία επίλυσης αποριών, ανοίγουμε αυτό το θέμα στο οποίο μπορείτε να υποβάλετε τις ερωτήσεις σας και να διαβάζετε τις ερωτήσεις και τις απαντήσεις των άλλων.

Μην ξεχνάτε ΠΡΙΝ ρωτήσετε να διαβάζετε κατά σειρά: τις οδηγίες για το μάθημα και τις ασκήσεις του, τις σημειώσεις σας και ότι έχει ήδη ερωτηθεί εδώ ή αναρτηθεί εδώ.

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

Λέφτερα,
Ch Iossif

chiossif

  • Posts: 334
Καλή ώρα σε όλες και όλους :-)

Υποδείξεις στο 3ο ερώτημα :

-> Τι παραδίδουμε (με σειρά σπουδαιότητας):

Τα αποτελέσματα εκτέλεσης του προγράμματος και για τις δύο εικόνες,
σύγκριση και αιτιολόγηση των τιμών οι οποίες προκύπτουν εδώ και από άλλα ερωτήματα νωρίτερα ή εντολές (gdalinfo),
οδηγίες για την εκτέλεση του προγράμματος,
μια περιγραφή των αλλαγών που κάναμε στο πρόγραμμα,
μια σύντομη αλλά πλήρη περιγραφή του κώδικα του προγράμματος.

-> Επιλέγουμε ΜΙΑ έκδοση, μία γλώσσα προγραμματισμού C ή C++ ή Python. Συνιστάται η Python εξαιτίας της ευκολίας αλλά και της συμβατότητάς της με το QGIS.

Αναλυτικότερα οι αλλαγές οι οποίες απαιτούνται:

-> Αλλαγή για το όνομα της εικόνας:
στην γραμμή 18 για τις C, C++ στην γραμμή 13 για την Python αλλάζουμε το πλήρες όνομα της εικόνας μας από το "/home/chiossif" σε "C:\\FT\\ASK3\\kommeni1.tif" για παράδειγμα σύμφωνα με τις οδηγίες που έχουμε δώσει.

-> Αλλαγή για το πλήθος των καναλιών:
στην γραμμή 54 στην C, στην γραμμή 52 στην C++ στην γραμμή 31 στην Python αλλάζουμε το εύρος των καναλιών. Δεν χρειάζεται αλλαγή η πρώτη τιμή (1) η οποία περιλαμβάνεται πάντα ενώ η 2η τιμή θέλει μείωση κατά 1 σε όλες τις γλώσσες. Προσέξτε ότι η range στην Python είναι ανοιχτή στο δεύτερο άκρο :-)

-> Αλλαγή για τα αρχεία επικεφαλίδες στις C C++ μόνο:
είναι πιθανή λόγω διαφορετικών εγκαταστάσεων qgis gdal και μεταγλωτιστών C C++ στον υπολογιστή σας διότι εξαρτάται από όλα αυτά στα παράθυρα της μικρομαλακής του Βασιλάκη του Πύλη ;-)

C C++ : ακολουθείστε τις οδηγίες / διαδικασία μεταγλωτίσεων όπως κάνατε στα ανάλογα μαθήματα προηγούμενων εξαμήνων.

Python : (επιπλέον οδηγίες με σειρά ενδιαφέροντος)
-> χρησιμοποιείστε το πρόγραμμα όπως προτείνεται στις τελευταίες γραμμές του κώδιμα δηλαδή μέσα από το QGIS. Στην 1η εντολή φορτώνουμε την βιβλιοθήκη στην 2η αλλάζουμε στον φάκελο εργασίας μας πχ C:\\FT\\ASK3\\kommeni1.tif και στην 3η (για το QGIS2) ή 4η (για το QGIS3) η εκτέλεση του προγράμματος. Τα αποτελέσματα βγαίνουν αυτόματα στο ίδιο μέρος. Τα βλέπετε με κύλιση προς τα πάνω.
-> διαβάστε την εισαγωγή στην Python
-> εγκαταστείστε το λογισμικό Sololearn Python στο κινητό σας για εύκολη εκμάθηση της γλώσσας
-> αναζητείστε στο διαδίκτυο την πληθώρα εγχειριδίων μαθημάτων τα οποία δίνονται ελεύθερα, ανοιχτά ή απλά δωρεάν.

Στην διάθεσή σας για περισσότερες οδηγίες ή υλικό αυτοδιδασκαλίας προγραμματισμού εντός του πλαισίου δραστηριοτήτων της κοινότητας.

Λέφτερα,
Ch Iossif

mariat

  • Posts: 5
καλησπερα
θα ηθελα να ρωτησω για τα καναλια red nir και swir που πρεπει να εμφανισουμε
νομιζω μας ειχατε πει οτι το κοκκινο καναλι ειναι το 3ο της εικονας μας αλλα με τη σειρα RGB και απο το ιστογραμμα δεν θα επρεπε να ειναι το 1ο?
ευχαριστω πολυ

andronis

  • Posts: 20
@mariat :
Δες τα κανάλια του Landsat 8 εδώ.
Όπως έκοψες την εικόνα σου, κράτησες τα κανάλια 2 έως 7 της αρχικής που κατέβασες. Οπότε η κομμένη εικόνα σου (εικόνες σου)  περιέχει τα κανάλια που περιγράφονται στο pdf που βλέπεις. Το κανάλι RED λοιπόν είναι το τρίτο ανάλι της κομμένης εικόνας σου.
Για τη  θεωρία των χρωμάτων και το χρωματικό μοντέλο RGB είπατε στη θεωρία. Δεν πρέπει να το συγχέεις με τα κανάλια του δορυφόρου.  Βοηθητικά δες και εδώ

mariat

  • Posts: 5
καλησπερα
θα ηθελα να ρωτησω για την φωτεινοτητα και αντιθεση των καναλιων πως θα χρησιμοποιησουμε τον μεσο ορο των τιμων ως αιτιολογηση?
αρκει μονο ο μεσος ορος ως αιτιολογηση?
ευχαριστω

andronis

  • Posts: 20
Θεωρώ ότι επεξηγήθηκε πλήρως στην άσκηση αλλά και στη θεωρία.
Σκέψου τι σημαίνει για  μια εικόνα η ελάχιστη και η μέγιστη τιμή και τι ο μέσος όρος των τιμών της.
Πότε φαίνεται μαύρο ένα pixel, πότε φαίνεται λευκό και πότε γκρι.
Πότε φαίνεται φωτεινό και πότε σκοτεινό ;
Επίσης σκέψου τι δείχνει σε ένα ιστόγραμμα το εύρος των τιμών της εικόνας. Τι σημαίνει μεγάλο εύρος τιμών και τι μικρό ;


καλησπερα
θα ηθελα να ρωτησω για την φωτεινοτητα και αντιθεση των καναλιων πως θα χρησιμοποιησουμε τον μεσο ορο των τιμων ως αιτιολογηση?
αρκει μονο ο μεσος ορος ως αιτιολογηση?
ευχαριστω

Dimitris812

  • Posts: 2
Επειδη εχω μακ δυσκολευομαι να τρεξω τις εντολες της δευτερης ασκησης στο τερματικο οπως πχ την gdalinfo που την χρειαζομαστε και σε αυτη την ασκηση τι μπορω να κανω ?

chiossif

  • Posts: 334
Επειδη εχω μακ δυσκολευομαι να τρεξω τις εντολες της δευτερης ασκησης στο τερματικο οπως πχ την gdalinfo που την χρειαζομαστε και σε αυτη την ασκηση τι μπορω να κανω ?

Καλησπέρα :-)

Αν είναι φορητός ο Μάκης φέρτον στις ώρες γραφείου να το δούμε μαζί. Αν όχι έλα έτσι κι αλλιώς στο γραφείο διότι το γκνου/λίνουξ είναι 99% ίδιο με το γιούνιξ του τερματικού του Μάκη. Οπότε θα το δεις ... Ήδη ένας άλλος φοιτητής με Μάκη έλυσε με επίσκεψη στο γραφείο τα προβλήματά του :-)

Λέφτερα,
Ch iossif

chiossif

  • Posts: 334
Κυρίες Δεσποσύνες και Κύριοι,

καλή σας ώρα :-)

Δείτε

σε C:

Code: [Select]
/*
    RSEx3 - Remote Sensing Lab. - SRSE NTUA
    
    GDAL Getting Dataset Information
    
    Reference: https://www.gdal.org/gdal_tutorial.html
    
    License: https://www.gnu.org/licenses/gpl.txt
*/

#include <stdio.h> /* for printf */
#include <stdlib.h> /* for exit */
#include "gdal/gdal.h"
#include "gdal/cpl_conv.h" /* for CPLMalloc() */

int main(void){
/* Enter your cropped.tif image full path filename here  */
    char *pszFilename="/home/chiossif/Data/LC08_L1TP_184033_20180622_20180703_01_T1/clipped.tif";
    
    GDALDatasetH  hDataset;
    GDALAllRegister();
    hDataset = GDALOpen( pszFilename, GA_ReadOnly );
    if( hDataset == NULL ) {
printf("Unable to open %s file\n", pszFilename);
        exit(1);
    }

    GDALDriverH   hDriver;
    double        adfGeoTransform[6];
    hDriver = GDALGetDatasetDriver( hDataset );
    printf( "Driver: %s/%s\n",
            GDALGetDriverShortName( hDriver ),
            GDALGetDriverLongName( hDriver ) );
    printf( "Size is %dx%dx%d\n",
            GDALGetRasterXSize( hDataset ),
            GDALGetRasterYSize( hDataset ),
            GDALGetRasterCount( hDataset ) );
    if( GDALGetProjectionRef( hDataset ) != NULL )
        printf( "Projection is `%s'\n",
            GDALGetProjectionRef( hDataset ) );
    if( GDALGetGeoTransform( hDataset, adfGeoTransform ) == CE_None ){
        printf( "Origin = (%.6f,%.6f)\n",
            adfGeoTransform[0], adfGeoTransform[3] );
        printf( "Pixel Size = (%.6f,%.6f)\n",
            adfGeoTransform[1], adfGeoTransform[5] );
    }

    GDALRasterBandH hBand;
    int             nBlockXSize, nBlockYSize;
    int             bGotMin, bGotMax;
    double          adfMinMax[2];
    double          myMin, myMax, myMean, myStdD, mySum, mySum2;
    float *pafScanline, tmpfloat;
    int             bands, bandno, myNum;
    register int i, r;
    
    bands = GDALGetRasterCount( hDataset );
    for (bandno=1;bandno<=bands;bandno++){
     printf("Band=%d\n",bandno);
     hBand = GDALGetRasterBand( hDataset, bandno );
     GDALGetBlockSize( hBand, &nBlockXSize, &nBlockYSize );
     printf( "Block=%dx%d Type=%s, ColorInterp=%s\n",
            nBlockXSize, nBlockYSize,
            GDALGetDataTypeName(GDALGetRasterDataType(hBand)),
        GDALGetColorInterpretationName( GDALGetRasterColorInterpretation(hBand)) );
        adfMinMax[0] = GDALGetRasterMinimum( hBand, &bGotMin );
        adfMinMax[1] = GDALGetRasterMaximum( hBand, &bGotMax );
        if( ! (bGotMin && bGotMax) )
            GDALComputeRasterMinMax( hBand, TRUE, adfMinMax );
        printf( "Min=%.3f, Max=%.3f\n", adfMinMax[0], adfMinMax[1] );
        if( GDALGetOverviewCount(hBand) > 0 )
            printf( "Band has %d overviews.\n", GDALGetOverviewCount(hBand));
        if( GDALGetRasterColorTable( hBand ) != NULL )
            printf( "Band has a color table with %d entries.\n",
        GDALGetColorEntryCount( GDALGetRasterColorTable( hBand ) ) );
        
        pafScanline = (float *) CPLMalloc(sizeof(float)*nBlockXSize);
        if (pafScanline!=NULL){
myMin = 65536.0;
myMax = -65537.0;
mySum  = 0.0;
mySum2 = 0.0;
myNum = 0;
for (r=0;r < GDALGetRasterYSize( hDataset );r++){
GDALRasterIO( hBand, GF_Read, 0, r, nBlockXSize, 1, pafScanline, nBlockXSize, 1, GDT_Float32, 0, 0 );
for (i=0;i<nBlockXSize;i++){
myNum++;
tmpfloat = pafScanline[i];
mySum  += tmpfloat;
mySum2 += tmpfloat * tmpfloat;
if (tmpfloat > myMax)
 myMax = tmpfloat;
if (tmpfloat < myMin)
 myMin = tmpfloat;
}
}
   myMean = mySum  / myNum;
   myStdD = sqrt(mySum2 / myNum - myMean * myMean);
          printf( "Min=%.3lf, Max=%.3lf, Mean=%.3lf, StdDev=%.3lf\n", myMin, myMax, myMean, myStdD);
 CPLFree(pafScanline);
 }
        }

    return 0;
}
/*
Driver: GTiff/GeoTIFF
Size is 2840x1541x7
Projection is `PROJCS["WGS 84 / UTM zone 34N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",21],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32634"]]'
Origin = (498255.000000,4282875.000000)
Pixel Size = (30.000000,-30.000000)
Band=1
Block=2840x1 Type=UInt16, ColorInterp=Gray
Min=9260.000, Max=48540.000
Min=9260.000, Max=48540.000, Mean=11449.392, StdDev=2289.782
Band=2
Block=2840x1 Type=UInt16, ColorInterp=Undefined
Min=8217.000, Max=51117.000
Min=8217.000, Max=51117.000, Mean=10476.116, StdDev=2487.918
Band=3
Block=2840x1 Type=UInt16, ColorInterp=Undefined
Min=6982.000, Max=52368.000
Min=6982.000, Max=52368.000, Mean=9604.907, StdDev=2668.937
Band=4
Block=2840x1 Type=UInt16, ColorInterp=Undefined
Min=6159.000, Max=55830.000
Min=6159.000, Max=55830.000, Mean=8906.102, StdDev=3117.561
Band=5
Block=2840x1 Type=UInt16, ColorInterp=Undefined
Min=5497.000, Max=53726.000
Min=5052.000, Max=63093.000, Mean=15602.425, StdDev=6334.375
Band=6
Block=2840x1 Type=UInt16, ColorInterp=Undefined
Min=5097.000, Max=60158.000
Min=4714.000, Max=63849.000, Mean=12332.979, StdDev=5270.812
Band=7
Block=2840x1 Type=UInt16, ColorInterp=Undefined
Min=4978.000, Max=53259.000
Min=4892.000, Max=65535.000, Mean=9391.264, StdDev=3782.803
*/

...

chiossif

  • Posts: 334
...
και σε Python:

Code: [Select]
#
# RSEx3 - Remote Sensing Lab. - SRSE NTUA
#
# GDAL Getting Dataset Information
#
# Reference: https://www.gdal.org/gdal_tutorial.html
#
# License: https://www.gnu.org/licenses/gpl.txt
#
from osgeo import gdal
from math import sqrt
import struct

# Enter your cropped.tif image full path filename here
filename="/home/chiossif/Data/LC08_L1TP_184033_20180622_20180703_01_T1/clipped.tif"

dataset = gdal.Open(filename, gdal.GA_ReadOnly)
if not dataset:
    print("Unable to open ",filename," file")
    sys.exit(1)

print("Driver: {}/{}".format(dataset.GetDriver().ShortName,
                             dataset.GetDriver().LongName))
print("Size is {} x {} x {}".format(dataset.RasterXSize,
                                    dataset.RasterYSize,
                                    dataset.RasterCount))
print("Projection is {}".format(dataset.GetProjection()))
geotransform = dataset.GetGeoTransform()
if geotransform:
    print("Origin = ({}, {})".format(geotransform[0], geotransform[3]))
    print("Pixel Size = ({}, {})".format(geotransform[1], geotransform[5]))

bands = dataset.RasterCount
for bandno in range(1,bands+1):
    print("Band={}".format(bandno))
    band = dataset.GetRasterBand(bandno)
    print("Band Type={}".format(gdal.GetDataTypeName(band.DataType)))

    min = band.GetMinimum()
    max = band.GetMaximum()
    if not min or not max:
        (min,max) = band.ComputeRasterMinMax(True)
    print("Min={:.3f}, Max={:.3f}".format(min,max))
      
    if band.GetOverviewCount() > 0:
        print("Band has {} overviews".format(band.GetOverviewCount()))
      
    if band.GetRasterColorTable():
        print("Band has a color table with {} entries".format(band.GetRasterColorTable().GetCount()))
    
    myMin =  65536
    myMax = -65537
    mySum = mySum2 = 0.0
    myNum = 0
    for y in range(0,dataset.RasterYSize):
        scanline = band.ReadRaster(xoff=0, yoff=y, xsize=band.XSize, ysize=1, buf_xsize=band.XSize, buf_ysize=1, buf_type=gdal.GDT_Float32)
        tuple_of_floats = struct.unpack('f' * band.XSize, scanline)
        for x in tuple_of_floats:
            myNum  += 1
            mySum  += x
            mySum2 += x * x
            if x > myMax:
                myMax = x
            if x < myMin:
                myMin = x

    myMean = mySum /myNum
    myStdD = sqrt(mySum2 / myNum - myMean * myMean)
    print("Min={:.3f}, Max={:.3f},  Mean={:.3f}, StdDev={:.3f}".format(myMin,myMax,myMean,myStdD))
  

# Display QGIS message
print(
"""
# Run inside QGIS at Plugins->Open python console >>> prompt
import os # make os functions avaliable here
os.chdir('C:\Tilepiskopisi\Ex3') # change to working folder
execfile('gdal_tutorial.py') # run your code Python2 QGIS 2.x OR
exec(open('gdal_tutorial.py').read()) # run your code Python3 QGIS 3.x
"""
)
# Driver: GTiff/GeoTIFF
# Size is 2840 x 1541 x 7
# Projection is PROJCS["WGS 84 / UTM zone 34N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",21],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32634"]]
# Origin = (498255.0, 4282875.0)
# Pixel Size = (30.0, -30.0)
# Band=1
# Band Type=UInt16
# Min=9260.000, Max=48540.000
# Min=9260.000, Max=48540.000,  Mean=11449.392, StdDev=2289.782
# Band=2
# Band Type=UInt16
# Min=8217.000, Max=51117.000
# Min=8217.000, Max=51117.000,  Mean=10476.116, StdDev=2487.918
# Band=3
# Band Type=UInt16
# Min=6982.000, Max=52368.000
# Min=6982.000, Max=52368.000,  Mean=9604.907, StdDev=2668.937
# Band=4
# Band Type=UInt16
# Min=6159.000, Max=55830.000
# Min=6159.000, Max=55830.000,  Mean=8906.102, StdDev=3117.561
# Band=5
# Band Type=UInt16
# Min=5497.000, Max=53726.000
# Min=5052.000, Max=63093.000,  Mean=15602.425, StdDev=6334.375
# Band=6
# Band Type=UInt16
# Min=5097.000, Max=60158.000
# Min=4714.000, Max=63849.000,  Mean=12332.979, StdDev=5270.812
# Band=7
# Band Type=UInt16
# Min=4978.000, Max=53259.000
# Min=4892.000, Max=65535.000,  Mean=9391.264, StdDev=3782.803

# # Run inside QGIS at Plugins->Open python console >>> prompt
# import os # make os functions avaliable here
# os.chdir('C:\Tilepiskopisi\Ex3') # change to working folder
# execfile('gdal_tutorial.py') # run your code Python2 QGIS 2.x OR
# exec(open('gdal_tutorial.py').read()) # run your code Python3 QGIS 3.x

εκτεταμένα προγράμματα και με υπολογισμό στατιστικών πίξελ προς πίξελ. Αλήθεια τώρα βλέπετε γιατί διαφέρουν τα στατιστικά ; :-)

Λέφτερα,
Ch Iossif

υγ
Όσες και όσοι εκτελέσετε και τα δύο προγράμματα παρατηρείστε την διαφορά ταχύτητας εκτέλεσης μεταξύ της C και της Python :-)


 

Copyright © topografoi.com