1 Στόχοι της άσκησης

  • Φορτωμά δεδομένων vector απο gpkg βάση δεδομένων
  • Μαθηματικές πράξεις σε υποσυνολα vector δεδομένων
  • Απλή οπτικοποίηση χωρικών δεδομένων
  • Καταμέτρηση (Απόλητη και Σχετική συχνότητα) κελιών Raster ανα κατηγορία τιμών
  • Υπολογισμός εμβαδού ανα κατηγορία τιμών σε Raster
  • Επαναδειγματοληψία (resample) Raster δεδομένων
  • Λογική συνθεση δυο Raster δεδομένων
  • Υπολογισμό ΜΟ τιμών Raster ανα πολύγωνο Vector
  • Οπτικοποιηση αποτελεσμάτων με απλούς θεματικούς χάρτες

2 Πρόβλημα

Μια απελευθέρωση ρύπανσης από ένα ερευνητικό εργοστάσιο στο νησί της Λέσβου έχει οδηγήσει σε ένα λοφίο ρύπανσης (pollution plume) στο νησί. Οι επιστήμονες έχουν ένα μοντέλο για την πυκνότητα της ρύπανσης και θέλουν να μάθουν ποιες πόλεις του νησιού θα εκτεθούν σε επίπεδο ρύπανσης πάνω από 8000 ng (νανογραμμάρια) ανά τετραγωνικό μέτρο.

3 Ερωτήματα

  1. Πόσα άτομα επηρεάζονται από την έκθεση σε επίπεδο 9000 ng; Πού βρίσκονται οι πληγείσες πόλεις;
  2. Πόση αγροτική έκταση (\(km^2\)) επηρεάζεται από την έκθεση σε επίπεδο 8000 ng; Πού βρίσκονται αυτές οι περιοχές;
  3. Ποιες Διοικητικές Ενότητες έχουν μέσο όρο μόλυνσης μεγαλύτερο απο το όριο;

4 Δεδομένα

Ολα τα διανυσματικά δεδομένα ειναι μέσα στο αρχείο: ‘Vectors.gpkg’ Αυτό (*.gpkg) είναι μια μικρή βάση δεδομένων την οποία μπορείτε να ανοίξετε με οποιοδηποτε GIS λογισμικό. Εμείς θα διαβάσουμε τα διανυσματικά δεδομένα στην R.

Φορτώνουμε βιβλιοθηκες

library(sf)
library(terra)
library(dplyr)
#install.packages("remotes")
#library(remotes)
#install_github("r-tmap/tmap")

library(tmap)
library(ggpubr)

Ποιά διανυσματικά δεδομένα (vector) έχει μέσα η μικρή βάση δεδομένων μας;

st_layers("data/Exported/Vectors.gpkg")
Driver: GPKG 
Available layers:
  layer_name geometry_type features fields            crs_name
1       acto Multi Polygon        1      6 GGRS87 / Greek Grid
2        oik         Point      201     11 GGRS87 / Greek Grid
3   enotites Multi Polygon       73      2 GGRS87 / Greek Grid
4    nuclear         Point        1      6 GGRS87 / Greek Grid
5   stations         Point       30      1 GGRS87 / Greek Grid

διαβάσουμε τα vector δεδομένα στην R

acto = st_read("data/Exported/Vectors.gpkg", layer="acto", quiet = T )
oik = st_read("data/Exported/Vectors.gpkg", layer="oik", quiet = T )
nuclear = st_read("data/Exported/Vectors.gpkg", layer="nuclear", quiet = T )
#stations = st_read("data/Exported/Vectors.gpkg", layer="stations", quiet = T )
enotites= st_read("data/Exported/Vectors.gpkg", layer="enotites", quiet = T )

διαβάσουμε τα raster δεδομένα στην R

corine = rast("data/Exported/Corine.tif")
dem = rast("data/Exported/dem.tif") # Δεν θα χρησημοποιηθεί
ActualPollution = rast("data/Exported/ActualPollution.tif")

4.1 Xρήσεις γης (CORINE)

CORINE land cover types of the island grouped in 5 categories.

par(mar=c(1,1,1,3) )
mycol = c("grey30","orange","green","purple","blue")
plot(corine, legend = T, col = mycol)
#legend(x='bottomleft', xpd=T, inset=c(0.98, 0.40),
#       legend=levels(corine)[[1]]$value,
#       fill = mycol,bg =NULL, box.lwd = 0)
plot(st_geometry(acto), add=T, border="grey40")

