diff --git a/.gitignore b/.gitignore index f08278d..d672ee6 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1,20 @@ -*.pdf \ No newline at end of file +*.pdf +*.TIF +*.TIF.aux.xml +*.db +*.xml +*.jp2 +*.o +a.out +*.out +data/ +data1/ +data2/ +data3/ +S2A_MSIL2A_20170527T102031_N9999_R065_T33UUU_20191015T100203.SAFE/ +S2A_MSIL2A_20241118T143741_N0511_R096_T18GXU_20241118T201252.SAFE/ +*.zip +*.TIF +*.tif +*.png +*.qgz \ No newline at end of file diff --git a/2024_Grundlagen_Betriebssysteme_Rechnernetze/uebung-3-pthreads b/2024_Grundlagen_Betriebssysteme_Rechnernetze/uebung-3-pthreads new file mode 160000 index 0000000..61606e3 --- /dev/null +++ b/2024_Grundlagen_Betriebssysteme_Rechnernetze/uebung-3-pthreads @@ -0,0 +1 @@ +Subproject commit 61606e311750df81ed51e0ec24b484060cffd3b1 diff --git a/2024_Praxis_der_Programmierung/ue4/char_array.c b/2024_Praxis_der_Programmierung/ue4/char_array.c index 3efb60c..a5c038b 100644 --- a/2024_Praxis_der_Programmierung/ue4/char_array.c +++ b/2024_Praxis_der_Programmierung/ue4/char_array.c @@ -15,11 +15,14 @@ int main() { printf("\nEingabe: %s", eingabe); - for (index = 0; eingabe[index] != '\0'; index++) - if (eingabe[index] == 'a') + for (index = 0; eingabe[index] != '\0'; index++) + if (eingabe[index] == 'a') break; + if (eingabe[index] == '\0') + printf("Der String enthält kein ’a’."); + else + printf("Das erste 'a' ist an Position %d.\n", index + 1); return 0; -} - +} \ No newline at end of file diff --git a/2024_Praxis_der_Programmierung/ue4/fehler.c b/2024_Praxis_der_Programmierung/ue4/fehler.c index 9474cba..23b6994 100644 --- a/2024_Praxis_der_Programmierung/ue4/fehler.c +++ b/2024_Praxis_der_Programmierung/ue4/fehler.c @@ -10,17 +10,18 @@ int main() { int num1, num2; num1 = 0; - nun2 = 1; + num2 = 1; - printf("Der Quotient der Variablen ist: ") + printf("Der Quotient der Variablen ist: "); printf("%d\n", num1/num2); - printf(\n); // Leerzeile + printf("\n"); // Leerzeile printf("Jetzt werden die Variablenwerte vertauscht.\n"); // Dies muss korrigiert werden. (logischer Fehler!!!) + int temp = num1; num1 = num2; - num2 = num1; + num2 = temp; printf("Der Quotient der Vaiablen ist nun: "); printf("%d\n", num1/num2); diff --git a/2024_Praxis_der_Programmierung/ue4/hallo_pdp.c b/2024_Praxis_der_Programmierung/ue4/hallo_pdp.c new file mode 100644 index 0000000..62fb887 --- /dev/null +++ b/2024_Praxis_der_Programmierung/ue4/hallo_pdp.c @@ -0,0 +1,31 @@ +#include +#include +#include + +int main() { + char statS[] = "Hallo, PdP!"; + + char *dynS = malloc( 50 * sizeof(char)); + if (dynS == NULL){ + printf("Speicher konnte nicht reserviert werden"); + return 1; + } + strcpy(dynS, "Hallo, PdP!"); + + printf("%s\n", statS); + printf("%s\n", dynS); + + statS[1] = 'e'; + dynS[1] = 'e'; + printf("%s\n", statS); + printf("%s\n", dynS); + + strcpy(statS, "neuer String"); + strcpy(dynS, "neuer String"); + + printf("%s\n", statS); + printf("%s\n", dynS); + + free(dynS); + return 0; +} diff --git a/2024_Praxis_der_Programmierung/ue4/zeichenketten.c b/2024_Praxis_der_Programmierung/ue4/zeichenketten.c new file mode 100644 index 0000000..3d3be7a --- /dev/null +++ b/2024_Praxis_der_Programmierung/ue4/zeichenketten.c @@ -0,0 +1,62 @@ +#include +#include +#include + +#define MAX 40 + +int main(){ + + char s1 [MAX]; + char s2 [MAX]; + + printf("Bitte Vorname eingeben (max. %d Zeichen): ", MAX-1); + fgets(s1, MAX, stdin); + + size_t len = strlen(s1); + if (len > 0 && s1[len - 1] == '\n') { + s1[len - 1] = '\0'; + } + + printf("Bitte Nachname eingeben (max. %d Zeichen): ", MAX-1); + fgets(s2, MAX, stdin); + + len = strlen(s2); + if (len > 0 && s2[len - 1] == '\n') { + s2[len - 1] = '\0'; + } + + if (strcmp(s1, s2) == 0) { + printf("Vorname und Nachname sind gleich.\n"); + } else { + printf("Vorname und Nachname sind nicht gleich.\n"); + } + + int index = 0; + while(s2[index] != '\0'){ + if (s2[index] >= 'a' && s2[index] <= 'z'){ + s2[index] -= 'a'-'A'; + } + index++; + } + + printf("Einzeln: %s %s \n", s1,s2); + + + char *name = malloc( 80 * sizeof(char)); + if (name==NULL){ + printf("Konnte nicht reserviert werden"); + return 1; + } + + strcpy(name, s1); + strcat(name, " "); + strcat(name, s2); + + + printf("Zusammen: %s \n", name); + printf("Laenge name: %zu\n", strlen(name)); + + free(name); + + return 0; +} \ No newline at end of file diff --git a/2024_Praxis_der_Programmierung/ue5/pkub.c b/2024_Praxis_der_Programmierung/ue5/pkub.c new file mode 100644 index 0000000..f791f4f --- /dev/null +++ b/2024_Praxis_der_Programmierung/ue5/pkub.c @@ -0,0 +1,37 @@ +#include +#include +#include + +int main(int argc, char *argv[]){ + + double value; + double cubic_value; + char *endptr; + + if (argc != 2){ + printf("Nur einen Parameter eingeben\n"); + return 1; + } + + errno = 0; + + value = strtod(argv[1], &endptr); + + if (errno == ERANGE) { + printf("Der Wert ist außerhalb des darstellbaren Bereichs.\n"); + return EXIT_FAILURE; + } + + if (endptr == argv[1]) { + printf("Keine gültige Zahl gefunden.\n"); + return EXIT_FAILURE; + } + + cubic_value = value * value * value; + + printf("Value: %f\n", value); + printf("Cubic Value: %f\n", cubic_value); + + return EXIT_SUCCESS; + +} \ No newline at end of file diff --git a/2024_Praxis_der_Programmierung/ue5/point_2.c b/2024_Praxis_der_Programmierung/ue5/point_2.c new file mode 100644 index 0000000..5f6920e --- /dev/null +++ b/2024_Praxis_der_Programmierung/ue5/point_2.c @@ -0,0 +1,48 @@ +/* point_1.c + * + * Datentyp eines zweidimensionalen Punktes + * verschiedene Methoden der Komponentenselektion + */ + +#include + +struct point { + float x; + float y; +}; + +void move_point(struct point * ptr, float delta_a, float delta_b){ + if (ptr != NULL){ + ptr->x += delta_a; + ptr->y += delta_b; + } +} + +int main() { + struct point p1 = {3.0f,4.0f}; + struct point * ptr = &p1; + struct point p2 = {0.0f,0.0f}; + + printf("\nx-Koordinate von p1: %f", ptr->x); + printf("\ny-Koordinate von p1: %f", ptr->y); + + p2 = p1; + + printf("\nx-Koordinate von p2: %f", p2.x); + printf("\ny-Koordinate von p2: %f", p2.y); + + printf("\nAdresse von p1: %p", (void*)&p1); + printf("\nAdresse von p2: %p", (void*)&p2); + + + printf("\n\n"); + + move_point(ptr, 4.0f , 5.0f); + + printf("\nx-Koordinate von p1: %f", ptr->x); + printf("\ny-Koordinate von p1: %f", ptr->y); + + printf("\n\n"); + + return 0; +} diff --git a/2024_Praxis_der_Programmierung/ue5/point_3.c b/2024_Praxis_der_Programmierung/ue5/point_3.c index c7475b1..6805b58 100644 --- a/2024_Praxis_der_Programmierung/ue5/point_3.c +++ b/2024_Praxis_der_Programmierung/ue5/point_3.c @@ -33,15 +33,16 @@ int main() { printf("\ny-Koordinate von p: %f", ptr->y); free(ptr); + ptr = NULL; - printf("\n"); - printf("\nnach free():"); + // printf("\n"); + // printf("\nnach free():"); - // Dangling Pointer! - printf("\nx-Koordinate von p: %f", ptr->x); - printf("\ny-Koordinate von p: %f", ptr->y); + // // Dangling Pointer! + // printf("\nx-Koordinate von p: %f", ptr->x); + // printf("\ny-Koordinate von p: %f", ptr->y); - printf("\n\n"); + // printf("\n\n"); } return 0; } diff --git a/2024_Praxis_der_Programmierung/ue5/polygonenzug.c b/2024_Praxis_der_Programmierung/ue5/polygonenzug.c new file mode 100644 index 0000000..54cb4ce --- /dev/null +++ b/2024_Praxis_der_Programmierung/ue5/polygonenzug.c @@ -0,0 +1,225 @@ +#include +#include +#include + +// Structs + +typedef struct point { + float x; + float y; +} Point; + +typedef struct Node { + Point *pt; + struct Node *next; +} Node; + +typedef struct Polygonenzug { + Node *head; +} Polygonenzug; + +// Methods + +Point *create_point(float a, float b){ + Point *pt = (Point *)malloc(sizeof(Point)); + if (pt == NULL) { + fprintf(stderr, "Memory allocation failed for Point.\n"); + exit(EXIT_FAILURE); + } + pt->x=a; + pt->y=b; + return pt; +} + +Polygonenzug* create_list() { + Polygonenzug* pzug = (Polygonenzug*)malloc(sizeof(Polygonenzug)); + if (pzug == NULL) { + fprintf(stderr, "Memory allocation failed for Polygonzug.\n"); + exit(EXIT_FAILURE); + } + pzug->head = NULL; + return pzug; +} + +void append(Polygonenzug *pzug, Point *pt){ + Node *new_node = (Node*)malloc(sizeof(Node)); + if (new_node == NULL) { + fprintf(stderr, "Memory allocation failed for new_node in append.\n"); + exit(EXIT_FAILURE); + } + new_node->pt = pt; + new_node->next = NULL; + + if (pzug->head == NULL){ + pzug->head = new_node; + } else { + Node * temp = pzug->head; + while(temp->next != NULL){ + temp = temp->next; + } + temp->next = new_node; + } +} + +void shorten(Polygonenzug *pzug){ + + // Wenn Liste leer + if (pzug->head == NULL){ + printf("The list is empty.\n"); + return; + } + + // Wenn nur ein Knoten + if (pzug->head->next == NULL){ + free(pzug->head->pt); + pzug->head->pt = NULL; + free(pzug->head); + pzug->head = NULL; + printf("Die Liste wurde geloescht.\n"); + return; + } + + + // Wenn mehrere Knoten + Node * temp = pzug->head; + while(temp->next && temp->next->next != NULL){ + temp = temp->next; + } + free(temp->next->pt); + free(temp->next); + temp->next = NULL; +} + +void insert(Polygonenzug *pzug, Point *pt, int index){ + + Node *new_node = (Node*)malloc(sizeof(Node)); + if (new_node == NULL) { + fprintf(stderr, "Memory allocation failed for new_node in insert.\n"); + exit(EXIT_FAILURE); + } + new_node->pt = pt; + new_node->next = NULL; + + // Wenn index 0 + if(index == 0){ + new_node->next = pzug->head; + pzug->head = new_node; + return; + } + + // Groesse der liste checken + int i = 0; + Node *temp = pzug->head; + + if(temp == NULL){ + printf("Cannot insert at index %d in an empty list.\n", index); + free(new_node); + return; + } + + while(temp != NULL && i < index - 1){ + temp = temp->next; + i++; + } + + // Wenn index > 0 + if(temp != NULL){ + new_node->next = temp->next; + temp->next = new_node; + } else{ + printf("Index is outside the list bounds.\n"); + free(new_node); + return; + } +} + +Polygonenzug *mirror(Polygonenzug *pzug){ + Polygonenzug *pzug2 = create_list(); + + Node *temp = pzug->head; + while(temp != NULL){ + Node * new_node = (Node*)malloc(sizeof(Node)); + if (new_node == NULL) { + fprintf(stderr, "Memory allocation failed for new_node in mirror.\n"); + exit(EXIT_FAILURE); + } + Point *new_pt = create_point(temp->pt->x, temp->pt->y); + new_node->pt = new_pt; + new_node->next = pzug2->head; + pzug2->head = new_node; + temp = temp->next; + } + + return pzug2; +} + + +void pretty_print(Polygonenzug *pzug){ + + Node *index = pzug->head; + + while(index != NULL){ + printf("%f , %f\n", index->pt->x, index->pt->y); + index = index->next; + } +} + +void free_list(Polygonenzug *pzug){ + Node *current = pzug->head; + while(current != NULL){ + Node *next_node = current->next; + free(current->pt); + free(current); + current = next_node; + } + free(pzug); +} + + +int main(){ + Point *p1 = create_point(3.0f,4.0f); + Point *p2 = create_point(5.0f,6.0f); + Point *p3 = create_point(8.0f, 9.0f); + + Polygonenzug *pzug = create_list(); + printf("Pzug after Append p1:\n"); + append(pzug, p1); + pretty_print(pzug); + + printf("Pzug after Append p2:\n"); + append(pzug, p2); + pretty_print(pzug); + + printf("Pzug after Insert p3:\n"); + insert(pzug, p3, 1); + pretty_print(pzug); + + Polygonenzug *pzug2 = mirror(pzug); + + // Test shorten + printf("Pzug after Shorten:\n"); + shorten(pzug); + pretty_print(pzug); + printf("Pzug after Shorten:\n"); + shorten(pzug); + pretty_print(pzug); + + // Test shorten liste loeschen + printf("Pzug after Shorten:\n"); + shorten(pzug); + pretty_print(pzug); + + // Test shorten leer + printf("Pzug after Shorten:\n"); + shorten(pzug); + pretty_print(pzug); + + // Test mirror + printf("Pzug2:\n"); + pretty_print(pzug2); + + free_list(pzug); + free_list(pzug2); + + return 0; +} \ No newline at end of file diff --git a/2024_Praxis_der_Programmierung/ue5/typfehler.c b/2024_Praxis_der_Programmierung/ue5/typfehler.c index fc71eac..b5bf585 100644 --- a/2024_Praxis_der_Programmierung/ue5/typfehler.c +++ b/2024_Praxis_der_Programmierung/ue5/typfehler.c @@ -9,19 +9,19 @@ int main() { float n; int rvalue; + int c; - printf("Geben Sie eine Zahl ein: "); - rvalue = scanf("%f", &n); - - // to be debugged: (Tip: Siehe VL4 F18) - printf("Rueckgabewert von scanf: %d\n", rvalue); - - if (rvalue == 0) { - printf("Sie haben keine Zahl eingegeben.\n\n"); - exit(EXIT_FAILURE); - + while (1){ + stdout("Geben Sie eine Zahl ein: "); + rvalue = scanf("%f", &n); + stdout("Rueckgabewert von scanf: %d\n", rvalue); + if (rvalue == 1){ + break; + } else { + while ((c = getchar()) != '\n' && c != EOF) {} + } } - printf("Die Zahl ist %f\n", n); + stdout("Die Zahl ist %f\n", n); return 0; } diff --git a/2024_Praxis_der_Programmierung/ue6/date.c b/2024_Praxis_der_Programmierung/ue6/date.c new file mode 100644 index 0000000..f4989d5 --- /dev/null +++ b/2024_Praxis_der_Programmierung/ue6/date.c @@ -0,0 +1,13 @@ +#include +#include +#include "date.h" + +void setDate(Date* date, int year){ + date->day = 1; + date->month = 1; + date->year = year; +} + +void prettyPrintDate(Date *date){ + printf("Day: %d Month: %d Year: %d \n", date->day,date->month,date->year); +} \ No newline at end of file diff --git a/2024_Praxis_der_Programmierung/ue6/date.h b/2024_Praxis_der_Programmierung/ue6/date.h new file mode 100644 index 0000000..3b60a17 --- /dev/null +++ b/2024_Praxis_der_Programmierung/ue6/date.h @@ -0,0 +1,14 @@ +#ifndef DATE_H +#define DATE_H + +typedef struct date{ + int day; + int month; + int year; +} Date; + +void setDate(Date *date, int year); + +void prettyPrintDate(Date *date); + +#endif \ No newline at end of file diff --git a/2024_Praxis_der_Programmierung/ue6/highscore.c b/2024_Praxis_der_Programmierung/ue6/highscore.c new file mode 100644 index 0000000..41241e6 --- /dev/null +++ b/2024_Praxis_der_Programmierung/ue6/highscore.c @@ -0,0 +1,12 @@ +#include +#include "date.h" +#include "highscore.h" + +void setHighscore(Highscore *highscore, int score, int year){ + setDate(&highscore->date, year); + highscore->score = score; +} + +void prettyPrintHighscore(Highscore *highscore){ + printf("The Highscore on the %d.%d.%d is %d \n", highscore->date.day,highscore->date.month,highscore->date.year,highscore->score); +} diff --git a/2024_Praxis_der_Programmierung/ue6/highscore.h b/2024_Praxis_der_Programmierung/ue6/highscore.h new file mode 100644 index 0000000..9e820d8 --- /dev/null +++ b/2024_Praxis_der_Programmierung/ue6/highscore.h @@ -0,0 +1,14 @@ +#ifndef HIGHSCORE_H +#define HIGHSCORE_H +#include "date.h" + +typedef struct highscore { + Date date; + int score; +} Highscore; + +void setHighscore(Highscore *highscore, int score, int year); + +void prettyPrintHighscore(Highscore *highscore); + +#endif \ No newline at end of file diff --git a/2024_Praxis_der_Programmierung/ue6/rec_count.c b/2024_Praxis_der_Programmierung/ue6/rec_count.c index 21432e8..adca7e9 100644 --- a/2024_Praxis_der_Programmierung/ue6/rec_count.c +++ b/2024_Praxis_der_Programmierung/ue6/rec_count.c @@ -2,6 +2,11 @@ #include +int count = 1; + + +void rec_out(int n); + void decr(int n) { rec_out(--n); @@ -9,7 +14,6 @@ void decr(int n) void rec_out(int n) { - int count = 1; printf("Die %d. Ausgabe.\n", count); count++; if (n > 1) diff --git a/2024_Praxis_der_Programmierung/ue6/scores.c b/2024_Praxis_der_Programmierung/ue6/scores.c new file mode 100644 index 0000000..0980956 --- /dev/null +++ b/2024_Praxis_der_Programmierung/ue6/scores.c @@ -0,0 +1,15 @@ +#include +#include +#include "date.h" +#include "highscore.h" + +int main(){ + + int year1 = 2024; + int score1 = 10; + + Highscore highscore; + setHighscore(&highscore, score1, year1); + prettyPrintHighscore(&highscore); + +} \ No newline at end of file diff --git a/2024_Remote_Sensing/lab2/Lab2.md b/2024_Remote_Sensing/lab2/Lab2.md new file mode 100644 index 0000000..8dc6812 --- /dev/null +++ b/2024_Remote_Sensing/lab2/Lab2.md @@ -0,0 +1,31 @@ +# Lab Assignment 2 + +Name: Joaquin Gottlebe +Matrikelnummer: 829101 + +## Question 1 + +![](maps/ndvi_ndwi.png) + +## Question 2 + +![](maps/surface_temp_b10.png) +![](maps/surface_temp_b11.png) + +The original Raster file included values from -125 C which are not realistic in germany. Also Maximum temperatures from 44 degrees are not common, but could be explainable through realy heat conductiv material on very small areas. + +Between Water bodies and land bodies. and also the north and south of the area. + +## Question 3 + +![](maps/true_color.png) +![](maps/false_color.png) + +The main differences of the map are the Bands and their display of it. For example the true color composite consists of band 2,3 and 4 in the order of RGB as 4,3,2. The false color composite consists of band 3,4 and 5 in the order ich which Vegetation is the most visible (5,4,3). + +## Question 4 + +![](maps/ndvi_chiloe_2013.png) +![](maps/ndvi_chiloe_2024.png) + +I choose this study area because i was on this island Chiloe in Chile for a bit and noticed a lot of National parks and vegeation there and i wanted to know how this changed in these past years. These two images show difference in NDVI between 2013 and 2024 of the island Chiloe in Chile and surroundings. I chose NDVI as it shows the difference in vegetation very well. \ No newline at end of file diff --git a/2024_Remote_Sensing/lab2/Makefile b/2024_Remote_Sensing/lab2/Makefile new file mode 100644 index 0000000..071fdbd --- /dev/null +++ b/2024_Remote_Sensing/lab2/Makefile @@ -0,0 +1,23 @@ +.PHONY: all run_scripts + +# Path to the Conda environment and scripts +CONDA_ENV_NAME = gdal_env +ACTIVATE_SCRIPT = ~/miniforge3/bin/activate +SCRIPTS = top_of_atmosphere.py ndvi.py ndwi.py +PANDOC=pandoc + +all: run_scripts + +run_scripts: + @echo "Activating Conda environment and running scripts..." + @source $(ACTIVATE_SCRIPT) && conda activate $(CONDA_ENV_NAME) && \ + python3 toa.py data/LC08_L1TP_193023_20170602_20170615_01_T1_MTL.txt toa data/LC08_L1TP_193023_20170602_20170615_01_T1_B2.TIF data/LC08_L1TP_193023_20170602_20170615_01_T1_B4.TIF data/LC08_L1TP_193023_20170602_20170615_01_T1_B3.TIF data/LC08_L1TP_193023_20170602_20170615_01_T1_B5.TIF && \ + python3 toa_radiance.py data/LC08_L1TP_193023_20170602_20170615_01_T1_MTL.txt toa_radiance data/LC08_L1TP_193023_20170602_20170615_01_T1_B10.TIF data/LC08_L1TP_193023_20170602_20170615_01_T1_B11.TIF && \ + python3 ndvi.py toa/LC08_L1TP_193023_20170602_20170615_01_T1_B4_toa.TIF toa/LC08_L1TP_193023_20170602_20170615_01_T1_B5_toa.TIF ndvi && \ + python3 ndwi.py toa/LC08_L1TP_193023_20170602_20170615_01_T1_B3_toa.TIF toa/LC08_L1TP_193023_20170602_20170615_01_T1_B5_toa.TIF ndwi && \ + python3 surface_temperature.py data/LC08_L1TP_193023_20170602_20170615_01_T1_MTL.txt toa_radiance/LC08_L1TP_193023_20170602_20170615_01_T1_B10_toa_radiance.TIF toa_radiance/LC08_L1TP_193023_20170602_20170615_01_T1_B11_toa_radiance.TIF surface_temperature && \ + python3 toa.py data2/LC08_L1TP_233089_20131224_20200912_02_T1_MTL.txt toa2 data2/LC08_L1TP_233089_20131224_20200912_02_T1_B2.TIF data2/LC08_L1TP_233089_20131224_20200912_02_T1_B4.TIF data2/LC08_L1TP_233089_20131224_20200912_02_T1_B3.TIF data2/LC08_L1TP_233089_20131224_20200912_02_T1_B5.TIF && \ + python3 toa.py data3/LC08_L1TP_233089_20240410_20240419_02_T1_MTL.txt toa3 data3/LC08_L1TP_233089_20240410_20240419_02_T1_B2.TIF data3/LC08_L1TP_233089_20240410_20240419_02_T1_B4.TIF data3/LC08_L1TP_233089_20240410_20240419_02_T1_B3.TIF data3/LC08_L1TP_233089_20240410_20240419_02_T1_B5.TIF && \ + python3 ndvi.py toa2/LC08_L1TP_233089_20131224_20200912_02_T1_B4_toa.TIF toa2/LC08_L1TP_233089_20131224_20200912_02_T1_B5_toa.TIF ndvi2 && \ + python3 ndvi.py toa3/LC08_L1TP_233089_20240410_20240419_02_T1_B4_toa.TIF toa3/LC08_L1TP_233089_20240410_20240419_02_T1_B5_toa.TIF ndvi3 && \ + $(PANDOC) Lab2.md -o rcm01_2425_lab2_gottlebe_829101.pdf \ No newline at end of file diff --git a/2024_Remote_Sensing/lab2/crop_tiff.py b/2024_Remote_Sensing/lab2/crop_tiff.py new file mode 100644 index 0000000..e69de29 diff --git a/2024_Remote_Sensing/lab2/false_color_map.py b/2024_Remote_Sensing/lab2/false_color_map.py new file mode 100644 index 0000000..0f02881 --- /dev/null +++ b/2024_Remote_Sensing/lab2/false_color_map.py @@ -0,0 +1,119 @@ +# Cleanup +project = QgsProject.instance() +layer_ids = list(project.mapLayers().keys()) +for layer_id in layer_ids: + project.removeMapLayer(layer_id) + +# Loading B10 Layer +path_st_b10 = "/home/huaqo/dev/courses/2024_Remote_Sensing/lab2/false_color/false_color.tif" +st_b10_layer = QgsRasterLayer(path_st_b10, "False Color Composite") +if not st_b10_layer.isValid(): + print("st_b10 layer failed to load!") +else: + QgsProject.instance().addMapLayer(st_b10_layer) + print("st_b10 Layer loaded!") + +# Styling Layer +st_b10_style_path = '/home/huaqo/dev/courses/2024_Remote_Sensing/styles/false_color.qml' +if not st_b10_layer.loadNamedStyle(st_b10_style_path): + print("Failed to load st_b10 style!") +else: + print("st_b10 style loaded!") + +# Load layout +layout = QgsPrintLayout(project) +layout.initializeDefaults() + +# Map +map_item = QgsLayoutItemMap(layout) +map_item.attemptMove(QgsLayoutPoint(20, 30, QgsUnitTypes.LayoutMillimeters)) # Added space between top and map +map_item.attemptResize(QgsLayoutSize(200, 150, QgsUnitTypes.LayoutMillimeters)) # Reduced height for better spacing +map_item.zoomToExtent(st_b10_layer.extent()) +layout.addLayoutItem(map_item) + +# Add Coordinate Grid +grid = map_item.grid() +grid.setEnabled(True) +grid.setIntervalX(100000) # Set interval for grid lines in map units +grid.setIntervalY(100000) +grid.setAnnotationEnabled(True) +grid.setAnnotationPrecision(0) # Set annotation precision for grid labels +grid.setFrameStyle(QgsLayoutItemMapGrid.Zebra) # Optional grid style + +# Label +title_label = QgsLayoutItemLabel(layout) +title_label.setText("False Color Composite - Lab 2") +title_label.setFont(QFont("Arial", 16)) +title_label.setHAlign(Qt.AlignCenter) +title_label.attemptMove(QgsLayoutPoint(105, 10, QgsUnitTypes.LayoutMillimeters)) +title_label.adjustSizeToText() +layout.addLayoutItem(title_label) + +# Legend +legend_item = QgsLayoutItemLegend(layout) +legend_item.setLinkedMap(map_item) +legend_item.setTitle("Legend") +legend_item.setFrameEnabled(True) +legend_item.attemptMove(QgsLayoutPoint(245, 30, QgsUnitTypes.LayoutMillimeters)) # Adjusted for better alignment and spacing +layout.addLayoutItem(legend_item) + +# Scalebar +scalebar_item = QgsLayoutItemScaleBar(layout) +scalebar_item.setStyle('Single Box') # Change to 'Single Box' or 'Double Box' for a proper scalebar +scalebar_item.setLinkedMap(map_item) +scalebar_item.setUnitLabel('m') # Add unit label to make it informative +scalebar_item.setNumberOfSegments(4) # Specify the number of segments in the scalebar +scalebar_item.setNumberOfSegmentsLeft(0) # Segments to the left of zero (if any) +scalebar_item.setUnitsPerSegment(50000) # Set distance per segment in map units +scalebar_item.setFont(QFont('Arial', 10)) # Set font size for better readability +scalebar_item.setHeight(5) # Set the height of the scalebar +scalebar_item.applyDefaultSize() +scalebar_item.attemptMove(QgsLayoutPoint(20, 190, QgsUnitTypes.LayoutMillimeters)) # Moved below the map for better organization +layout.addLayoutItem(scalebar_item) + +# North Arrow +north_arrow_item = QgsLayoutItemPicture(layout) +north_arrow_item.setPicturePath('/home/huaqo/dev/courses/2024_Remote_Sensing/styles/north_arrow.svg') +if not north_arrow_item.picturePath(): + print("North arrow image failed to load!") +north_arrow_item.setReferencePoint(QgsLayoutItemPicture.UpperLeft) +north_arrow_item.attemptMove(QgsLayoutPoint(30, 40, QgsUnitTypes.LayoutMillimeters)) # Adjusted placement for better visibility +north_arrow_item.attemptResize(QgsLayoutSize(15, 15, QgsUnitTypes.LayoutMillimeters)) # Increased size for better visibility +layout.addLayoutItem(north_arrow_item) + +# Author +info_label = QgsLayoutItemLabel(layout) +info_label.setText("Author: Joaquin Gottlebe") +info_label.setFont(QFont("Arial", 12)) +info_label.setHAlign(Qt.AlignLeft) +info_label.attemptMove(QgsLayoutPoint(245, 100, QgsUnitTypes.LayoutMillimeters)) +info_label.adjustSizeToText() +layout.addLayoutItem(info_label) + +# Author +info_label = QgsLayoutItemLabel(layout) +info_label.setText("Date: 16.11.2025") +info_label.setFont(QFont("Arial", 12)) +info_label.setHAlign(Qt.AlignLeft) +info_label.attemptMove(QgsLayoutPoint(245, 110, QgsUnitTypes.LayoutMillimeters)) +info_label.adjustSizeToText() +layout.addLayoutItem(info_label) + +# Coordinate System +info_label = QgsLayoutItemLabel(layout) +info_label.setText("CRS: EPSG 32633") +info_label.setFont(QFont("Arial", 12)) +info_label.setHAlign(Qt.AlignLeft) +info_label.attemptMove(QgsLayoutPoint(245, 120, QgsUnitTypes.LayoutMillimeters)) +info_label.adjustSizeToText() +layout.addLayoutItem(info_label) + + +# Export +export_path = "/home/huaqo/dev/courses/2024_Remote_Sensing/lab2/maps/false_color.png" +exporter = QgsLayoutExporter(layout) +export_result = exporter.exportToImage(export_path, QgsLayoutExporter.ImageExportSettings()) +if export_result != QgsLayoutExporter.Success: + print("Failed to export map!") +else: + print("Map exported to:", export_path) diff --git a/2024_Remote_Sensing/lab2/missfont.log b/2024_Remote_Sensing/lab2/missfont.log new file mode 100644 index 0000000..2566c32 --- /dev/null +++ b/2024_Remote_Sensing/lab2/missfont.log @@ -0,0 +1,2 @@ +mktextfm ecrm1000 +mktextfm ecrm1000 diff --git a/2024_Remote_Sensing/lab2/ndvi.py b/2024_Remote_Sensing/lab2/ndvi.py new file mode 100644 index 0000000..0d58fa4 --- /dev/null +++ b/2024_Remote_Sensing/lab2/ndvi.py @@ -0,0 +1,73 @@ +from osgeo import gdal +import numpy as np +import sys +import os +import re + +def tif_to_array(path): + data = gdal.Open(path) + if data is None: + raise FileNotFoundError(f"Cannot open file: {path}") + band = data.GetRasterBand(1) + array = band.ReadAsArray().astype(np.float32) + return data, array + +def handle_nodata(array): + clean_array = np.where((array < 0) | (array > 1), np.nan, array) + return clean_array + +def ndvi_calc(NIR,RED): + ndvi_array = (NIR - RED) / (NIR + RED + 1e-10) + return ndvi_array + +def save_tif(array,data,output_path): + driver = gdal.GetDriverByName('GTiff') + rows, cols = array.shape + out_data = driver.Create(output_path, cols, rows, 1, gdal.GDT_Float32) + out_data.SetGeoTransform(data.GetGeoTransform()) + out_data.SetProjection(data.GetProjection()) + out_band = out_data.GetRasterBand(1) + out_band.WriteArray(array) + out_band.SetNoDataValue(-9999) + out_band.FlushCache() + out_data = None + +def get_output_path(output_dir,input_path): + basename = os.path.basename(input_path) + basename_no_ext = os.path.splitext(basename)[0] + basename_no_ext = re.sub(r'_B[1-9]|_B10|_B11', '', basename_no_ext) + path = os.path.join(output_dir, f"{basename_no_ext}_ndvi.TIF") + return path + +def print_min_max(array, label): + print(label," Min:", np.nanmin(array), " Max:",np.nanmax(array)) + +def main(): + if len(sys.argv) != 4: + print("Usage: python3 script.py ") + sys.exit(1) + + red_path = sys.argv[1] + nir_path = sys.argv[2] + output_dir = sys.argv[3] + + red_data, red_array = tif_to_array(red_path) + nir_data, nir_array = tif_to_array(nir_path) + red_array = handle_nodata(red_array) + nir_array = handle_nodata(nir_array) + + ndvi_array = ndvi_calc(nir_array, red_array) + + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + output_path = get_output_path(output_dir,red_path) + + save_tif(ndvi_array, red_data, output_path) + + print_min_max(red_array, "RED") + print_min_max(nir_array, "NIR") + print_min_max(ndvi_array, "NDVI") + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/2024_Remote_Sensing/lab2/ndvi_map.py b/2024_Remote_Sensing/lab2/ndvi_map.py new file mode 100644 index 0000000..c7523b8 --- /dev/null +++ b/2024_Remote_Sensing/lab2/ndvi_map.py @@ -0,0 +1,119 @@ +# Cleanup +project = QgsProject.instance() +layer_ids = list(project.mapLayers().keys()) +for layer_id in layer_ids: + project.removeMapLayer(layer_id) + +# Loading NDVI Layer +path_ndvi = "/home/huaqo/dev/courses/2024_Remote_Sensing/lab2/ndvi3/LC08_L1TP_233089_20240410_20240419_02_T1_toa_ndvi.TIF" +ndvi_layer = QgsRasterLayer(path_ndvi, "NDVI") +if not ndvi_layer.isValid(): + print("NDVI layer failed to load!") +else: + QgsProject.instance().addMapLayer(ndvi_layer) + print("NDVI Layer loaded!") + +# Styling NDVI Layer +ndvi_style_path = '/home/huaqo/dev/courses/2024_Remote_Sensing/styles/ndvi.qml' +if not ndvi_layer.loadNamedStyle(ndvi_style_path): + print("Failed to load NDVI style!") +else: + print("NDVI style loaded!") + +# Load layout +layout = QgsPrintLayout(project) +layout.initializeDefaults() + +# Map +map_item = QgsLayoutItemMap(layout) +map_item.attemptMove(QgsLayoutPoint(20, 30, QgsUnitTypes.LayoutMillimeters)) # Added space between top and map +map_item.attemptResize(QgsLayoutSize(200, 150, QgsUnitTypes.LayoutMillimeters)) # Reduced height for better spacing +map_item.zoomToExtent(ndvi_layer.extent()) +layout.addLayoutItem(map_item) + +# Add Coordinate Grid +grid = map_item.grid() +grid.setEnabled(True) +grid.setIntervalX(100000) # Set interval for grid lines in map units +grid.setIntervalY(100000) +grid.setAnnotationEnabled(True) +grid.setAnnotationPrecision(0) # Set annotation precision for grid labels +grid.setFrameStyle(QgsLayoutItemMapGrid.Zebra) # Optional grid style + +# Label +title_label = QgsLayoutItemLabel(layout) +title_label.setText("NDVI Chiloe 2024 - Lab 2") +title_label.setFont(QFont("Arial", 16)) +title_label.setHAlign(Qt.AlignCenter) +title_label.attemptMove(QgsLayoutPoint(105, 10, QgsUnitTypes.LayoutMillimeters)) +title_label.adjustSizeToText() +layout.addLayoutItem(title_label) + +# Legend +legend_item = QgsLayoutItemLegend(layout) +legend_item.setLinkedMap(map_item) +legend_item.setTitle("Legend") +legend_item.setFrameEnabled(True) +legend_item.attemptMove(QgsLayoutPoint(245, 30, QgsUnitTypes.LayoutMillimeters)) # Adjusted for better alignment and spacing +layout.addLayoutItem(legend_item) + +# Scalebar +scalebar_item = QgsLayoutItemScaleBar(layout) +scalebar_item.setStyle('Single Box') # Change to 'Single Box' or 'Double Box' for a proper scalebar +scalebar_item.setLinkedMap(map_item) +scalebar_item.setUnitLabel('m') # Add unit label to make it informative +scalebar_item.setNumberOfSegments(4) # Specify the number of segments in the scalebar +scalebar_item.setNumberOfSegmentsLeft(0) # Segments to the left of zero (if any) +scalebar_item.setUnitsPerSegment(50000) # Set distance per segment in map units +scalebar_item.setFont(QFont('Arial', 10)) # Set font size for better readability +scalebar_item.setHeight(5) # Set the height of the scalebar +scalebar_item.applyDefaultSize() +scalebar_item.attemptMove(QgsLayoutPoint(20, 190, QgsUnitTypes.LayoutMillimeters)) # Moved below the map for better organization +layout.addLayoutItem(scalebar_item) + +# North Arrow +north_arrow_item = QgsLayoutItemPicture(layout) +north_arrow_item.setPicturePath('/home/huaqo/dev/courses/2024_Remote_Sensing/styles/north_arrow.svg') +if not north_arrow_item.picturePath(): + print("North arrow image failed to load!") +north_arrow_item.setReferencePoint(QgsLayoutItemPicture.UpperLeft) +north_arrow_item.attemptMove(QgsLayoutPoint(30, 40, QgsUnitTypes.LayoutMillimeters)) # Adjusted placement for better visibility +north_arrow_item.attemptResize(QgsLayoutSize(15, 15, QgsUnitTypes.LayoutMillimeters)) # Increased size for better visibility +layout.addLayoutItem(north_arrow_item) + +# Author +info_label = QgsLayoutItemLabel(layout) +info_label.setText("Author: Joaquin Gottlebe") +info_label.setFont(QFont("Arial", 12)) +info_label.setHAlign(Qt.AlignLeft) +info_label.attemptMove(QgsLayoutPoint(245, 100, QgsUnitTypes.LayoutMillimeters)) +info_label.adjustSizeToText() +layout.addLayoutItem(info_label) + +# Author +info_label = QgsLayoutItemLabel(layout) +info_label.setText("Date: 16.11.2025") +info_label.setFont(QFont("Arial", 12)) +info_label.setHAlign(Qt.AlignLeft) +info_label.attemptMove(QgsLayoutPoint(245, 110, QgsUnitTypes.LayoutMillimeters)) +info_label.adjustSizeToText() +layout.addLayoutItem(info_label) + +# Coordinate System +info_label = QgsLayoutItemLabel(layout) +info_label.setText("CRS: EPSG 32633") +info_label.setFont(QFont("Arial", 12)) +info_label.setHAlign(Qt.AlignLeft) +info_label.attemptMove(QgsLayoutPoint(245, 120, QgsUnitTypes.LayoutMillimeters)) +info_label.adjustSizeToText() +layout.addLayoutItem(info_label) + + +# Export +export_path = "/home/huaqo/dev/courses/2024_Remote_Sensing/lab2/maps/ndvi_chiloe_2024.png" +exporter = QgsLayoutExporter(layout) +export_result = exporter.exportToImage(export_path, QgsLayoutExporter.ImageExportSettings()) +if export_result != QgsLayoutExporter.Success: + print("Failed to export map!") +else: + print("Map exported to:", export_path) diff --git a/2024_Remote_Sensing/lab2/ndvi_ndwi_map.py b/2024_Remote_Sensing/lab2/ndvi_ndwi_map.py new file mode 100644 index 0000000..ed95ac0 --- /dev/null +++ b/2024_Remote_Sensing/lab2/ndvi_ndwi_map.py @@ -0,0 +1,139 @@ +# Cleanup +project = QgsProject.instance() +layer_ids = list(project.mapLayers().keys()) +for layer_id in layer_ids: + project.removeMapLayer(layer_id) + +# Loading NDVI Layer +path_ndvi = "/home/huaqo/dev/courses/2024_Remote_Sensing/lab2/ndvi/LC08_L1TP_193023_20170602_20170615_01_T1_toa_ndvi.TIF" +ndvi_layer = QgsRasterLayer(path_ndvi, "NDVI") +if not ndvi_layer.isValid(): + print("NDVI layer failed to load!") +else: + QgsProject.instance().addMapLayer(ndvi_layer) + print("NDVI Layer loaded!") + +# Loading Additional Layer +path_overlay = "/home/huaqo/dev/courses/2024_Remote_Sensing/lab2/ndwi/LC08_L1TP_193023_20170602_20170615_01_T1_toa_ndwi.TIF" # Replace with your actual layer path +ndwi_layer = QgsRasterLayer(path_overlay, "NDWI") +if not ndwi_layer.isValid(): + print("NDWI layer failed to load!") +else: + QgsProject.instance().addMapLayer(ndwi_layer) + print("NDWI Layer loaded!") + +# Styling NDVI Layer +ndvi_style_path = '/home/huaqo/dev/courses/2024_Remote_Sensing/styles/ndvi.qml' +if not ndvi_layer.loadNamedStyle(ndvi_style_path): + print("Failed to load NDVI style!") +else: + print("NDVI style loaded!") + +ndwi_style_path = '/home/huaqo/dev/courses/2024_Remote_Sensing/styles/ndwi.qml' +if not ndwi_layer.loadNamedStyle(ndwi_style_path): + print("Failed to load NDWI style!") +else: + print("NDWI style loaded!") + +# Reference system (optional, uncomment if needed) +# crs = QgsCoordinateReferenceSystem("EPSG:4326") +# crsSrc = QgsCoordinateReferenceSystem("EPSG:32633") +# crsDest = QgsCoordinateReferenceSystem("EPSG:4326") + +# Load layout +layout = QgsPrintLayout(project) +layout.initializeDefaults() + +# Map +map_item = QgsLayoutItemMap(layout) +map_item.attemptMove(QgsLayoutPoint(20, 30, QgsUnitTypes.LayoutMillimeters)) # Added space between top and map +map_item.attemptResize(QgsLayoutSize(200, 150, QgsUnitTypes.LayoutMillimeters)) # Reduced height for better spacing +map_item.zoomToExtent(ndvi_layer.extent()) +layout.addLayoutItem(map_item) + +# Add Coordinate Grid +grid = map_item.grid() +grid.setEnabled(True) +grid.setIntervalX(100000) # Set interval for grid lines in map units +grid.setIntervalY(100000) +grid.setAnnotationEnabled(True) +grid.setAnnotationPrecision(0) # Set annotation precision for grid labels +grid.setFrameStyle(QgsLayoutItemMapGrid.Zebra) # Optional grid style + +# Label +title_label = QgsLayoutItemLabel(layout) +title_label.setText("NDVI & NDWI - Lab 2") +title_label.setFont(QFont("Arial", 16)) +title_label.setHAlign(Qt.AlignCenter) +title_label.attemptMove(QgsLayoutPoint(105, 10, QgsUnitTypes.LayoutMillimeters)) +title_label.adjustSizeToText() +layout.addLayoutItem(title_label) + +# Legend +legend_item = QgsLayoutItemLegend(layout) +legend_item.setLinkedMap(map_item) +legend_item.setTitle("Legend") +legend_item.setFrameEnabled(True) +legend_item.attemptMove(QgsLayoutPoint(245, 30, QgsUnitTypes.LayoutMillimeters)) # Adjusted for better alignment and spacing +layout.addLayoutItem(legend_item) + +# Scalebar +scalebar_item = QgsLayoutItemScaleBar(layout) +scalebar_item.setStyle('Single Box') # Change to 'Single Box' or 'Double Box' for a proper scalebar +scalebar_item.setLinkedMap(map_item) +scalebar_item.setUnitLabel('m') # Add unit label to make it informative +scalebar_item.setNumberOfSegments(4) # Specify the number of segments in the scalebar +scalebar_item.setNumberOfSegmentsLeft(0) # Segments to the left of zero (if any) +scalebar_item.setUnitsPerSegment(50000) # Set distance per segment in map units +scalebar_item.setFont(QFont('Arial', 10)) # Set font size for better readability +scalebar_item.setHeight(5) # Set the height of the scalebar +scalebar_item.applyDefaultSize() +scalebar_item.attemptMove(QgsLayoutPoint(20, 190, QgsUnitTypes.LayoutMillimeters)) # Moved below the map for better organization +layout.addLayoutItem(scalebar_item) + +# North Arrow +north_arrow_item = QgsLayoutItemPicture(layout) +north_arrow_item.setPicturePath('/home/huaqo/dev/courses/2024_Remote_Sensing/styles/north_arrow.svg') +if not north_arrow_item.picturePath(): + print("North arrow image failed to load!") +north_arrow_item.setReferencePoint(QgsLayoutItemPicture.UpperLeft) +north_arrow_item.attemptMove(QgsLayoutPoint(30, 40, QgsUnitTypes.LayoutMillimeters)) # Adjusted placement for better visibility +north_arrow_item.attemptResize(QgsLayoutSize(15, 15, QgsUnitTypes.LayoutMillimeters)) # Increased size for better visibility +layout.addLayoutItem(north_arrow_item) + +# Author +info_label = QgsLayoutItemLabel(layout) +info_label.setText("Author: Joaquin Gottlebe") +info_label.setFont(QFont("Arial", 12)) +info_label.setHAlign(Qt.AlignLeft) +info_label.attemptMove(QgsLayoutPoint(245, 100, QgsUnitTypes.LayoutMillimeters)) +info_label.adjustSizeToText() +layout.addLayoutItem(info_label) + +# Author +info_label = QgsLayoutItemLabel(layout) +info_label.setText("Date: 16.11.2025") +info_label.setFont(QFont("Arial", 12)) +info_label.setHAlign(Qt.AlignLeft) +info_label.attemptMove(QgsLayoutPoint(245, 110, QgsUnitTypes.LayoutMillimeters)) +info_label.adjustSizeToText() +layout.addLayoutItem(info_label) + +# Coordinate System +info_label = QgsLayoutItemLabel(layout) +info_label.setText("CRS: EPSG 32633") +info_label.setFont(QFont("Arial", 12)) +info_label.setHAlign(Qt.AlignLeft) +info_label.attemptMove(QgsLayoutPoint(245, 120, QgsUnitTypes.LayoutMillimeters)) +info_label.adjustSizeToText() +layout.addLayoutItem(info_label) + + +# Export +export_path = "/home/huaqo/dev/courses/2024_Remote_Sensing/lab2/maps/ndvi_ndwi.png" +exporter = QgsLayoutExporter(layout) +export_result = exporter.exportToImage(export_path, QgsLayoutExporter.ImageExportSettings()) +if export_result != QgsLayoutExporter.Success: + print("Failed to export map!") +else: + print("Map exported to:", export_path) diff --git a/2024_Remote_Sensing/lab2/ndwi.py b/2024_Remote_Sensing/lab2/ndwi.py new file mode 100644 index 0000000..5b07faf --- /dev/null +++ b/2024_Remote_Sensing/lab2/ndwi.py @@ -0,0 +1,81 @@ +from osgeo import gdal +import numpy as np +import sys +import os +import re + + +def tif_to_array(path): + data = gdal.Open(path) + if data is None: + raise FileNotFoundError(f"Cannot open file: {path}") + band = data.GetRasterBand(1) + array = band.ReadAsArray().astype(np.float32) + return data, array + + +def handle_nodata(array): + clean_array = np.where((array < 0) | (array > 1), np.nan, array) + return clean_array + + +def ndwi_calc(GREEN, NIR): + ndwi_array = (GREEN - NIR) / (GREEN + NIR + 1e-10) + return ndwi_array + + +def save_tif(array, data, output_path): + driver = gdal.GetDriverByName('GTiff') + rows, cols = array.shape + out_data = driver.Create(output_path, cols, rows, 1, gdal.GDT_Float32) + out_data.SetGeoTransform(data.GetGeoTransform()) + out_data.SetProjection(data.GetProjection()) + out_band = out_data.GetRasterBand(1) + out_band.WriteArray(array) + out_band.SetNoDataValue(-9999) + out_band.FlushCache() + out_data = None + + +def get_output_path(output_dir, input_path): + basename = os.path.basename(input_path) + basename_no_ext = os.path.splitext(basename)[0] + basename_no_ext = re.sub(r'_B[1-9]|_B10|_B11', '', basename_no_ext) + path = os.path.join(output_dir, f"{basename_no_ext}_ndwi.TIF") + return path + + +def print_min_max(array, label): + print(label, " Min:", np.nanmin(array), " Max:", np.nanmax(array)) + + +def main(): + if len(sys.argv) != 4: + print("Usage: python3 script.py ") + sys.exit(1) + + green_path = sys.argv[1] + nir_path = sys.argv[2] + output_dir = sys.argv[3] + + green_data, green_array = tif_to_array(green_path) + nir_data, nir_array = tif_to_array(nir_path) + green_array = handle_nodata(green_array) + nir_array = handle_nodata(nir_array) + + ndwi_array = ndwi_calc(green_array, nir_array) + + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + output_path = get_output_path(output_dir, green_path) + + save_tif(ndwi_array, green_data, output_path) + + print_min_max(nir_array, "NIR") + print_min_max(green_array, "GREEN") + print_min_max(ndwi_array, "NDWI") + + +if __name__ == "__main__": + main() diff --git a/2024_Remote_Sensing/lab2/surface_temperature.py b/2024_Remote_Sensing/lab2/surface_temperature.py new file mode 100644 index 0000000..8047475 --- /dev/null +++ b/2024_Remote_Sensing/lab2/surface_temperature.py @@ -0,0 +1,81 @@ +import sys +import os +from osgeo import gdal +import numpy as np +import math + +def load_metadata(metadata_path): + K_CONSTANTS = {} + with open(metadata_path, "r") as file: + metadata_lines = file.readlines() + for line in metadata_lines: + if "K1_CONSTANT_BAND_10" in line or "K2_CONSTANT_BAND_10" in line or \ + "K1_CONSTANT_BAND_11" in line or "K2_CONSTANT_BAND_11" in line: + variable, value = line.split(" = ") + K_CONSTANTS[variable.strip()] = float(value.strip()) + return K_CONSTANTS + +def load_band(band_path): + band_data = gdal.Open(band_path) + band = band_data.GetRasterBand(1) + band_array = band.ReadAsArray().astype(np.float32) + no_data_value = band.GetNoDataValue() + band_array = np.where(band_array == no_data_value, np.nan, band_array) + return band_data, band_array + +def calculate_surface_temperature(band_array, K1_CONSTANT, K2_CONSTANT): + epsilon = 1e-10 + result_array = (K2_CONSTANT / np.log((K1_CONSTANT / (band_array + epsilon)) + 1)) - 273.15 + return result_array + +def save_result(output_directory, band_path, result_array, band_data): + if not os.path.exists(output_directory): + os.makedirs(output_directory) + driver = gdal.GetDriverByName('GTiff') + rows, cols = result_array.shape + band_filename = os.path.basename(band_path) + output_filename = os.path.splitext(band_filename)[0] + '_surfacetemp.TIF' + output_path = os.path.join(output_directory, output_filename) + out_data = driver.Create(output_path, cols, rows, 1, gdal.GDT_Float32) + out_data.SetGeoTransform(band_data.GetGeoTransform()) + out_data.SetProjection(band_data.GetProjection()) + out_band = out_data.GetRasterBand(1) + out_band.WriteArray(result_array) + out_band.SetNoDataValue(0) + out_band.FlushCache() + out_data = None + +def main(): + metadata_path = sys.argv[1] + band_paths = { + 'BAND_10': sys.argv[2], + 'BAND_11': sys.argv[3] + } + output_directory = sys.argv[4] + + # Load metadata + K_CONSTANTS = load_metadata(metadata_path) + print(K_CONSTANTS) + + # Process each band + for band_key, band_path in band_paths.items(): + # Load band + band_data, band_array = load_band(band_path) + + # Get K1 and K2 constants + K1_CONSTANT = K_CONSTANTS[f'K1_CONSTANT_{band_key}'] + K2_CONSTANT = K_CONSTANTS[f'K2_CONSTANT_{band_key}'] + + # Calculate surface temperature + result_array = calculate_surface_temperature(band_array, K1_CONSTANT, K2_CONSTANT) + + # Save the result + save_result(output_directory, band_path, result_array, band_data) + + # Print min and max of the arrays + print(f"{band_key} Min, Max:", np.nanmin(band_array), np.nanmax(band_array)) + print(f"{band_key} Surface Temp Min, Max:", np.nanmin(result_array), np.nanmax(result_array)) + +if __name__ == "__main__": + main() + diff --git a/2024_Remote_Sensing/lab2/surface_temperature_map.py b/2024_Remote_Sensing/lab2/surface_temperature_map.py new file mode 100644 index 0000000..341426f --- /dev/null +++ b/2024_Remote_Sensing/lab2/surface_temperature_map.py @@ -0,0 +1,119 @@ +# Cleanup +project = QgsProject.instance() +layer_ids = list(project.mapLayers().keys()) +for layer_id in layer_ids: + project.removeMapLayer(layer_id) + +# Loading B10 Layer +path_st_b10 = "/home/huaqo/dev/courses/2024_Remote_Sensing/lab2/surface_temperature/LC08_L1TP_193023_20170602_20170615_01_T1_B11_toa_radiance_surfacetemp.TIF" +st_b10_layer = QgsRasterLayer(path_st_b10, "Surf. Temp. B11") +if not st_b10_layer.isValid(): + print("st_b10 layer failed to load!") +else: + QgsProject.instance().addMapLayer(st_b10_layer) + print("st_b10 Layer loaded!") + +# Styling Layer +st_b10_style_path = '/home/huaqo/dev/courses/2024_Remote_Sensing/styles/surface_temperature.qml' +if not st_b10_layer.loadNamedStyle(st_b10_style_path): + print("Failed to load st_b10 style!") +else: + print("st_b10 style loaded!") + +# Load layout +layout = QgsPrintLayout(project) +layout.initializeDefaults() + +# Map +map_item = QgsLayoutItemMap(layout) +map_item.attemptMove(QgsLayoutPoint(20, 30, QgsUnitTypes.LayoutMillimeters)) # Added space between top and map +map_item.attemptResize(QgsLayoutSize(200, 150, QgsUnitTypes.LayoutMillimeters)) # Reduced height for better spacing +map_item.zoomToExtent(st_b10_layer.extent()) +layout.addLayoutItem(map_item) + +# Add Coordinate Grid +grid = map_item.grid() +grid.setEnabled(True) +grid.setIntervalX(100000) # Set interval for grid lines in map units +grid.setIntervalY(100000) +grid.setAnnotationEnabled(True) +grid.setAnnotationPrecision(0) # Set annotation precision for grid labels +grid.setFrameStyle(QgsLayoutItemMapGrid.Zebra) # Optional grid style + +# Label +title_label = QgsLayoutItemLabel(layout) +title_label.setText("Surace Temperautre B11 - Lab 2") +title_label.setFont(QFont("Arial", 16)) +title_label.setHAlign(Qt.AlignCenter) +title_label.attemptMove(QgsLayoutPoint(105, 10, QgsUnitTypes.LayoutMillimeters)) +title_label.adjustSizeToText() +layout.addLayoutItem(title_label) + +# Legend +legend_item = QgsLayoutItemLegend(layout) +legend_item.setLinkedMap(map_item) +legend_item.setTitle("Legend") +legend_item.setFrameEnabled(True) +legend_item.attemptMove(QgsLayoutPoint(245, 30, QgsUnitTypes.LayoutMillimeters)) # Adjusted for better alignment and spacing +layout.addLayoutItem(legend_item) + +# Scalebar +scalebar_item = QgsLayoutItemScaleBar(layout) +scalebar_item.setStyle('Single Box') # Change to 'Single Box' or 'Double Box' for a proper scalebar +scalebar_item.setLinkedMap(map_item) +scalebar_item.setUnitLabel('m') # Add unit label to make it informative +scalebar_item.setNumberOfSegments(4) # Specify the number of segments in the scalebar +scalebar_item.setNumberOfSegmentsLeft(0) # Segments to the left of zero (if any) +scalebar_item.setUnitsPerSegment(50000) # Set distance per segment in map units +scalebar_item.setFont(QFont('Arial', 10)) # Set font size for better readability +scalebar_item.setHeight(5) # Set the height of the scalebar +scalebar_item.applyDefaultSize() +scalebar_item.attemptMove(QgsLayoutPoint(20, 190, QgsUnitTypes.LayoutMillimeters)) # Moved below the map for better organization +layout.addLayoutItem(scalebar_item) + +# North Arrow +north_arrow_item = QgsLayoutItemPicture(layout) +north_arrow_item.setPicturePath('/home/huaqo/dev/courses/2024_Remote_Sensing/styles/north_arrow.svg') +if not north_arrow_item.picturePath(): + print("North arrow image failed to load!") +north_arrow_item.setReferencePoint(QgsLayoutItemPicture.UpperLeft) +north_arrow_item.attemptMove(QgsLayoutPoint(30, 40, QgsUnitTypes.LayoutMillimeters)) # Adjusted placement for better visibility +north_arrow_item.attemptResize(QgsLayoutSize(15, 15, QgsUnitTypes.LayoutMillimeters)) # Increased size for better visibility +layout.addLayoutItem(north_arrow_item) + +# Author +info_label = QgsLayoutItemLabel(layout) +info_label.setText("Author: Joaquin Gottlebe") +info_label.setFont(QFont("Arial", 12)) +info_label.setHAlign(Qt.AlignLeft) +info_label.attemptMove(QgsLayoutPoint(245, 100, QgsUnitTypes.LayoutMillimeters)) +info_label.adjustSizeToText() +layout.addLayoutItem(info_label) + +# Author +info_label = QgsLayoutItemLabel(layout) +info_label.setText("Date: 16.11.2025") +info_label.setFont(QFont("Arial", 12)) +info_label.setHAlign(Qt.AlignLeft) +info_label.attemptMove(QgsLayoutPoint(245, 110, QgsUnitTypes.LayoutMillimeters)) +info_label.adjustSizeToText() +layout.addLayoutItem(info_label) + +# Coordinate System +info_label = QgsLayoutItemLabel(layout) +info_label.setText("CRS: EPSG 32633") +info_label.setFont(QFont("Arial", 12)) +info_label.setHAlign(Qt.AlignLeft) +info_label.attemptMove(QgsLayoutPoint(245, 120, QgsUnitTypes.LayoutMillimeters)) +info_label.adjustSizeToText() +layout.addLayoutItem(info_label) + + +# Export +export_path = "/home/huaqo/dev/courses/2024_Remote_Sensing/lab2/maps/surface_temp_b11.png" +exporter = QgsLayoutExporter(layout) +export_result = exporter.exportToImage(export_path, QgsLayoutExporter.ImageExportSettings()) +if export_result != QgsLayoutExporter.Success: + print("Failed to export map!") +else: + print("Map exported to:", export_path) diff --git a/2024_Remote_Sensing/lab2/toa.py b/2024_Remote_Sensing/lab2/toa.py new file mode 100644 index 0000000..1450540 --- /dev/null +++ b/2024_Remote_Sensing/lab2/toa.py @@ -0,0 +1,122 @@ +from osgeo import gdal +import numpy as np +import sys +import os +import re +import math + +def tif_to_array(path): + data = gdal.Open(path) + if data is None: + raise FileNotFoundError(f"Cannot open file: {path}") + band = data.GetRasterBand(1) + array = band.ReadAsArray().astype(np.float32) + return data, array + +def handle_nodata(array): + # Replace only negative values with NaN + clean_array = np.where(array < 0, np.nan, array) + return clean_array + +def reflectance_calc(array, reflectance_mult, reflectance_add, sun_elevation): + # Debugging print statements + print(f"Reflectance Mult: {reflectance_mult}, Reflectance Add: {reflectance_add}, Sun Elevation: {sun_elevation}") + print(f"Array Sample Before: {array[:5, :5]}") # Print a small portion of the array for inspection + + result = ((array * reflectance_mult) + reflectance_add) / math.sin(math.radians(sun_elevation)) + + # Debugging print statements + print(f"Array Sample After: {result[:5, :5]}") # Print a small portion of the result for inspection + return result + +def save_tif(array, data, output_path): + driver = gdal.GetDriverByName('GTiff') + rows, cols = array.shape + out_data = driver.Create(output_path, cols, rows, 1, gdal.GDT_Float32) + out_data.SetGeoTransform(data.GetGeoTransform()) + out_data.SetProjection(data.GetProjection()) + out_band = out_data.GetRasterBand(1) + out_band.WriteArray(array) + out_band.SetNoDataValue(-9999) + out_band.FlushCache() + out_data = None + +def get_output_path(output_dir, input_path, suffix): + basename = os.path.basename(input_path) + basename_no_ext = os.path.splitext(basename)[0] +# basename_no_ext = re.sub(r'_B[1-9]|_B10|_B11', '', basename_no_ext) + path = os.path.join(output_dir, f"{basename_no_ext}_{suffix}.TIF") + return path + +def print_min_max(array, label): + print(label, " Min:", np.nanmin(array), " Max:", np.nanmax(array)) + +def get_reflectance_params(band_number, metadata_lines): + reflectance_mult = None + reflectance_add = None + for line in metadata_lines: + if f"REFLECTANCE_MULT_BAND_{band_number}" in line: + variable, value = line.split(" = ") + reflectance_mult = float(value.strip()) + if f"REFLECTANCE_ADD_BAND_{band_number}" in line: + variable, value = line.split(" = ") + reflectance_add = float(value.strip()) + if reflectance_mult is None or reflectance_add is None: + raise ValueError(f"Reflectance parameters not found for band {band_number} in metadata.") + return reflectance_mult, reflectance_add + +def process_bands(band_paths, metadata_path, output_dir): + # Read metadata + with open(metadata_path, "r") as file: + metadata_lines = file.readlines() + + # Extract sun elevation + sun_elevation = None + for line in metadata_lines: + if "SUN_ELEVATION" in line: + variable, value = line.split(" = ") + sun_elevation = float(value.strip()) + break + + if sun_elevation is None: + raise ValueError("Sun elevation not found in metadata.") + + # Debugging sun elevation value + print(f"Sun Elevation: {sun_elevation}, Sun Elevation Radians: {math.radians(sun_elevation)}") + + # Process each band + for band_path in band_paths: + band_number = int(re.search(r'_B(\d+)', band_path).group(1)) + data, array = tif_to_array(band_path) + array = handle_nodata(array) + + # Get reflectance parameters for the band + reflectance_mult, reflectance_add = get_reflectance_params(band_number, metadata_lines) + + # Calculate reflectance for the band + reflectance = reflectance_calc(array, reflectance_mult, reflectance_add, sun_elevation) + + # Create output directory if it doesn't exist + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + # Save Top of Atmosphere (TOA) reflectance for the band + output_path = get_output_path(output_dir, band_path, f"toa") + save_tif(reflectance, data, output_path) + + # Print min and max for the reflectance band + print_min_max(reflectance, f"Band {band_number} Reflectance") + +def main(): + if len(sys.argv) < 4: + print("Usage: python3 script.py ") + sys.exit(1) + + metadata_path = sys.argv[1] + output_dir = sys.argv[2] + band_paths = sys.argv[3:] + + process_bands(band_paths, metadata_path, output_dir) + +if __name__ == "__main__": + main() diff --git a/2024_Remote_Sensing/lab2/toa_radiance.py b/2024_Remote_Sensing/lab2/toa_radiance.py new file mode 100644 index 0000000..b3ec38b --- /dev/null +++ b/2024_Remote_Sensing/lab2/toa_radiance.py @@ -0,0 +1,107 @@ +from osgeo import gdal +import numpy as np +import sys +import os +import re +import math + +def tif_to_array(path): + data = gdal.Open(path) + if data is None: + raise FileNotFoundError(f"Cannot open file: {path}") + band = data.GetRasterBand(1) + array = band.ReadAsArray().astype(np.float32) + return data, array + +def get_radiance_params(band_number, metadata_lines): + radiance_mult = None + radiance_add = None + for line in metadata_lines: + if f"RADIANCE_MULT_BAND_{band_number}" in line: + _, value = line.split(" = ") + radiance_mult = float(value.strip()) + if f"RADIANCE_ADD_BAND_{band_number}" in line: + _, value = line.split(" = ") + radiance_add = float(value.strip()) + if radiance_mult is None or radiance_add is None: + raise ValueError(f"Radiance parameters not found for band {band_number} in metadata.") + return radiance_mult, radiance_add + +def save_tif(array, data, output_path): + driver = gdal.GetDriverByName('GTiff') + rows, cols = array.shape + out_data = driver.Create(output_path, cols, rows, 1, gdal.GDT_Float32) + out_data.SetGeoTransform(data.GetGeoTransform()) + out_data.SetProjection(data.GetProjection()) + out_band = out_data.GetRasterBand(1) + out_band.WriteArray(array) + out_band.SetNoDataValue(-9999) + out_band.FlushCache() + out_data = None + +def get_output_path(output_dir, input_path, suffix): + basename = os.path.basename(input_path) + basename_no_ext = os.path.splitext(basename)[0] + path = os.path.join(output_dir, f"{basename_no_ext}_{suffix}.TIF") + return path + +def print_min_max(array, label): + print(label, " Min:", np.nanmin(array), " Max:", np.nanmax(array)) + +def process_bands(bands, metadata_path, output_dir): + # Read metadata + with open(metadata_path, "r") as file: + metadata_lines = file.readlines() + + # Extract sun elevation (optional for radiance, included for debugging) + sun_elevation = None + for line in metadata_lines: + if "SUN_ELEVATION" in line: + _, value = line.split(" = ") + sun_elevation = float(value.strip()) + break + + # Debugging sun elevation value + if sun_elevation is not None: + print(f"Sun Elevation: {sun_elevation}, Sun Elevation Radians: {math.radians(sun_elevation)}") + + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + # Process each band + for band_data in bands: + band_number, band_path = band_data[0], band_data[1] + + # Load the band data + data, array = tif_to_array(band_path) + + # Get radiance parameters for the band + radiance_mult, radiance_add = get_radiance_params(band_number, metadata_lines) + + # Calculate radiance + radiance_array = (array * radiance_mult) + radiance_add + + # Save the Top of Atmosphere (TOA) radiance + output_path = get_output_path(output_dir, band_path, "toa_radiance") + save_tif(radiance_array, data, output_path) + + # Print min and max for the radiance band + print_min_max(radiance_array, f"Band {band_number} Radiance") + +def main(): + print(f"Number of arguments received: {len(sys.argv) - 1}") + print("Arguments:", sys.argv) + if len(sys.argv) < 4: + print("Usage: python3 script.py ") + sys.exit(1) + + metadata_path = sys.argv[1] + output_dir = sys.argv[2] + band_paths = sys.argv[3:] + + bands = [(int(re.search(r'B(\d+)', path).group(1)), path) for path in band_paths] + + process_bands(bands, metadata_path, output_dir) + +if __name__ == "__main__": + main() diff --git a/2024_Remote_Sensing/lab2/true_color_map.py b/2024_Remote_Sensing/lab2/true_color_map.py new file mode 100644 index 0000000..5e6986d --- /dev/null +++ b/2024_Remote_Sensing/lab2/true_color_map.py @@ -0,0 +1,119 @@ +# Cleanup +project = QgsProject.instance() +layer_ids = list(project.mapLayers().keys()) +for layer_id in layer_ids: + project.removeMapLayer(layer_id) + +# Loading B10 Layer +path_st_b10 = "/home/huaqo/dev/courses/2024_Remote_Sensing/lab2/true_color/true_color.tif" +st_b10_layer = QgsRasterLayer(path_st_b10, "True Color Composite") +if not st_b10_layer.isValid(): + print("st_b10 layer failed to load!") +else: + QgsProject.instance().addMapLayer(st_b10_layer) + print("st_b10 Layer loaded!") + +# Styling Layer +st_b10_style_path = '/home/huaqo/dev/courses/2024_Remote_Sensing/styles/true_color.qml' +if not st_b10_layer.loadNamedStyle(st_b10_style_path): + print("Failed to load st_b10 style!") +else: + print("st_b10 style loaded!") + +# Load layout +layout = QgsPrintLayout(project) +layout.initializeDefaults() + +# Map +map_item = QgsLayoutItemMap(layout) +map_item.attemptMove(QgsLayoutPoint(20, 30, QgsUnitTypes.LayoutMillimeters)) # Added space between top and map +map_item.attemptResize(QgsLayoutSize(200, 150, QgsUnitTypes.LayoutMillimeters)) # Reduced height for better spacing +map_item.zoomToExtent(st_b10_layer.extent()) +layout.addLayoutItem(map_item) + +# Add Coordinate Grid +grid = map_item.grid() +grid.setEnabled(True) +grid.setIntervalX(100000) # Set interval for grid lines in map units +grid.setIntervalY(100000) +grid.setAnnotationEnabled(True) +grid.setAnnotationPrecision(0) # Set annotation precision for grid labels +grid.setFrameStyle(QgsLayoutItemMapGrid.Zebra) # Optional grid style + +# Label +title_label = QgsLayoutItemLabel(layout) +title_label.setText("True Color Composite - Lab 2") +title_label.setFont(QFont("Arial", 16)) +title_label.setHAlign(Qt.AlignCenter) +title_label.attemptMove(QgsLayoutPoint(105, 10, QgsUnitTypes.LayoutMillimeters)) +title_label.adjustSizeToText() +layout.addLayoutItem(title_label) + +# Legend +legend_item = QgsLayoutItemLegend(layout) +legend_item.setLinkedMap(map_item) +legend_item.setTitle("Legend") +legend_item.setFrameEnabled(True) +legend_item.attemptMove(QgsLayoutPoint(245, 30, QgsUnitTypes.LayoutMillimeters)) # Adjusted for better alignment and spacing +layout.addLayoutItem(legend_item) + +# Scalebar +scalebar_item = QgsLayoutItemScaleBar(layout) +scalebar_item.setStyle('Single Box') # Change to 'Single Box' or 'Double Box' for a proper scalebar +scalebar_item.setLinkedMap(map_item) +scalebar_item.setUnitLabel('m') # Add unit label to make it informative +scalebar_item.setNumberOfSegments(4) # Specify the number of segments in the scalebar +scalebar_item.setNumberOfSegmentsLeft(0) # Segments to the left of zero (if any) +scalebar_item.setUnitsPerSegment(50000) # Set distance per segment in map units +scalebar_item.setFont(QFont('Arial', 10)) # Set font size for better readability +scalebar_item.setHeight(5) # Set the height of the scalebar +scalebar_item.applyDefaultSize() +scalebar_item.attemptMove(QgsLayoutPoint(20, 190, QgsUnitTypes.LayoutMillimeters)) # Moved below the map for better organization +layout.addLayoutItem(scalebar_item) + +# North Arrow +north_arrow_item = QgsLayoutItemPicture(layout) +north_arrow_item.setPicturePath('/home/huaqo/dev/courses/2024_Remote_Sensing/styles/north_arrow.svg') +if not north_arrow_item.picturePath(): + print("North arrow image failed to load!") +north_arrow_item.setReferencePoint(QgsLayoutItemPicture.UpperLeft) +north_arrow_item.attemptMove(QgsLayoutPoint(30, 40, QgsUnitTypes.LayoutMillimeters)) # Adjusted placement for better visibility +north_arrow_item.attemptResize(QgsLayoutSize(15, 15, QgsUnitTypes.LayoutMillimeters)) # Increased size for better visibility +layout.addLayoutItem(north_arrow_item) + +# Author +info_label = QgsLayoutItemLabel(layout) +info_label.setText("Author: Joaquin Gottlebe") +info_label.setFont(QFont("Arial", 12)) +info_label.setHAlign(Qt.AlignLeft) +info_label.attemptMove(QgsLayoutPoint(245, 100, QgsUnitTypes.LayoutMillimeters)) +info_label.adjustSizeToText() +layout.addLayoutItem(info_label) + +# Author +info_label = QgsLayoutItemLabel(layout) +info_label.setText("Date: 16.11.2025") +info_label.setFont(QFont("Arial", 12)) +info_label.setHAlign(Qt.AlignLeft) +info_label.attemptMove(QgsLayoutPoint(245, 110, QgsUnitTypes.LayoutMillimeters)) +info_label.adjustSizeToText() +layout.addLayoutItem(info_label) + +# Coordinate System +info_label = QgsLayoutItemLabel(layout) +info_label.setText("CRS: EPSG 32633") +info_label.setFont(QFont("Arial", 12)) +info_label.setHAlign(Qt.AlignLeft) +info_label.attemptMove(QgsLayoutPoint(245, 120, QgsUnitTypes.LayoutMillimeters)) +info_label.adjustSizeToText() +layout.addLayoutItem(info_label) + + +# Export +export_path = "/home/huaqo/dev/courses/2024_Remote_Sensing/lab2/maps/true_color.png" +exporter = QgsLayoutExporter(layout) +export_result = exporter.exportToImage(export_path, QgsLayoutExporter.ImageExportSettings()) +if export_result != QgsLayoutExporter.Success: + print("Failed to export map!") +else: + print("Map exported to:", export_path) diff --git a/2024_Remote_Sensing/lab3/Lab3.md b/2024_Remote_Sensing/lab3/Lab3.md new file mode 100644 index 0000000..667ea89 --- /dev/null +++ b/2024_Remote_Sensing/lab3/Lab3.md @@ -0,0 +1,27 @@ +# Lab 3 + +Name: Joaquin Gottlebe + +Matrikelnummer: 829101 + +## Question 1 + +![](maps/ndwi.png) + +## Question 2 + +![](maps/diff_ndwi.png) + +![](plots/diff_hist.png) + +Landsat 8 uses broader NIR and Green bands with a 30m resolution, while Sentinel-2 uses narrower bands with a 10m resolution. Differences in atmospheric correction methods and time of capture could also be factors. + +## Question 3 + +For the third question i choose the Chacao Channel as it is also a very interesting region because of its national parks agriculture and abundance of waterbodies. + +![](maps/ndvi2.png) + +![](maps/ndwi2.png) + +Here the symbology was reduced to 0-0.2 because the calculations resultet in a peak around this area and near to zero values till 1. Which doenst make much sense because the area has a lot of waterbodies. I couldnt find and error. \ No newline at end of file diff --git a/2024_Remote_Sensing/lab3/Makefile b/2024_Remote_Sensing/lab3/Makefile new file mode 100644 index 0000000..cd598bd --- /dev/null +++ b/2024_Remote_Sensing/lab3/Makefile @@ -0,0 +1,44 @@ +.PHONY: all preproc1 preproc2 proc proc2 maps export + +ACTIVATE_SCRIPT = ~/miniforge3/bin/activate +CONDA_ENV_NAME = gdal_env +PANDOC=pandoc + +all: preproc1 preproc2 proc proc2 maps export + +preproc1: + gdal_translate -of GTiff /home/huaqo/dev/courses/2024_Remote_Sensing/lab3/S2A_MSIL2A_20170527T102031_N9999_R065_T33UUU_20191015T100203.SAFE/GRANULE/L2A_T33UUU_A010072_20170527T102301/IMG_DATA/R10m/T33UUU_20170527T102031_B02_10m.jp2 /home/huaqo/dev/courses/2024_Remote_Sensing/lab3/sen2tiff/B02.tif && \ + gdal_translate -of GTiff /home/huaqo/dev/courses/2024_Remote_Sensing/lab3/S2A_MSIL2A_20170527T102031_N9999_R065_T33UUU_20191015T100203.SAFE/GRANULE/L2A_T33UUU_A010072_20170527T102301/IMG_DATA/R10m/T33UUU_20170527T102031_B03_10m.jp2 /home/huaqo/dev/courses/2024_Remote_Sensing/lab3/sen2tiff/B03.tif && \ + gdal_translate -of GTiff /home/huaqo/dev/courses/2024_Remote_Sensing/lab3/S2A_MSIL2A_20170527T102031_N9999_R065_T33UUU_20191015T100203.SAFE/GRANULE/L2A_T33UUU_A010072_20170527T102301/IMG_DATA/R10m/T33UUU_20170527T102031_B04_10m.jp2 /home/huaqo/dev/courses/2024_Remote_Sensing/lab3/sen2tiff/B04.tif && \ + gdal_translate -of GTiff /home/huaqo/dev/courses/2024_Remote_Sensing/lab3/S2A_MSIL2A_20170527T102031_N9999_R065_T33UUU_20191015T100203.SAFE/GRANULE/L2A_T33UUU_A010072_20170527T102301/IMG_DATA/R10m/T33UUU_20170527T102031_B08_10m.jp2 /home/huaqo/dev/courses/2024_Remote_Sensing/lab3/sen2tiff/B08.tif + +preproc2: + gdal_translate -of GTiff /home/huaqo/dev/courses/2024_Remote_Sensing/lab3/S2A_MSIL2A_20170527T102031_N9999_R065_T33UUU_20191015T100203.SAFE/GRANULE/L2A_T33UUU_A010072_20170527T102301/IMG_DATA/R10m/T33UUU_20170527T102031_B02_10m.jp2 /home/huaqo/dev/courses/2024_Remote_Sensing/lab3/sen2tiff2/B02.tif && \ + gdal_translate -of GTiff /home/huaqo/dev/courses/2024_Remote_Sensing/lab3/S2A_MSIL2A_20241118T143741_N0511_R096_T18GXU_20241118T201252.SAFE/GRANULE/L2A_T18GXU_A049142_20241118T144600/IMG_DATA/R10m/T18GXU_20241118T143741_B03_10m.jp2 /home/huaqo/dev/courses/2024_Remote_Sensing/lab3/sen2tiff2/B03.tif && \ + gdal_translate -of GTiff /home/huaqo/dev/courses/2024_Remote_Sensing/lab3/S2A_MSIL2A_20241118T143741_N0511_R096_T18GXU_20241118T201252.SAFE/GRANULE/L2A_T18GXU_A049142_20241118T144600/IMG_DATA/R10m/T18GXU_20241118T143741_B04_10m.jp2 /home/huaqo/dev/courses/2024_Remote_Sensing/lab3/sen2tiff2/B04.tif && \ + gdal_translate -of GTiff /home/huaqo/dev/courses/2024_Remote_Sensing/lab3/S2A_MSIL2A_20241118T143741_N0511_R096_T18GXU_20241118T201252.SAFE/GRANULE/L2A_T18GXU_A049142_20241118T144600/IMG_DATA/R10m/T18GXU_20241118T143741_B08_10m.jp2 /home/huaqo/dev/courses/2024_Remote_Sensing/lab3/sen2tiff2/B08.tif +proc: + @source $(ACTIVATE_SCRIPT) && conda activate $(CONDA_ENV_NAME) && \ + python3 ndvi.py sen2tiff/B04.tif sen2tiff/B08.tif ndvi && \ + python3 ndwi.py sen2tiff/B03.tif sen2tiff/B08.tif ndwi && \ + python clip.py /home/huaqo/dev/courses/2024_Remote_Sensing/lab3/l8_ndvi/l8_ndvi.TIF /home/huaqo/dev/courses/2024_Remote_Sensing/lab3/ndvi/B04_ndvi.TIF /home/huaqo/dev/courses/2024_Remote_Sensing/lab3/clip/l8_ndvi.TIF /home/huaqo/dev/courses/2024_Remote_Sensing/lab3/clip/sen_ndvi.TIF && \ + python diff.py /home/huaqo/dev/courses/2024_Remote_Sensing/lab3/clip/sen_ndvi.TIF /home/huaqo/dev/courses/2024_Remote_Sensing/lab3/clip/l8_ndvi.TIF /home/huaqo/dev/courses/2024_Remote_Sensing/lab3/diff/aligned.TIF /home/huaqo/dev/courses/2024_Remote_Sensing/lab3/diff/diff.TIF && \ + python clip.py /home/huaqo/dev/courses/2024_Remote_Sensing/lab3/l8_ndwi/l8_ndwi.TIF /home/huaqo/dev/courses/2024_Remote_Sensing/lab3/ndwi/B03_ndwi.TIF /home/huaqo/dev/courses/2024_Remote_Sensing/lab3/clip/l8_ndwi.TIF /home/huaqo/dev/courses/2024_Remote_Sensing/lab3/clip/sen_ndwi.TIF && \ + python diff.py /home/huaqo/dev/courses/2024_Remote_Sensing/lab3/clip/sen_ndwi.TIF /home/huaqo/dev/courses/2024_Remote_Sensing/lab3/clip/l8_ndwi.TIF /home/huaqo/dev/courses/2024_Remote_Sensing/lab3/diff/aligned_ndwi.TIF /home/huaqo/dev/courses/2024_Remote_Sensing/lab3/diff/diff_ndwi.TIF && \ + python hist.py diff/diff_ndwi.TIF --output plots/diff_hist.png + +proc2: + @source $(ACTIVATE_SCRIPT) && conda activate $(CONDA_ENV_NAME) && \ + python3 ndvi.py sen2tiff2/B04.tif sen2tiff2/B08.tif ndvi2 && \ + python3 ndwi.py sen2tiff2/B03.tif sen2tiff2/B08.tif ndwi2 + +maps: + python3 ndvi_map.py + python3 ndvi_map2.py + python3 ndwi_map.py + python3 ndwi_map2.py + python3 diff_map.py + +export: + + $(PANDOC) Lab3.md -o rcm01_2425_lab3_gottlebe_829101.pdf \ No newline at end of file diff --git a/2024_Remote_Sensing/lab3/clip.py b/2024_Remote_Sensing/lab3/clip.py new file mode 100644 index 0000000..41c53b4 --- /dev/null +++ b/2024_Remote_Sensing/lab3/clip.py @@ -0,0 +1,62 @@ +import rasterio +from rasterio.mask import mask +from shapely.geometry import box +import argparse + +def get_common_extent(raster1_path, raster2_path): + """Calculate the common extent (intersection) of two rasters.""" + with rasterio.open(raster1_path) as src1, rasterio.open(raster2_path) as src2: + bounds1 = src1.bounds + bounds2 = src2.bounds + + # Calculate the intersection of the extents + common_bounds = box( + max(bounds1.left, bounds2.left), + max(bounds1.bottom, bounds2.bottom), + min(bounds1.right, bounds2.right), + min(bounds1.top, bounds2.top), + ) + return common_bounds + +def clip_raster_to_extent(raster_path, extent, output_path): + """Clip the raster to the provided extent and save it to a new file.""" + with rasterio.open(raster_path) as src: + # Convert extent to GeoJSON-like dict for rasterio.mask + geojson_extent = [extent.__geo_interface__] + + # Clip the raster to the extent + out_image, out_transform = mask(src, geojson_extent, crop=True) + out_meta = src.meta.copy() + + # Update metadata with new dimensions, transform, and bounds + out_meta.update({ + "driver": "GTiff", + "height": out_image.shape[1], + "width": out_image.shape[2], + "transform": out_transform + }) + + # Save the clipped raster + with rasterio.open(output_path, "w", **out_meta) as dest: + dest.write(out_image) + +def main(): + parser = argparse.ArgumentParser(description="Clip two rasters to their common extent.") + parser.add_argument("raster1", help="Path to the first raster.") + parser.add_argument("raster2", help="Path to the second raster.") + parser.add_argument("output1", help="Path to save the clipped first raster.") + parser.add_argument("output2", help="Path to save the clipped second raster.") + + args = parser.parse_args() + + # Determine the common extent + common_extent = get_common_extent(args.raster1, args.raster2) + + # Clip both rasters to the common extent + clip_raster_to_extent(args.raster1, common_extent, args.output1) + clip_raster_to_extent(args.raster2, common_extent, args.output2) + + print("Rasters clipped successfully!") + +if __name__ == "__main__": + main() diff --git a/2024_Remote_Sensing/lab3/diff.py b/2024_Remote_Sensing/lab3/diff.py new file mode 100644 index 0000000..3d4db05 --- /dev/null +++ b/2024_Remote_Sensing/lab3/diff.py @@ -0,0 +1,85 @@ +import rasterio +from rasterio.warp import reproject, Resampling +import numpy as np +import argparse + +def resample_raster_to_match(reference_path, target_path, output_path): + """Resample the target raster to match the resolution, extent, and CRS of the reference raster.""" + with rasterio.open(reference_path) as ref: + ref_transform = ref.transform + ref_crs = ref.crs + ref_width = ref.width + ref_height = ref.height + ref_nodata = ref.nodata or -9999 + + with rasterio.open(target_path) as target: + target_nodata = target.nodata or -9999 + profile = target.profile + profile.update( + transform=ref_transform, + crs=ref_crs, + width=ref_width, + height=ref_height, + nodata=ref_nodata + ) + + # Resample target raster + with rasterio.open(output_path, "w", **profile) as output: + for i in range(1, target.count + 1): + reproject( + source=rasterio.band(target, i), + destination=rasterio.band(output, i), + src_transform=target.transform, + src_crs=target.crs, + dst_transform=ref_transform, + dst_crs=ref_crs, + resampling=Resampling.bilinear, + ) + +def calculate_difference(raster1_path, raster2_resampled_path, output_path): + """Calculate the difference between two aligned rasters.""" + with rasterio.open(raster1_path) as src1: + raster1 = src1.read(1) + nodata1 = src1.nodata or -9999 + profile = src1.profile + + with rasterio.open(raster2_resampled_path) as src2: + raster2 = src2.read(1) + nodata2 = src2.nodata or -9999 + + # Replace nodata values with NaN for calculation + raster1 = np.where(raster1 == nodata1, np.nan, raster1) + raster2 = np.where(raster2 == nodata2, np.nan, raster2) + + # Compute difference + difference = raster1 - raster2 + + # Replace NaN with nodata for saving + difference[np.isnan(difference)] = -9999 + + # Update metadata to include nodata + profile.update(nodata=-9999) + + # Save the difference raster + with rasterio.open(output_path, "w", **profile) as dst: + dst.write(difference, 1) + +def main(): + parser = argparse.ArgumentParser(description="Resample a raster and calculate the difference.") + parser.add_argument("reference_raster", help="Path to the reference raster.") + parser.add_argument("target_raster", help="Path to the target raster to be resampled.") + parser.add_argument("aligned_raster", help="Path to save the resampled (aligned) raster.") + parser.add_argument("difference_raster", help="Path to save the difference raster.") + + args = parser.parse_args() + + # Resample raster 2 to match raster 1 + resample_raster_to_match(args.reference_raster, args.target_raster, args.aligned_raster) + + # Calculate the difference + calculate_difference(args.reference_raster, args.aligned_raster, args.difference_raster) + + print("Aligned raster and difference raster created successfully!") + +if __name__ == "__main__": + main() diff --git a/2024_Remote_Sensing/lab3/diff_map.py b/2024_Remote_Sensing/lab3/diff_map.py new file mode 100644 index 0000000..78c2891 --- /dev/null +++ b/2024_Remote_Sensing/lab3/diff_map.py @@ -0,0 +1,171 @@ +from qgis.core import * +from PyQt5.QtGui import * +from PyQt5.QtCore import * +from PyQt5 import * +from datetime import datetime +import os + +QGIS_PREFIX_PATH = "/usr/share/qgis" +LAYER_PATH = "/home/huaqo/dev/courses/2024_Remote_Sensing/lab3/diff/diff_ndwi.TIF" +STYLE_PATH = "/home/huaqo/dev/courses/2024_Remote_Sensing/styles/diff.qml" +EXPORT_PATH = "/home/huaqo/dev/courses/2024_Remote_Sensing/lab3/maps/diff_ndwi.png" +NORTH_ARROW_PATH = "/home/huaqo/dev/courses/2024_Remote_Sensing/styles/north_arrow.svg" + +map_pos_x = 20 +map_pos_y = 30 +map_size_x = 200 +map_size_y = 150 + +grid_interval_x = 100000 +grid_interval_y = 100000 + +LAYER_NAME = "DIFF" +legend_pos_x = 245 +legend_pos_y = 30 + +scalebar_pos_x = 20 +scalebar_pos_y = 190 + +north_arrow_pos_x = 30 +north_arrow_pos_y = 40 +north_arrow_size_x = 15 +north_arrow_size_y = 15 + +title_text = "Difference NDWI Landsat 8 and Sentinel 2 - Lab 3" +title_fontsize = 16 +title_height = 10 + +author_text = "Author: Joaquin Gottlebe" +author_pos_y = 100 + +date_text = f"Date: {datetime.now().strftime('%d.%m.%Y')}" +date_pos_y = 110 + +crs_text = "CRS: EPSG 32633" +crs_pos_y = 120 + +fontfamily = "Arial" +info_font_size = 12 +info_pos_x = legend_pos_x + +# Initialize QGIS Application +QgsApplication.setPrefixPath(QGIS_PREFIX_PATH, True) +qgs = QgsApplication([], False) +qgs.initQgis() + +# Cleanup +project = QgsProject.instance() +layer_ids = list(project.mapLayers().keys()) +for layer_id in layer_ids: + project.removeMapLayer(layer_id) + +# Load Layer +raster_layer = QgsRasterLayer(LAYER_PATH, LAYER_NAME) +if not raster_layer.isValid(): + raise Exception(f"{LAYER_NAME} layer failed to load!") + +QgsProject.instance().addMapLayer(raster_layer) +if not raster_layer.loadNamedStyle(STYLE_PATH): + raise Exception(f"Style {STYLE_PATH} failed to load") + +# Create Layout +layout = QgsPrintLayout(project) +layout.initializeDefaults() + +# Map +map_item = QgsLayoutItemMap(layout) +map_item.attemptMove(QgsLayoutPoint(map_pos_x,map_pos_y, QgsUnitTypes.LayoutMillimeters)) +map_item.attemptResize(QgsLayoutSize(map_size_x, map_size_y, QgsUnitTypes.LayoutMillimeters)) +map_item.zoomToExtent(raster_layer.extent()) +layout.addLayoutItem(map_item) + +# Add Coordinate Grid +grid = map_item.grid() +grid.setEnabled(True) +grid.setIntervalX(grid_interval_x) +grid.setIntervalY(grid_interval_y) +grid.setAnnotationEnabled(True) +grid.setAnnotationPrecision(0) +grid.setFrameStyle(QgsLayoutItemMapGrid.Zebra) + +# Label +title_label = QgsLayoutItemLabel(layout) +title_label.setText(title_text) +title_label.setFont(QFont(fontfamily, title_fontsize)) +title_label.setHAlign(Qt.AlignCenter) +title_label.adjustSizeToText() +layout_width = layout.pageCollection().page(0).pageSize().width() +title_width = title_label.rectWithFrame().width() +center_x = (layout_width - title_width) / 2 +title_label.attemptMove(QgsLayoutPoint(center_x, title_height, QgsUnitTypes.LayoutMillimeters)) +layout.addLayoutItem(title_label) + +# Legend +legend_item = QgsLayoutItemLegend(layout) +legend_item.setLinkedMap(map_item) +legend_item.setTitle("Legend") +legend_item.setFrameEnabled(True) +legend_item.attemptMove(QgsLayoutPoint(legend_pos_x,legend_pos_y, QgsUnitTypes.LayoutMillimeters)) +layout.addLayoutItem(legend_item) + +# Scalebar +scalebar_item = QgsLayoutItemScaleBar(layout) +scalebar_item.setStyle('Single Box') +scalebar_item.setLinkedMap(map_item) +scalebar_item.setUnitLabel('m') +scalebar_item.setNumberOfSegments(4) +scalebar_item.setNumberOfSegmentsLeft(0) +scalebar_item.setUnitsPerSegment(50000) +scalebar_item.applyDefaultSize() +scalebar_item.attemptMove(QgsLayoutPoint(scalebar_pos_x,scalebar_pos_y, QgsUnitTypes.LayoutMillimeters)) +layout.addLayoutItem(scalebar_item) + +# North Arrow +north_arrow_item = QgsLayoutItemPicture(layout) +north_arrow_item.setPicturePath(NORTH_ARROW_PATH) +if not north_arrow_item.picturePath(): + print("North arrow image failed to load!") +north_arrow_item.setReferencePoint(QgsLayoutItemPicture.UpperLeft) +north_arrow_item.attemptMove(QgsLayoutPoint(north_arrow_pos_x,north_arrow_pos_y,QgsUnitTypes.LayoutMillimeters)) +north_arrow_item.attemptResize(QgsLayoutSize(north_arrow_size_x,north_arrow_size_y,QgsUnitTypes.LayoutMillimeters)) +layout.addLayoutItem(north_arrow_item) + +# Author Info +info_label = QgsLayoutItemLabel(layout) +info_label.setText(author_text) +info_label.setFont(QFont(fontfamily, info_font_size)) +info_label.setHAlign(Qt.AlignLeft) +info_label.attemptMove(QgsLayoutPoint(info_pos_x, author_pos_y, QgsUnitTypes.LayoutMillimeters)) +info_label.adjustSizeToText() +layout.addLayoutItem(info_label) + +# Date Info +info_label = QgsLayoutItemLabel(layout) +info_label.setText(date_text) +info_label.setFont(QFont(fontfamily, info_font_size)) +info_label.setHAlign(Qt.AlignLeft) +info_label.attemptMove(QgsLayoutPoint(info_pos_x, date_pos_y, QgsUnitTypes.LayoutMillimeters)) +info_label.adjustSizeToText() +layout.addLayoutItem(info_label) + +# CRS Info +info_label = QgsLayoutItemLabel(layout) +info_label.setText(crs_text) +info_label.setFont(QFont(fontfamily, info_font_size)) +info_label.setHAlign(Qt.AlignLeft) +info_label.attemptMove(QgsLayoutPoint(info_pos_x, crs_pos_y, QgsUnitTypes.LayoutMillimeters)) +info_label.adjustSizeToText() +layout.addLayoutItem(info_label) + +# Export Map +exporter = QgsLayoutExporter(layout) +if os.path.exists(EXPORT_PATH): + os.remove(EXPORT_PATH) +export_result = exporter.exportToImage(EXPORT_PATH, QgsLayoutExporter.ImageExportSettings()) +if export_result != QgsLayoutExporter.Success: + print(f"Failed to export {LAYER_NAME} map!") +else: + print(f"Map exported to: {EXPORT_PATH}") + +# Cleanup QGIS +qgs.exitQgis() diff --git a/2024_Remote_Sensing/lab3/hist.py b/2024_Remote_Sensing/lab3/hist.py new file mode 100644 index 0000000..c806fc5 --- /dev/null +++ b/2024_Remote_Sensing/lab3/hist.py @@ -0,0 +1,46 @@ +import argparse +import rasterio +import numpy as np +import matplotlib.pyplot as plt + +def plot_histogram(tiff_path, output_path=None): + """Plot a linear histogram for a given TIFF file and save it if an output path is provided.""" + try: + with rasterio.open(tiff_path) as src: + data = src.read(1) # Read the first band + + # Mask nodata values + if src.nodata is not None: + data = np.ma.masked_equal(data, src.nodata) + + # Flatten the array to 1D for histogram plotting + data = data.compressed() if np.ma.isMaskedArray(data) else data.flatten() + + # Plot the histogram + plt.figure(figsize=(10, 6)) + plt.hist(data, bins=256, range=(np.min(data), np.max(data)), edgecolor='black', alpha=0.7) + plt.title(f"Histogram of {tiff_path}") + plt.xlabel("Pixel Value") + plt.ylabel("Frequency") + plt.grid(axis='y', alpha=0.75) + + # Save the plot if an output path is specified + if output_path: + plt.savefig(output_path, dpi=300) + print(f"Histogram saved to {output_path}") + else: + plt.show() + except Exception as e: + print(f"Error processing file {tiff_path}: {e}") + +def main(): + parser = argparse.ArgumentParser(description="Plot a histogram for a TIFF file.") + parser.add_argument("tiff_path", help="Path to the TIFF file.") + parser.add_argument("--output", help="Path to save the histogram plot (optional).", default=None) + + args = parser.parse_args() + + plot_histogram(args.tiff_path, args.output) + +if __name__ == "__main__": + main() diff --git a/2024_Remote_Sensing/lab3/ndvi.py b/2024_Remote_Sensing/lab3/ndvi.py new file mode 100644 index 0000000..cbf15b4 --- /dev/null +++ b/2024_Remote_Sensing/lab3/ndvi.py @@ -0,0 +1,77 @@ +from osgeo import gdal +import numpy as np +import sys +import os +import re + +def tif_to_array(path): + data = gdal.Open(path) + if data is None: + raise FileNotFoundError(f"Cannot open file: {path}") + band = data.GetRasterBand(1) + array = band.ReadAsArray().astype(np.float32) + return data, array + +def handle_nodata(array): + # Normalize the reflectance values to [0, 1] + normalized_array = array / 10000.0 + # Replace out-of-range values with NaN + clean_array = np.where((normalized_array < 0) | (normalized_array > 1), np.nan, normalized_array) + return clean_array + + +def ndvi_calc(NIR,RED): + ndvi_array = (NIR - RED) / (NIR + RED + 1e-10) + return ndvi_array + +def save_tif(array,data,output_path): + driver = gdal.GetDriverByName('GTiff') + rows, cols = array.shape + out_data = driver.Create(output_path, cols, rows, 1, gdal.GDT_Float32) + out_data.SetGeoTransform(data.GetGeoTransform()) + out_data.SetProjection(data.GetProjection()) + out_band = out_data.GetRasterBand(1) + out_band.WriteArray(array) + out_band.SetNoDataValue(-9999) + out_band.FlushCache() + out_data = None + +def get_output_path(output_dir,input_path): + basename = os.path.basename(input_path) + basename_no_ext = os.path.splitext(basename)[0] + basename_no_ext = re.sub(r'_B[1-9]|_B10|_B11', '', basename_no_ext) + path = os.path.join(output_dir, f"{basename_no_ext}_ndvi.TIF") + return path + +def print_min_max(array, label): + print(label," Min:", np.nanmin(array), " Max:",np.nanmax(array)) + +def main(): + if len(sys.argv) != 4: + print("Usage: python3 script.py ") + sys.exit(1) + + red_path = sys.argv[1] + nir_path = sys.argv[2] + output_dir = sys.argv[3] + + red_data, red_array = tif_to_array(red_path) + nir_data, nir_array = tif_to_array(nir_path) + red_array = handle_nodata(red_array) + nir_array = handle_nodata(nir_array) + + ndvi_array = ndvi_calc(nir_array, red_array) + + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + output_path = get_output_path(output_dir,red_path) + + save_tif(ndvi_array, red_data, output_path) + + print_min_max(red_array, "RED") + print_min_max(nir_array, "NIR") + print_min_max(ndvi_array, "NDVI") + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/2024_Remote_Sensing/lab3/ndvi_map.py b/2024_Remote_Sensing/lab3/ndvi_map.py new file mode 100644 index 0000000..e898dc8 --- /dev/null +++ b/2024_Remote_Sensing/lab3/ndvi_map.py @@ -0,0 +1,171 @@ +from qgis.core import * +from PyQt5.QtGui import * +from PyQt5.QtCore import * +from PyQt5 import * +from datetime import datetime +import os + +QGIS_PREFIX_PATH = "/usr/share/qgis" +LAYER_PATH = "/home/huaqo/dev/courses/2024_Remote_Sensing/lab3/ndvi/B04_ndvi.TIF" +STYLE_PATH = "/home/huaqo/dev/courses/2024_Remote_Sensing/styles/ndvi.qml" +EXPORT_PATH = "/home/huaqo/dev/courses/2024_Remote_Sensing/lab3/maps/ndvi.png" +NORTH_ARROW_PATH = "/home/huaqo/dev/courses/2024_Remote_Sensing/styles/north_arrow.svg" + +map_pos_x = 20 +map_pos_y = 30 +map_size_x = 200 +map_size_y = 150 + +grid_interval_x = 100000 +grid_interval_y = 100000 + +LAYER_NAME = "NDVI" +legend_pos_x = 245 +legend_pos_y = 30 + +scalebar_pos_x = 20 +scalebar_pos_y = 190 + +north_arrow_pos_x = 30 +north_arrow_pos_y = 40 +north_arrow_size_x = 15 +north_arrow_size_y = 15 + +title_text = "NDVI - Lab 3" +title_fontsize = 16 +title_height = 10 + +author_text = "Author: Joaquin Gottlebe" +author_pos_y = 100 + +date_text = f"Date: {datetime.now().strftime('%d.%m.%Y')}" +date_pos_y = 110 + +crs_text = "CRS: EPSG 32633" +crs_pos_y = 120 + +fontfamily = "Arial" +info_font_size = 12 +info_pos_x = legend_pos_x + +# Initialize QGIS Application +QgsApplication.setPrefixPath(QGIS_PREFIX_PATH, True) +qgs = QgsApplication([], False) +qgs.initQgis() + +# Cleanup +project = QgsProject.instance() +layer_ids = list(project.mapLayers().keys()) +for layer_id in layer_ids: + project.removeMapLayer(layer_id) + +# Load Layer +raster_layer = QgsRasterLayer(LAYER_PATH, LAYER_NAME) +if not raster_layer.isValid(): + raise Exception(f"{LAYER_NAME} layer failed to load!") + +QgsProject.instance().addMapLayer(raster_layer) +if not raster_layer.loadNamedStyle(STYLE_PATH): + raise Exception(f"Style {STYLE_PATH} failed to load") + +# Create Layout +layout = QgsPrintLayout(project) +layout.initializeDefaults() + +# Map +map_item = QgsLayoutItemMap(layout) +map_item.attemptMove(QgsLayoutPoint(map_pos_x,map_pos_y, QgsUnitTypes.LayoutMillimeters)) +map_item.attemptResize(QgsLayoutSize(map_size_x, map_size_y, QgsUnitTypes.LayoutMillimeters)) +map_item.zoomToExtent(raster_layer.extent()) +layout.addLayoutItem(map_item) + +# Add Coordinate Grid +grid = map_item.grid() +grid.setEnabled(True) +grid.setIntervalX(grid_interval_x) +grid.setIntervalY(grid_interval_y) +grid.setAnnotationEnabled(True) +grid.setAnnotationPrecision(0) +grid.setFrameStyle(QgsLayoutItemMapGrid.Zebra) + +# Label +title_label = QgsLayoutItemLabel(layout) +title_label.setText(title_text) +title_label.setFont(QFont(fontfamily, title_fontsize)) +title_label.setHAlign(Qt.AlignCenter) +title_label.adjustSizeToText() +layout_width = layout.pageCollection().page(0).pageSize().width() +title_width = title_label.rectWithFrame().width() +center_x = (layout_width - title_width) / 2 +title_label.attemptMove(QgsLayoutPoint(center_x, title_height, QgsUnitTypes.LayoutMillimeters)) +layout.addLayoutItem(title_label) + +# Legend +legend_item = QgsLayoutItemLegend(layout) +legend_item.setLinkedMap(map_item) +legend_item.setTitle("Legend") +legend_item.setFrameEnabled(True) +legend_item.attemptMove(QgsLayoutPoint(legend_pos_x,legend_pos_y, QgsUnitTypes.LayoutMillimeters)) +layout.addLayoutItem(legend_item) + +# Scalebar +scalebar_item = QgsLayoutItemScaleBar(layout) +scalebar_item.setStyle('Single Box') +scalebar_item.setLinkedMap(map_item) +scalebar_item.setUnitLabel('m') +scalebar_item.setNumberOfSegments(4) +scalebar_item.setNumberOfSegmentsLeft(0) +scalebar_item.setUnitsPerSegment(50000) +scalebar_item.applyDefaultSize() +scalebar_item.attemptMove(QgsLayoutPoint(scalebar_pos_x,scalebar_pos_y, QgsUnitTypes.LayoutMillimeters)) +layout.addLayoutItem(scalebar_item) + +# North Arrow +north_arrow_item = QgsLayoutItemPicture(layout) +north_arrow_item.setPicturePath(NORTH_ARROW_PATH) +if not north_arrow_item.picturePath(): + print("North arrow image failed to load!") +north_arrow_item.setReferencePoint(QgsLayoutItemPicture.UpperLeft) +north_arrow_item.attemptMove(QgsLayoutPoint(north_arrow_pos_x,north_arrow_pos_y,QgsUnitTypes.LayoutMillimeters)) +north_arrow_item.attemptResize(QgsLayoutSize(north_arrow_size_x,north_arrow_size_y,QgsUnitTypes.LayoutMillimeters)) +layout.addLayoutItem(north_arrow_item) + +# Author Info +info_label = QgsLayoutItemLabel(layout) +info_label.setText(author_text) +info_label.setFont(QFont(fontfamily, info_font_size)) +info_label.setHAlign(Qt.AlignLeft) +info_label.attemptMove(QgsLayoutPoint(info_pos_x, author_pos_y, QgsUnitTypes.LayoutMillimeters)) +info_label.adjustSizeToText() +layout.addLayoutItem(info_label) + +# Date Info +info_label = QgsLayoutItemLabel(layout) +info_label.setText(date_text) +info_label.setFont(QFont(fontfamily, info_font_size)) +info_label.setHAlign(Qt.AlignLeft) +info_label.attemptMove(QgsLayoutPoint(info_pos_x, date_pos_y, QgsUnitTypes.LayoutMillimeters)) +info_label.adjustSizeToText() +layout.addLayoutItem(info_label) + +# CRS Info +info_label = QgsLayoutItemLabel(layout) +info_label.setText(crs_text) +info_label.setFont(QFont(fontfamily, info_font_size)) +info_label.setHAlign(Qt.AlignLeft) +info_label.attemptMove(QgsLayoutPoint(info_pos_x, crs_pos_y, QgsUnitTypes.LayoutMillimeters)) +info_label.adjustSizeToText() +layout.addLayoutItem(info_label) + +# Export Map +exporter = QgsLayoutExporter(layout) +if os.path.exists(EXPORT_PATH): + os.remove(EXPORT_PATH) +export_result = exporter.exportToImage(EXPORT_PATH, QgsLayoutExporter.ImageExportSettings()) +if export_result != QgsLayoutExporter.Success: + print(f"Failed to export {LAYER_NAME} map!") +else: + print(f"Map exported to: {EXPORT_PATH}") + +# Cleanup QGIS +qgs.exitQgis() diff --git a/2024_Remote_Sensing/lab3/ndvi_map2.py b/2024_Remote_Sensing/lab3/ndvi_map2.py new file mode 100644 index 0000000..f030c22 --- /dev/null +++ b/2024_Remote_Sensing/lab3/ndvi_map2.py @@ -0,0 +1,171 @@ +from qgis.core import * +from PyQt5.QtGui import * +from PyQt5.QtCore import * +from PyQt5 import * +from datetime import datetime +import os + +QGIS_PREFIX_PATH = "/usr/share/qgis" +LAYER_PATH = "/home/huaqo/dev/courses/2024_Remote_Sensing/lab3/ndvi2/B04_ndvi.TIF" +STYLE_PATH = "/home/huaqo/dev/courses/2024_Remote_Sensing/styles/ndvi2.qml" +EXPORT_PATH = "/home/huaqo/dev/courses/2024_Remote_Sensing/lab3/maps/ndvi2.png" +NORTH_ARROW_PATH = "/home/huaqo/dev/courses/2024_Remote_Sensing/styles/north_arrow.svg" + +map_pos_x = 20 +map_pos_y = 30 +map_size_x = 200 +map_size_y = 150 + +grid_interval_x = 100000 +grid_interval_y = 100000 + +LAYER_NAME = "NDVI" +legend_pos_x = 245 +legend_pos_y = 30 + +scalebar_pos_x = 20 +scalebar_pos_y = 190 + +north_arrow_pos_x = 30 +north_arrow_pos_y = 40 +north_arrow_size_x = 15 +north_arrow_size_y = 15 + +title_text = "Chacao Channel 18.11.24 NDVI - Lab 3" +title_fontsize = 16 +title_height = 10 + +author_text = "Author: Joaquin Gottlebe" +author_pos_y = 100 + +date_text = f"Date: {datetime.now().strftime('%d.%m.%Y')}" +date_pos_y = 110 + +crs_text = "CRS: EPSG 32633" +crs_pos_y = 120 + +fontfamily = "Arial" +info_font_size = 12 +info_pos_x = legend_pos_x + +# Initialize QGIS Application +QgsApplication.setPrefixPath(QGIS_PREFIX_PATH, True) +qgs = QgsApplication([], False) +qgs.initQgis() + +# Cleanup +project = QgsProject.instance() +layer_ids = list(project.mapLayers().keys()) +for layer_id in layer_ids: + project.removeMapLayer(layer_id) + +# Load Layer +raster_layer = QgsRasterLayer(LAYER_PATH, LAYER_NAME) +if not raster_layer.isValid(): + raise Exception(f"{LAYER_NAME} layer failed to load!") + +QgsProject.instance().addMapLayer(raster_layer) +if not raster_layer.loadNamedStyle(STYLE_PATH): + raise Exception(f"Style {STYLE_PATH} failed to load") + +# Create Layout +layout = QgsPrintLayout(project) +layout.initializeDefaults() + +# Map +map_item = QgsLayoutItemMap(layout) +map_item.attemptMove(QgsLayoutPoint(map_pos_x,map_pos_y, QgsUnitTypes.LayoutMillimeters)) +map_item.attemptResize(QgsLayoutSize(map_size_x, map_size_y, QgsUnitTypes.LayoutMillimeters)) +map_item.zoomToExtent(raster_layer.extent()) +layout.addLayoutItem(map_item) + +# Add Coordinate Grid +grid = map_item.grid() +grid.setEnabled(True) +grid.setIntervalX(grid_interval_x) +grid.setIntervalY(grid_interval_y) +grid.setAnnotationEnabled(True) +grid.setAnnotationPrecision(0) +grid.setFrameStyle(QgsLayoutItemMapGrid.Zebra) + +# Label +title_label = QgsLayoutItemLabel(layout) +title_label.setText(title_text) +title_label.setFont(QFont(fontfamily, title_fontsize)) +title_label.setHAlign(Qt.AlignCenter) +title_label.adjustSizeToText() +layout_width = layout.pageCollection().page(0).pageSize().width() +title_width = title_label.rectWithFrame().width() +center_x = (layout_width - title_width) / 2 +title_label.attemptMove(QgsLayoutPoint(center_x, title_height, QgsUnitTypes.LayoutMillimeters)) +layout.addLayoutItem(title_label) + +# Legend +legend_item = QgsLayoutItemLegend(layout) +legend_item.setLinkedMap(map_item) +legend_item.setTitle("Legend") +legend_item.setFrameEnabled(True) +legend_item.attemptMove(QgsLayoutPoint(legend_pos_x,legend_pos_y, QgsUnitTypes.LayoutMillimeters)) +layout.addLayoutItem(legend_item) + +# Scalebar +scalebar_item = QgsLayoutItemScaleBar(layout) +scalebar_item.setStyle('Single Box') +scalebar_item.setLinkedMap(map_item) +scalebar_item.setUnitLabel('m') +scalebar_item.setNumberOfSegments(4) +scalebar_item.setNumberOfSegmentsLeft(0) +scalebar_item.setUnitsPerSegment(50000) +scalebar_item.applyDefaultSize() +scalebar_item.attemptMove(QgsLayoutPoint(scalebar_pos_x,scalebar_pos_y, QgsUnitTypes.LayoutMillimeters)) +layout.addLayoutItem(scalebar_item) + +# North Arrow +north_arrow_item = QgsLayoutItemPicture(layout) +north_arrow_item.setPicturePath(NORTH_ARROW_PATH) +if not north_arrow_item.picturePath(): + print("North arrow image failed to load!") +north_arrow_item.setReferencePoint(QgsLayoutItemPicture.UpperLeft) +north_arrow_item.attemptMove(QgsLayoutPoint(north_arrow_pos_x,north_arrow_pos_y,QgsUnitTypes.LayoutMillimeters)) +north_arrow_item.attemptResize(QgsLayoutSize(north_arrow_size_x,north_arrow_size_y,QgsUnitTypes.LayoutMillimeters)) +layout.addLayoutItem(north_arrow_item) + +# Author Info +info_label = QgsLayoutItemLabel(layout) +info_label.setText(author_text) +info_label.setFont(QFont(fontfamily, info_font_size)) +info_label.setHAlign(Qt.AlignLeft) +info_label.attemptMove(QgsLayoutPoint(info_pos_x, author_pos_y, QgsUnitTypes.LayoutMillimeters)) +info_label.adjustSizeToText() +layout.addLayoutItem(info_label) + +# Date Info +info_label = QgsLayoutItemLabel(layout) +info_label.setText(date_text) +info_label.setFont(QFont(fontfamily, info_font_size)) +info_label.setHAlign(Qt.AlignLeft) +info_label.attemptMove(QgsLayoutPoint(info_pos_x, date_pos_y, QgsUnitTypes.LayoutMillimeters)) +info_label.adjustSizeToText() +layout.addLayoutItem(info_label) + +# CRS Info +info_label = QgsLayoutItemLabel(layout) +info_label.setText(crs_text) +info_label.setFont(QFont(fontfamily, info_font_size)) +info_label.setHAlign(Qt.AlignLeft) +info_label.attemptMove(QgsLayoutPoint(info_pos_x, crs_pos_y, QgsUnitTypes.LayoutMillimeters)) +info_label.adjustSizeToText() +layout.addLayoutItem(info_label) + +# Export Map +exporter = QgsLayoutExporter(layout) +if os.path.exists(EXPORT_PATH): + os.remove(EXPORT_PATH) +export_result = exporter.exportToImage(EXPORT_PATH, QgsLayoutExporter.ImageExportSettings()) +if export_result != QgsLayoutExporter.Success: + print(f"Failed to export {LAYER_NAME} map!") +else: + print(f"Map exported to: {EXPORT_PATH}") + +# Cleanup QGIS +qgs.exitQgis() diff --git a/2024_Remote_Sensing/lab3/ndwi.py b/2024_Remote_Sensing/lab3/ndwi.py new file mode 100644 index 0000000..1113d87 --- /dev/null +++ b/2024_Remote_Sensing/lab3/ndwi.py @@ -0,0 +1,84 @@ +from osgeo import gdal +import numpy as np +import sys +import os +import re + + +def tif_to_array(path): + data = gdal.Open(path) + if data is None: + raise FileNotFoundError(f"Cannot open file: {path}") + band = data.GetRasterBand(1) + array = band.ReadAsArray().astype(np.float32) + return data, array + + +def handle_nodata(array): + # Normalize the reflectance values to [0, 1] + normalized_array = array / 10000.0 + # Replace out-of-range values with NaN + clean_array = np.where((normalized_array < 0) | (normalized_array > 1), np.nan, normalized_array) + return clean_array + + +def ndwi_calc(GREEN, NIR): + ndwi_array = (GREEN - NIR) / (GREEN + NIR + 1e-10) + return ndwi_array + + +def save_tif(array, data, output_path): + driver = gdal.GetDriverByName('GTiff') + rows, cols = array.shape + out_data = driver.Create(output_path, cols, rows, 1, gdal.GDT_Float32) + out_data.SetGeoTransform(data.GetGeoTransform()) + out_data.SetProjection(data.GetProjection()) + out_band = out_data.GetRasterBand(1) + out_band.WriteArray(array) + out_band.SetNoDataValue(-9999) + out_band.FlushCache() + out_data = None + + +def get_output_path(output_dir, input_path): + basename = os.path.basename(input_path) + basename_no_ext = os.path.splitext(basename)[0] + basename_no_ext = re.sub(r'_B[1-9]|_B10|_B11', '', basename_no_ext) + path = os.path.join(output_dir, f"{basename_no_ext}_ndwi.TIF") + return path + + +def print_min_max(array, label): + print(label, " Min:", np.nanmin(array), " Max:", np.nanmax(array)) + + +def main(): + if len(sys.argv) != 4: + print("Usage: python3 script.py ") + sys.exit(1) + + green_path = sys.argv[1] + nir_path = sys.argv[2] + output_dir = sys.argv[3] + + green_data, green_array = tif_to_array(green_path) + nir_data, nir_array = tif_to_array(nir_path) + green_array = handle_nodata(green_array) + nir_array = handle_nodata(nir_array) + + ndwi_array = ndwi_calc(green_array, nir_array) + + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + output_path = get_output_path(output_dir, green_path) + + save_tif(ndwi_array, green_data, output_path) + + print_min_max(nir_array, "NIR") + print_min_max(green_array, "GREEN") + print_min_max(ndwi_array, "NDWI") + + +if __name__ == "__main__": + main() diff --git a/2024_Remote_Sensing/lab3/ndwi_map.py b/2024_Remote_Sensing/lab3/ndwi_map.py new file mode 100644 index 0000000..50db323 --- /dev/null +++ b/2024_Remote_Sensing/lab3/ndwi_map.py @@ -0,0 +1,171 @@ +from qgis.core import * +from PyQt5.QtGui import * +from PyQt5.QtCore import * +from PyQt5 import * +from datetime import datetime +import os + +QGIS_PREFIX_PATH = "/usr/share/qgis" +LAYER_PATH = "/home/huaqo/dev/courses/2024_Remote_Sensing/lab3/ndwi/B03_ndwi.TIF" +STYLE_PATH = "/home/huaqo/dev/courses/2024_Remote_Sensing/styles/ndwi.qml" +EXPORT_PATH = "/home/huaqo/dev/courses/2024_Remote_Sensing/lab3/maps/ndwi.png" +NORTH_ARROW_PATH = "/home/huaqo/dev/courses/2024_Remote_Sensing/styles/north_arrow.svg" + +map_pos_x = 20 +map_pos_y = 30 +map_size_x = 200 +map_size_y = 150 + +grid_interval_x = 100000 +grid_interval_y = 100000 + +LAYER_NAME = "NDWI" +legend_pos_x = 245 +legend_pos_y = 30 + +scalebar_pos_x = 20 +scalebar_pos_y = 190 + +north_arrow_pos_x = 30 +north_arrow_pos_y = 40 +north_arrow_size_x = 15 +north_arrow_size_y = 15 + +title_text = "NDWI - Lab 3" +title_fontsize = 16 +title_height = 10 + +author_text = "Author: Joaquin Gottlebe" +author_pos_y = 100 + +date_text = f"Date: {datetime.now().strftime('%d.%m.%Y')}" +date_pos_y = 110 + +crs_text = "CRS: EPSG 32633" +crs_pos_y = 120 + +info_font_size = 12 +info_pos_x = legend_pos_x +fontfamily = "Arial" + +# Initialize QGIS Application +QgsApplication.setPrefixPath(QGIS_PREFIX_PATH, True) +qgs = QgsApplication([], False) +qgs.initQgis() + +# Cleanup +project = QgsProject.instance() +layer_ids = list(project.mapLayers().keys()) +for layer_id in layer_ids: + project.removeMapLayer(layer_id) + +# Load Layer +raster_layer = QgsRasterLayer(LAYER_PATH, LAYER_NAME) +if not raster_layer.isValid(): + raise Exception(f"{LAYER_NAME} layer failed to load!") + +QgsProject.instance().addMapLayer(raster_layer) +if not raster_layer.loadNamedStyle(STYLE_PATH): + raise Exception(f"Style {STYLE_PATH} failed to load") + +# Create Layout +layout = QgsPrintLayout(project) +layout.initializeDefaults() + +# Map +map_item = QgsLayoutItemMap(layout) +map_item.attemptMove(QgsLayoutPoint(map_pos_x,map_pos_y, QgsUnitTypes.LayoutMillimeters)) +map_item.attemptResize(QgsLayoutSize(map_size_x, map_size_y, QgsUnitTypes.LayoutMillimeters)) +map_item.zoomToExtent(raster_layer.extent()) +layout.addLayoutItem(map_item) + +# Add Coordinate Grid +grid = map_item.grid() +grid.setEnabled(True) +grid.setIntervalX(grid_interval_x) +grid.setIntervalY(grid_interval_y) +grid.setAnnotationEnabled(True) +grid.setAnnotationPrecision(0) +grid.setFrameStyle(QgsLayoutItemMapGrid.Zebra) + +# Label +title_label = QgsLayoutItemLabel(layout) +title_label.setText(title_text) +title_label.setFont(QFont(fontfamily, title_fontsize)) +title_label.setHAlign(Qt.AlignCenter) +title_label.adjustSizeToText() +layout_width = layout.pageCollection().page(0).pageSize().width() +title_width = title_label.rectWithFrame().width() +center_x = (layout_width - title_width) / 2 +title_label.attemptMove(QgsLayoutPoint(center_x, title_height, QgsUnitTypes.LayoutMillimeters)) +layout.addLayoutItem(title_label) + +# Legend +legend_item = QgsLayoutItemLegend(layout) +legend_item.setLinkedMap(map_item) +legend_item.setTitle("Legend") +legend_item.setFrameEnabled(True) +legend_item.attemptMove(QgsLayoutPoint(legend_pos_x,legend_pos_y, QgsUnitTypes.LayoutMillimeters)) +layout.addLayoutItem(legend_item) + +# Scalebar +scalebar_item = QgsLayoutItemScaleBar(layout) +scalebar_item.setStyle('Single Box') +scalebar_item.setLinkedMap(map_item) +scalebar_item.setUnitLabel('m') +scalebar_item.setNumberOfSegments(4) +scalebar_item.setNumberOfSegmentsLeft(0) +scalebar_item.setUnitsPerSegment(50000) +scalebar_item.applyDefaultSize() +scalebar_item.attemptMove(QgsLayoutPoint(scalebar_pos_x,scalebar_pos_y, QgsUnitTypes.LayoutMillimeters)) +layout.addLayoutItem(scalebar_item) + +# North Arrow +north_arrow_item = QgsLayoutItemPicture(layout) +north_arrow_item.setPicturePath(NORTH_ARROW_PATH) +if not north_arrow_item.picturePath(): + print("North arrow image failed to load!") +north_arrow_item.setReferencePoint(QgsLayoutItemPicture.UpperLeft) +north_arrow_item.attemptMove(QgsLayoutPoint(north_arrow_pos_x,north_arrow_pos_y,QgsUnitTypes.LayoutMillimeters)) +north_arrow_item.attemptResize(QgsLayoutSize(north_arrow_size_x,north_arrow_size_y,QgsUnitTypes.LayoutMillimeters)) +layout.addLayoutItem(north_arrow_item) + +# Author Info +info_label = QgsLayoutItemLabel(layout) +info_label.setText(author_text) +info_label.setFont(QFont(fontfamily, info_font_size)) +info_label.setHAlign(Qt.AlignLeft) +info_label.attemptMove(QgsLayoutPoint(info_pos_x, author_pos_y, QgsUnitTypes.LayoutMillimeters)) +info_label.adjustSizeToText() +layout.addLayoutItem(info_label) + +# Date Info +info_label = QgsLayoutItemLabel(layout) +info_label.setText(date_text) +info_label.setFont(QFont(fontfamily, info_font_size)) +info_label.setHAlign(Qt.AlignLeft) +info_label.attemptMove(QgsLayoutPoint(info_pos_x, date_pos_y, QgsUnitTypes.LayoutMillimeters)) +info_label.adjustSizeToText() +layout.addLayoutItem(info_label) + +# CRS Info +info_label = QgsLayoutItemLabel(layout) +info_label.setText(crs_text) +info_label.setFont(QFont(fontfamily, info_font_size)) +info_label.setHAlign(Qt.AlignLeft) +info_label.attemptMove(QgsLayoutPoint(info_pos_x, crs_pos_y, QgsUnitTypes.LayoutMillimeters)) +info_label.adjustSizeToText() +layout.addLayoutItem(info_label) + +# Export Map +exporter = QgsLayoutExporter(layout) +if os.path.exists(EXPORT_PATH): + os.remove(EXPORT_PATH) +export_result = exporter.exportToImage(EXPORT_PATH, QgsLayoutExporter.ImageExportSettings()) +if export_result != QgsLayoutExporter.Success: + print(f"Failed to export {LAYER_NAME} map!") +else: + print(f"Map exported to: {EXPORT_PATH}") + +# Cleanup QGIS +qgs.exitQgis() diff --git a/2024_Remote_Sensing/lab3/ndwi_map2.py b/2024_Remote_Sensing/lab3/ndwi_map2.py new file mode 100644 index 0000000..56f8f17 --- /dev/null +++ b/2024_Remote_Sensing/lab3/ndwi_map2.py @@ -0,0 +1,171 @@ +from qgis.core import * +from PyQt5.QtGui import * +from PyQt5.QtCore import * +from PyQt5 import * +from datetime import datetime +import os + +QGIS_PREFIX_PATH = "/usr/share/qgis" +LAYER_PATH = "/home/huaqo/dev/courses/2024_Remote_Sensing/lab3/ndwi2/B03_ndwi.TIF" +STYLE_PATH = "/home/huaqo/dev/courses/2024_Remote_Sensing/styles/ndwi2.qml" +EXPORT_PATH = "/home/huaqo/dev/courses/2024_Remote_Sensing/lab3/maps/ndwi2.png" +NORTH_ARROW_PATH = "/home/huaqo/dev/courses/2024_Remote_Sensing/styles/north_arrow.svg" + +map_pos_x = 20 +map_pos_y = 30 +map_size_x = 200 +map_size_y = 150 + +grid_interval_x = 100000 +grid_interval_y = 100000 + +LAYER_NAME = "NDWI" +legend_pos_x = 245 +legend_pos_y = 30 + +scalebar_pos_x = 20 +scalebar_pos_y = 190 + +north_arrow_pos_x = 30 +north_arrow_pos_y = 40 +north_arrow_size_x = 15 +north_arrow_size_y = 15 + +title_text = "Chacao Channel 18.11.24 NDWI - Lab 3" +title_fontsize = 16 +title_height = 10 + +author_text = "Author: Joaquin Gottlebe" +author_pos_y = 100 + +date_text = f"Date: {datetime.now().strftime('%d.%m.%Y')}" +date_pos_y = 110 + +crs_text = "CRS: EPSG 32633" +crs_pos_y = 120 + +info_font_size = 12 +info_pos_x = legend_pos_x +fontfamily = "Arial" + +# Initialize QGIS Application +QgsApplication.setPrefixPath(QGIS_PREFIX_PATH, True) +qgs = QgsApplication([], False) +qgs.initQgis() + +# Cleanup +project = QgsProject.instance() +layer_ids = list(project.mapLayers().keys()) +for layer_id in layer_ids: + project.removeMapLayer(layer_id) + +# Load Layer +raster_layer = QgsRasterLayer(LAYER_PATH, LAYER_NAME) +if not raster_layer.isValid(): + raise Exception(f"{LAYER_NAME} layer failed to load!") + +QgsProject.instance().addMapLayer(raster_layer) +if not raster_layer.loadNamedStyle(STYLE_PATH): + raise Exception(f"Style {STYLE_PATH} failed to load") + +# Create Layout +layout = QgsPrintLayout(project) +layout.initializeDefaults() + +# Map +map_item = QgsLayoutItemMap(layout) +map_item.attemptMove(QgsLayoutPoint(map_pos_x,map_pos_y, QgsUnitTypes.LayoutMillimeters)) +map_item.attemptResize(QgsLayoutSize(map_size_x, map_size_y, QgsUnitTypes.LayoutMillimeters)) +map_item.zoomToExtent(raster_layer.extent()) +layout.addLayoutItem(map_item) + +# Add Coordinate Grid +grid = map_item.grid() +grid.setEnabled(True) +grid.setIntervalX(grid_interval_x) +grid.setIntervalY(grid_interval_y) +grid.setAnnotationEnabled(True) +grid.setAnnotationPrecision(0) +grid.setFrameStyle(QgsLayoutItemMapGrid.Zebra) + +# Label +title_label = QgsLayoutItemLabel(layout) +title_label.setText(title_text) +title_label.setFont(QFont(fontfamily, title_fontsize)) +title_label.setHAlign(Qt.AlignCenter) +title_label.adjustSizeToText() +layout_width = layout.pageCollection().page(0).pageSize().width() +title_width = title_label.rectWithFrame().width() +center_x = (layout_width - title_width) / 2 +title_label.attemptMove(QgsLayoutPoint(center_x, title_height, QgsUnitTypes.LayoutMillimeters)) +layout.addLayoutItem(title_label) + +# Legend +legend_item = QgsLayoutItemLegend(layout) +legend_item.setLinkedMap(map_item) +legend_item.setTitle("Legend") +legend_item.setFrameEnabled(True) +legend_item.attemptMove(QgsLayoutPoint(legend_pos_x,legend_pos_y, QgsUnitTypes.LayoutMillimeters)) +layout.addLayoutItem(legend_item) + +# Scalebar +scalebar_item = QgsLayoutItemScaleBar(layout) +scalebar_item.setStyle('Single Box') +scalebar_item.setLinkedMap(map_item) +scalebar_item.setUnitLabel('m') +scalebar_item.setNumberOfSegments(4) +scalebar_item.setNumberOfSegmentsLeft(0) +scalebar_item.setUnitsPerSegment(50000) +scalebar_item.applyDefaultSize() +scalebar_item.attemptMove(QgsLayoutPoint(scalebar_pos_x,scalebar_pos_y, QgsUnitTypes.LayoutMillimeters)) +layout.addLayoutItem(scalebar_item) + +# North Arrow +north_arrow_item = QgsLayoutItemPicture(layout) +north_arrow_item.setPicturePath(NORTH_ARROW_PATH) +if not north_arrow_item.picturePath(): + print("North arrow image failed to load!") +north_arrow_item.setReferencePoint(QgsLayoutItemPicture.UpperLeft) +north_arrow_item.attemptMove(QgsLayoutPoint(north_arrow_pos_x,north_arrow_pos_y,QgsUnitTypes.LayoutMillimeters)) +north_arrow_item.attemptResize(QgsLayoutSize(north_arrow_size_x,north_arrow_size_y,QgsUnitTypes.LayoutMillimeters)) +layout.addLayoutItem(north_arrow_item) + +# Author Info +info_label = QgsLayoutItemLabel(layout) +info_label.setText(author_text) +info_label.setFont(QFont(fontfamily, info_font_size)) +info_label.setHAlign(Qt.AlignLeft) +info_label.attemptMove(QgsLayoutPoint(info_pos_x, author_pos_y, QgsUnitTypes.LayoutMillimeters)) +info_label.adjustSizeToText() +layout.addLayoutItem(info_label) + +# Date Info +info_label = QgsLayoutItemLabel(layout) +info_label.setText(date_text) +info_label.setFont(QFont(fontfamily, info_font_size)) +info_label.setHAlign(Qt.AlignLeft) +info_label.attemptMove(QgsLayoutPoint(info_pos_x, date_pos_y, QgsUnitTypes.LayoutMillimeters)) +info_label.adjustSizeToText() +layout.addLayoutItem(info_label) + +# CRS Info +info_label = QgsLayoutItemLabel(layout) +info_label.setText(crs_text) +info_label.setFont(QFont(fontfamily, info_font_size)) +info_label.setHAlign(Qt.AlignLeft) +info_label.attemptMove(QgsLayoutPoint(info_pos_x, crs_pos_y, QgsUnitTypes.LayoutMillimeters)) +info_label.adjustSizeToText() +layout.addLayoutItem(info_label) + +# Export Map +exporter = QgsLayoutExporter(layout) +if os.path.exists(EXPORT_PATH): + os.remove(EXPORT_PATH) +export_result = exporter.exportToImage(EXPORT_PATH, QgsLayoutExporter.ImageExportSettings()) +if export_result != QgsLayoutExporter.Success: + print(f"Failed to export {LAYER_NAME} map!") +else: + print(f"Map exported to: {EXPORT_PATH}") + +# Cleanup QGIS +qgs.exitQgis() diff --git a/2024_Remote_Sensing/lab4/Makefile b/2024_Remote_Sensing/lab4/Makefile new file mode 100644 index 0000000..e69de29 diff --git a/2024_Remote_Sensing/styles/diff.qml b/2024_Remote_Sensing/styles/diff.qml new file mode 100644 index 0000000..82b159b --- /dev/null +++ b/2024_Remote_Sensing/styles/diff.qml @@ -0,0 +1,223 @@ + + + + 1 + 1 + 1 + 0 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + None + WholeRaster + Estimated + 0.02 + 0.98 + 2 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + resamplingFilter + + 0 + diff --git a/2024_Remote_Sensing/styles/false_color.qml b/2024_Remote_Sensing/styles/false_color.qml new file mode 100644 index 0000000..2489b52 --- /dev/null +++ b/2024_Remote_Sensing/styles/false_color.qml @@ -0,0 +1,158 @@ + + + + 1 + 1 + 1 + 0 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + CumulativeCut + WholeRaster + Estimated + 0.02 + 0.98 + 2 + + + -0.119304 + 0.514984 + StretchToMinimumMaximum + + + -0.119304 + 0.168884 + StretchToMinimumMaximum + + + -0.119304 + 0.149555 + StretchToMinimumMaximum + + + + + + resamplingFilter + + 0 + diff --git a/2024_Remote_Sensing/styles/ndvi.qml b/2024_Remote_Sensing/styles/ndvi.qml new file mode 100644 index 0000000..0212059 --- /dev/null +++ b/2024_Remote_Sensing/styles/ndvi.qml @@ -0,0 +1,176 @@ + + + + 1 + 1 + 1 + 0 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + None + WholeRaster + Estimated + 0.02 + 0.98 + 2 + + + + + + + + + + + + + + + + + + + + + + + resamplingFilter + + 0 + diff --git a/2024_Remote_Sensing/styles/ndvi2.qml b/2024_Remote_Sensing/styles/ndvi2.qml new file mode 100644 index 0000000..f51ff65 --- /dev/null +++ b/2024_Remote_Sensing/styles/ndvi2.qml @@ -0,0 +1,176 @@ + + + + 1 + 1 + 1 + 0 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + None + WholeRaster + Estimated + 0.02 + 0.98 + 2 + + + + + + + + + + + + + + + + + + + + + + + resamplingFilter + + 0 + diff --git a/2024_Remote_Sensing/styles/ndwi.qml b/2024_Remote_Sensing/styles/ndwi.qml new file mode 100644 index 0000000..82e353e --- /dev/null +++ b/2024_Remote_Sensing/styles/ndwi.qml @@ -0,0 +1,172 @@ + + + + 1 + 1 + 1 + 0 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + None + WholeRaster + Estimated + 0.02 + 0.98 + 2 + + + + + + + + + + + + + + + + + + + + resamplingFilter + + 0 + diff --git a/2024_Remote_Sensing/styles/ndwi2.qml b/2024_Remote_Sensing/styles/ndwi2.qml new file mode 100644 index 0000000..2474513 --- /dev/null +++ b/2024_Remote_Sensing/styles/ndwi2.qml @@ -0,0 +1,180 @@ + + + + 1 + 1 + 1 + 0 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + None + WholeRaster + Estimated + 0.02 + 0.98 + 2 + + + + + + + + + + + + + + + + + + + + + + + + + + + resamplingFilter + + 0 + diff --git a/2024_Remote_Sensing/styles/north_arrow.svg b/2024_Remote_Sensing/styles/north_arrow.svg new file mode 100644 index 0000000..218a3e4 --- /dev/null +++ b/2024_Remote_Sensing/styles/north_arrow.svg @@ -0,0 +1,5 @@ + + + + + \ No newline at end of file diff --git a/2024_Remote_Sensing/styles/surface_temperature.qml b/2024_Remote_Sensing/styles/surface_temperature.qml new file mode 100644 index 0000000..e3ead3e --- /dev/null +++ b/2024_Remote_Sensing/styles/surface_temperature.qml @@ -0,0 +1,176 @@ + + + + 1 + 1 + 1 + 0 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + None + WholeRaster + Estimated + 0.02 + 0.98 + 2 + + + + + + + + + + + + + + + + + + + + + + + resamplingFilter + + 0 + diff --git a/2024_Remote_Sensing/styles/true_color.qml b/2024_Remote_Sensing/styles/true_color.qml new file mode 100644 index 0000000..3b74b31 --- /dev/null +++ b/2024_Remote_Sensing/styles/true_color.qml @@ -0,0 +1,158 @@ + + + + 1 + 1 + 1 + 0 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + CumulativeCut + WholeRaster + Estimated + 0.02 + 0.98 + 2 + + + -0.119304 + 0.168884 + StretchToMinimumMaximum + + + -0.119304 + 0.149555 + StretchToMinimumMaximum + + + -0.119304 + 0.143137 + StretchToMinimumMaximum + + + + + + resamplingFilter + + 0 +