Συζητήσεις > Μαθήματα και Εξετάσεις
Φωτοερμηνεία Τηλεπισκόπηση Εαρινό 2019 - Άσκηση 3
andronis:
Θεωρώ ότι επεξηγήθηκε πλήρως στην άσκηση αλλά και στη θεωρία.
Σκέψου τι σημαίνει για μια εικόνα η ελάχιστη και η μέγιστη τιμή και τι ο μέσος όρος των τιμών της.
Πότε φαίνεται μαύρο ένα pixel, πότε φαίνεται λευκό και πότε γκρι.
Πότε φαίνεται φωτεινό και πότε σκοτεινό ;
Επίσης σκέψου τι δείχνει σε ένα ιστόγραμμα το εύρος των τιμών της εικόνας. Τι σημαίνει μεγάλο εύρος τιμών και τι μικρό ;
--- Quote from: mariat on 01 Μαρ 2019, 14:52 ---καλησπερα
θα ηθελα να ρωτησω για την φωτεινοτητα και αντιθεση των καναλιων πως θα χρησιμοποιησουμε τον μεσο ορο των τιμων ως αιτιολογηση?
αρκει μονο ο μεσος ορος ως αιτιολογηση?
ευχαριστω
--- End quote ---
Dimitris812:
Επειδη εχω μακ δυσκολευομαι να τρεξω τις εντολες της δευτερης ασκησης στο τερματικο οπως πχ την gdalinfo που την χρειαζομαστε και σε αυτη την ασκηση τι μπορω να κανω ?
chiossif:
--- Quote from: Dimitris812 on 03 Μαρ 2019, 20:31 ---Επειδη εχω μακ δυσκολευομαι να τρεξω τις εντολες της δευτερης ασκησης στο τερματικο οπως πχ την gdalinfo που την χρειαζομαστε και σε αυτη την ασκηση τι μπορω να κανω ?
--- End quote ---
Καλησπέρα :-)
Αν είναι φορητός ο Μάκης φέρτον στις ώρες γραφείου να το δούμε μαζί. Αν όχι έλα έτσι κι αλλιώς στο γραφείο διότι το γκνου/λίνουξ είναι 99% ίδιο με το γιούνιξ του τερματικού του Μάκη. Οπότε θα το δεις ... Ήδη ένας άλλος φοιτητής με Μάκη έλυσε με επίσκεψη στο γραφείο τα προβλήματά του :-)
Λέφτερα,
Ch iossif
chiossif:
Κυρίες Δεσποσύνες και Κύριοι,
καλή σας ώρα :-)
Δείτε
σε C:
--- Code: ---/*
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
*/
--- End code ---
...
chiossif:
...
και σε Python:
--- Code: ---#
# 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
--- End code ---
εκτεταμένα προγράμματα και με υπολογισμό στατιστικών πίξελ προς πίξελ. Αλήθεια τώρα βλέπετε γιατί διαφέρουν τα στατιστικά ; :-)
Λέφτερα,
Ch Iossif
υγ
Όσες και όσοι εκτελέσετε και τα δύο προγράμματα παρατηρείστε την διαφορά ταχύτητας εκτέλεσης μεταξύ της C και της Python :-)
Πλοήγηση
[0] Λίστα μηνυμάτων
Go to full version