#library(rasterVis)
#rasterVis::levelplot(corine, col.regions=mycol, xlab="", ylab="")

4.2 Οικισμοί

par(mar=c(1,1,1,1) )
plot(st_geometry(oik), main="Towns", pch=19, cex=0.6, col="red")
plot(st_geometry(acto), add=T, border="grey40")
plot(st_geometry(nuclear), add=T, pch=23, col="black", bg="black"  )
legend(x='bottomleft', pch=c(19,23), legend=c("Towns","Plant"),
       col=c("red","black"), pt.bg="black", bg=NULL, box.lwd=0)
box()

4.3 Αέρια μόλυνση

plot((ActualPollution), main="Actual Pollution")
plot(st_geometry(acto), add=T, border="grey40")
#plot(st_geometry(stations), add=T, pch=6, col="blue")
plot(st_geometry(nuclear), add=T, pch=23, col="black", bg="black")

legend(x='bottomleft', xpd=T, inset=c(0.10, 0.03),
       pch=c(6,23),
       legend=c("Nuclear Plant"),
       col = c("black"), pt.bg="black",
       bg =NULL, box.lwd = 0)

4.4 Όριο αέριας μόλυνσης

Το όριο της αέριας μόλυνσης

TheThreshold=8000

4.5 Διοικητικά όρια

plot(st_geometry(enotites), main="Administrative sub-units", border="grey40" )
plot(st_geometry(acto), add=T  )
text(st_coordinates(st_centroid( enotites)), label=enotites$kalcode,  cex=0.4,col="blue" ) #
box()

5 Ερώτημα 1: Πληθυσμός

Τι πληθυσμός επιρεάζεται απο το όριο των 8000 ή περισσότερο? Ποιοι οικισμοί επιρεάζονται?

Επιλογή ολων των οικισμών με μόλυνση >= 8000

oik$ActualPollution = terra::extract(ActualPollution, oik, ID=F)[[1]]
oik2 = subset(oik, ActualPollution >= TheThreshold)
dim(oik2)
[1] 46 13

Υπάρχουν 46 οικισμοί με έκθεση σε μόλυνση >= 8000. Αθροισμα πληθυσμού:

sum(oik2$pop2011)
[1] 41663

Οπτικοποιηση των 46 οικισμών με συνολικό πληθυσμό 41663 κατοίκους.

par(mar=c(1,1,2,1) )
plot(st_geometry(acto), main=sprintf("Οικισμοί σε επικύνδινη ζώνη\n>= %s ng",TheThreshold) )
plot(st_geometry(oik2), main="Towns", pch=19, cex=0.6, col="red", add=T)
plot(st_geometry(nuclear), add=T, pch=23, bg="black", col="black")
legend(x='bottomleft', pch=c(19,23), legend=c("Οικισμοι","Εργοστάσιο"),
       col=c("red","black"), pt.bg="black", bg=NULL, box.lwd=0)
box()

par(mar=c(0,0,0,0) )
terra::plot(ActualPollution, axes=F, box=F, main="Οικισμοί σε επικύνδινη ζώνη", plg=list(loc="bottom") )
plot(st_geometry(acto), main=sprintf("Οικισμοί σε επικύνδινη ζώνη\n>= %s ng",TheThreshold), add=T)
plot(st_geometry(oik),   pch=19, cex=0.5, col=rgb(red=0.1,green=0,blue=0,alpha=0.1), add=T)
plot(st_geometry(oik2),  pch=19, cex=0.6, col="red", add=T)
plot(st_geometry(nuclear), add=T, pch=23, bg="black", col="black")
legend(xpd=T, x='bottomleft', pch=c(19,23), legend=c("Οικισμοι","Εργοστάσιο"),  col=c("red","black"), pt.bg="black", bg=NULL, box.lwd=0)
box()

5.1 Απάντηση 1

Υπάρχουν 46 οικισμοί με συνολικό πληθυσμό 41663 κατοικων, που επιρεάζονται απο την μολυνση των 8000 ng.

6 Ερώτημα 2: Αγροτικές εκτάσεις

Πόσες αγροτικές εκτάσεις (σε \(χλμ^2\)) επιρεάζονται απο τη μόλυνση των 8000ng? Πού είναι αυτές οι εκτάσεις?

6.1 Επιλογή αγροτικών εκτάσεων

Επιλογή μόνο των αγροτικών εκτάσεων απο τα δεδομένα CORINE και μετατρωπη σε πολύγωνο.

corine
class       : SpatRaster 
dimensions  : 469, 686, 1  (nrow, ncol, nlyr)
resolution  : 101.0898, 101.0898  (x, y)
extent      : 656997.2, 726344.8, 4314945, 4362356  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=tmerc +lat_0=0 +lon_0=24 +k=0.9996 +x_0=500000 +y_0=0 +ellps=GRS80 +towgs84=-199.87,74.79,246.62,0,0,0,0 +units=m +no_defs 
source      : Corine.tif 
categories  : lesvoscorine 
name        : lesvoscorine 
min value   :   Artificial 
max value   :        Water

Ποιές κατηγορίες περιλαμβάνει αυτή η εκδοση του CORINE ?

levels(corine)[[1]]
  value lesvoscorine
1     1   Artificial
2     2 Agricultural
3     3       Forest
4     4     Wetlands
5     5        Water

Ποιό είναι το πληθος κελιών ανα κατηγορία?

#pol = rasterToPolygons(corine, fun=function(x){x==2})
table(corine[])

    1     2     3     4     5 
 5934 70333 81951   772  1386

Ποιο ειναι το συνολικό εμβαδό των αγροτικών εκτάσεων?

cellSize = cellSize(corine[[1]])[][[1]] / 1000
ola = sum(table(corine[]))
agr = sum((corine==2)[], na.rm=T)
emvado_agrotika = round((agr*cellSize) ,1)
emvado_agrotika
[1] 718870.3

Το ποσοστό των συνολικών αγροτικών κελιών:

pososto = 100*(agr/ola)
pososto
[1] 43.85507
res(corine==2) %>% prod() # Εναλακτικός τρόπος υπολογισμού του cell size (m)
[1] 10219.15

Οπτικοποίηση μόνο των αγροτικών χρήσεων από το CORINE:

plot( (corine==2), main=sprintf("Αγροτικές περιοχές (%s%% των κελιών)\n(%s χλμ^2)",
                   round(pososto,1), emvado_agrotika ))
plot(st_geometry(acto), add=T )

6.2 Επιλογή μεγάλης μόλυνσης

Απο το raster της μόλυνσης, επιλογή μονο των κελιών που ειναι πάνω απο το όριο τ

bigPolution = ( ActualPollution >= TheThreshold )
#bigPolution = terra::crop(  bigPolution, vect(acto) ) %>% mask(vect(acto))
#bigPolution = classify(bigPolution, cbind(c(TRUE, FALSE), c(1,NA))  )

Πόσα κελιά είνα πάνω απο το όριο της μέγάλης μολυνσης?

table(bigPolution[])

 FALSE   TRUE 
214347  35653

Ποιό είναι το εμβαδό της μεγάλης μολυνσης?

emvado_bigPollution = res(bigPolution) %>% prod()
(emvado_bigPollution)
[1] 13109.38

Οπτικοποιηση raster μεγάλης μολυνσης:

plot(bigPolution, main=sprintf("Μεγάλη μόλυνση (%s χλμ^2)", round(emvado_bigPollution ,1) ), legend=F, col=c(NA,"orange"))
plot(st_geometry(acto), add=T)
plot(st_geometry(nuclear), add=T, pch=23, bg="black", col="black")
legend(x='bottomleft', pch=c(19,23), legend=c("Big Poolution","Plant"),
       col=c("green4","black"), pt.bg="black", bg=NULL, box.lwd=0)

6.3 Μεγάλη μόλυνση + Αγρότικές εκτάσεις

Οπτικοποιηση των δυο raster.

par(mfrow = c(2, 1))
plot(corine==2, main="CORINE (Αγροτικά)", col=c(NA,"orange"))
plot(acto, add=T, col=NA)

plot(bigPolution, main="Μεγάλη μόλυνση", col=c(NA,"green"))
plot(acto, add=T, col=NA)

Ποιο ειναι το μέγεθος κελιου των δυο raster (resolution). Είναι διαφορετικό. Αρα θα χρειαστεί να κάνουμε επαναδειγματοληψία (resample) ώστε να έχουμε στη διάθεση μας δυο raster με ίδιο μέγεθος κελιου.

res(corine)
[1] 101.0898 101.0898
res(bigPolution)
[1] 138.5155  94.6420

Επαναδειγματοληψία (resample) CORINE στις διαστάσεις του bigPolution (raster μεγάλης μόλυνσης):

r2 = resample((corine==2), bigPolution, method='near')
plot(r2, main="Αγροτικά + Mόλυνση", col=c(NA,rgb(1,0.64,0,0.7)), legend=F )
plot(bigPolution,  col=c(NA,rgb(0,1,0,0.5)), add=T, legend=F)
plot(st_geometry( acto), add=T, col=NA)

Πληροφορίες κελιών του raster: CORINE

table((corine==2)[])

FALSE  TRUE 
90043 70333

Πληροφορίες κελιών του raster: resampled CORINE

table(r2[])

    0     1 
70072 54959

Πληροφορίες κελιών του raster: της μεγάλης μολυνσης

table(bigPolution[])

 FALSE   TRUE 
214347  35653 
#r2resampled <- projectRaster((corine==2),bigPolution, method = 'ngb')
#plot(r2resampled)

Συνθεση των δυο raster: Resampled CORINE + μεγάλη μολυνση

result = (r2 + (bigPolution))==2
plot(result , main="Αποτέλεσμα: Αγροτικά (CORINE)+\n μεγάλη μολυνση")

table(result[])

 FALSE   TRUE 
112468  12563

Πλήθος θετικών κελιών

ncells = sum((result[]==1), na.rm=T)
ncells
[1] 12563

Εμβαδό του κάθε κελιού

emvado_cell = res(result) %>% prod()
emvado_cell
[1] 13109.38

Συνολικό εμβαδό των θετικών κελιών

result_emvado = (ncells * emvado_cell)/1000
result_emvado
[1] 164693.2

6.4 Απάντηση 2

Από τα 718870.3 \(χλμ^2\) της συνολικής αγροτικής εκτασης, υπάρχουν 164693.1905764 \(χλμ^2\) τα οποία βρίσκονται σε περιοχή πάνω απο το όριο των 8000ng μολυνσης. Δηλαδή είναι το 23% της συνολικής αγροτικής εκτασης.

7 Ερώτημα 3. Διοικητικές ενότητες

Ερώτημα 3. Ποιες Διοικητικές Ενότητες έχουν μέσο όρο μόλυνσης μεγαλύτερο απο το όριο;

Για κάθε μια απο τις 73 Διοικητικές Ενοτητες, υπολογισμός της μολυνσης εντός του κάθε πολυγώνου:

  • Μέσος Ορος: meanPolution

  • Συνολικό εμβαδό Διοικητής Ενότητας: area

enotites$meanPolution = raster::extract(ActualPollution, enotites, fun=mean )$layer
enotites$area = st_area(enotites)

Ποιες Διοικητικές Ενοτητες έχουν ΜΟ μεγαλύτερο απο το όριο των 8000 ?

enotites$danger = ifelse(enotites$meanPolution>=TheThreshold, TRUE, FALSE)
ggbarplot(enotites %>% as_tibble(), x = "lektiko", y = "meanPolution",
          fill = "danger",
          color = "danger",
          sort.val = "asc",          # Sort the value in dscending order
          rotate = T,
          x.text.angle = 90,           # Rotate vertically x axis texts
          ggtheme = theme_bw(),
          palette = "Set2", legend = "bottom"
)

Απλοί θεματικοί χάρτες

qtm(enotites, fill = "meanPolution" )+ tm_layout(legend.position = c("right", "top")  )

qtm(enotites, fill = "danger" )+ tm_layout(legend.position = c("right", "top")  )

7.1 Απάντηση 3

Υπάρχουν 21 Διοικητικές Ενότητες με μεσο όρο μεγαλύτερο απο το όριο 8000.

8 Εργασία

Επαναλάβετε το παραπάνω tutorial χρησημοποιόντας νέα δεδομένα απο τον φάκελο: data/Exported2/ και θέτοντας όριο 7500ng. Απαντηστε στα παραπάνω 3 ερωτηματα σε μια αναφορά (1 pdf) η οποια να περιλαμβάνει: χάρτες, πινακες και διαγράμματα και κώδικα.