In etsakl/DasyMapR: Geographical ETRS-LAEA Grid Based Dasymetric Mapping
library(knitr)
library(DasyMapR)
opts_chunk$set(
tidy = TRUE,
cache = TRUE,
warning = FALSE,
message = FALSE
)
Το πακέτο DasyMapR
Το πακέτο της R που αναπτύχθηκε περιέχει μία σειρά εργαλείων με σκοπό την υλοποίηση δασυμετρικών χαρτών με την μέθοδο του επιμερισμού των δεδομένων στον κανονικό ETRS-LAEA κάναβο. Πληροφορίες για τον γεωγραφικό κανάβο μπορούν να βρεθούν INSPIRE Specification on Geographical Grid Systems
Τι περιέχεται σε αυτό το κείμενο
- Εγκατάσταση του πακέτου
- Τι υπάρχει στο πακέτο
- Φόρτωση Δεδομένων
- Απόδοση τιμής στο κελί - Κατηγορικά(ποιοτικά) Δεδομένα - MaxArea
- Απόδοση τιμής στο κελί - Αριθμητικά (ποσοτικά) Δεδομένα - PropCal
- Επιφάνεια Πηγή
- Βοηθητική Επιφάνεια
- Οι δασυμετρικοί υπολογισμοί
- Η Δασυμετρική επιφάνεια σε Βοηθητική
Εγκατάσταση
Το πακέτο δεν φιλοξενείται στο CRAN (θα χρειαστεί πολλή δουλειά ακόμη για αυτό) αλλά στο προσωπικό αποθετήριο του σπουδαστή στο github. Για να εγκατασταθεί το DasyMapR θα πρέπει καταρχήν να εγκατασταθεί το πακέτο devtools και στην συνέχεια με την χρήση της συνάρτησης install_github() να γίνει η εγκατάσταση του πακέτου
install.packages("devtools")
library(devtools)
install_github("etsakl/DasyMapR",build_vignettes=TRUE)
library(DasyMapR)
DasyMapR τι περιέχεται
Το πακέτο αποτελείται αποτελείται από από S4 - Classes από τις μεθόδους τους methods, την τεκμηρίωσή τους files.Rd και τα δεδομένα data που θα χρησιμοποιηθούν σε αυτό το κείμενο με την τεκμηρίωσή τους και το ίδιο το κείμενο ως vignette που επίσης συνοδεύει το πακέτο.
library(DasyMapR)
library(knitr)
DasyMapR.contains<-sort(ls("package:DasyMapR"))
kable(as.data.frame(DasyMapR.contains))
Περισσότερες πληροφορίες για το πακέτο θα μπορούσε να βρει κάποιος στα file του πακέτου help("DasyMapR") ή στο github
Φόρτωση Δεδομένων
Το πακέτο συνοδεύεται από κάποια δεδομένα κατά κανόνα ανοικτά δημόσια δεδομένα που η προσπάθεια για να συμβάλλουμε στην αξιοποίηση τους έγινε και η αφορμή για ανάπτυξη του συγκεκριμένου πακέτου που σκοπεύει να συμβάλει στην οπτικοποίηση και ανάλυσή τους. Προφανώς για λόγους οικονομίας χρόνου μεταφόρτωσης και αποθηκευτικού χώρου τα δεδομένα αυτά είναι υποσύνολα των δεδομένων που προσφέρονται από τους οργανισμούς που τα διαθέτουν. Τεκμηρίωση και πληροφορίες υπάρχουν στα help files του πακέτου help("DasyMapR")` ή στο github
# Αν τρέξει ο κώδικας παρουσιάζονται τα διαθέσιμα δεδομένα
DasyMapR.data<-data(package ="DasyMapR")
print(DasyMapR.data)
Θα φορτώσουμε τα δεδομένα που μεταμορφώθηκαν από το GeoData.open.gov και θα επιλέξουμε τους νομούς της περιφέρειας Πελοποννήσου για τους δασυμετρικούς μας υπολογισμούς.
# Φόρτωσε τα δεδομένα
data("NUTS3_OCMG")
# Δίάλεξε κάποια
pp <- c("Ν. ΑΡΓΟΛΙΔΟΣ","Ν. ΚΟΡΙΝΘΙΑΣ","Ν. ΑΡΚΑΔΙΑΣ","Ν. ΜΕΣΣΗΝΙΑΣ","Ν. ΛΑΚΩΝΙΑΣ")
Π.ΠΕΛΟΠΟΝΝΗΣΟΥ <- NUTS3_OCMG[which(!is.na(match(NUTS3_OCMG[["NAME"]],pp))),]
par("mar"=c(.1,.1,.1,.1))
# και σχεδίασε τα
plot(NUTS3_OCMG)
plot(Π.ΠΕΛΟΠΟΝΝΗΣΟΥ,border=2,lwd=2,add=T)
plot(Π.ΠΕΛΟΠΟΝΝΗΣΟΥ,border=2,lwd=2)
ας δούμε και τα περιγραφικά δεδομένα των νομών που περιέχει το αρχείο.
# Τα περιγράφικά σε πίνακα
kable(head(Π.ΠΕΛΟΠΟΝΝΗΣΟΥ@data[,3:10]))
Κατηγορικά Δεδομένα - Η χρήση της MaxArea
Η Νομοί έχουν ένα ποιοτικό χαρακτηριστικό αυτό της "Ονομασίας" Π.ΠΕΛΟΠΟΝΝΗΣΟΥ[["ΝΑΜΕ"]] που τους προσδίδει την ιδιότητα του Νομού (περιοχή διοικητικής διαίρεσης). Θα ήταν χρήσιμο αυτή η ιδιότητα να αποδοθεί στα κελιά του κανάβου αφού με βάση αυτή την ιδιότητα στη συνέχεια θα μπορούσαν να του π.χ. αποδοθούν και άλλα χαρακτηριστικά που έχουν ως επιφάνεια απαρίθμησης, αυτή ακριβώς την ιδιότητα δηλαδή του νόμου. Το κελί σε αυτήν την περίπτωση των κατηγορικών δεδομένων παίρνει την τιμή της επιφάνειας που καταλαμβάνει το μεγαλύτερο μέρος του κελιού. Η μέθοδος etrsSurface που δέχεται ορίσματα την επιφάνεια input.surface την μέθοδο Max.Area και το μέγεθος του κελιού cell.size. Θα πρέπει να αναφερθεί εδώ ότι η επιφάνεια NUTS3_OCMG είναι σε σύστημα συντεταγμένων ΕΓΣΑ87(GGRS'87 (EPSG:2100)) όμως τα αποτελέσματα δηλαδή η τελικά υπολογιζόμενη etrsSurface θα έχει πλέον σύστημα αναφοράς το ETRS-LAEA καλώντας κατά τους υπολογισμούς την DasyMapR::etsrTransform() του πακέτου και στέλνει μήνυμα warning ενημερώνοντας για την μετατροπή.
library(DasyMapR)
par("mar"=c(.1,.1,.1,.1))
# Καλεί την EtrsTransForm
srf.grd.max <- etrsSurface(input.surface = Π.ΠΕΛΟΠΟΝΝΗΣΟΥ,over.method.type = "MaxArea",cell.size = 10000)
plot(Π.ΠΕΛΟΠΟΝΝΗΣΟΥ)
plot(srf.grd.max)
Το αντικείμενο που δημιουργήθηκε δηλαδή srf.grd.max είναι αντικείμενο κλάσης EtrsSurface και τα περιγραφικά δεδομένα στο slot (@data) είναι σε μορφή πίνακα :
kable(head(srf.grd.max@data))
Όπως φαίνεται στο πίνακα τα στοιχεία που περιέχονται είναι ο κωδικός $CELLCODEτου κελιού οι συντεταγμένες του κάτω άκρου του κελιού $EASTOFORIGIN,$NORTHOFORIGIN" και ο αριθμός αναγνώρισης της επιφάνειας $FEATURE από την οποία έχει πάρει την ιδιότητα το κελί. Ο κωδικός του κελιού που επαναλαμβάνεται στην αρχή είναι η ιδιότητα row.names() του data.frame @data και χρησιμεύει ως αναγνωριστικό σύνδεσης με το γεωμετρικό μέρος του αντικείμενου. (@polygons). Περισσότερες πληροφορίες για την δομή του αντικειμένου EtrsSurface αλλά και του SpatialPolygonsDataFrame παρέχονται στο κυρίως κείμενο και στα help files.
Θα ήταν χρήσιμο να συμπεριλάβει κάποιος ιδιότητες από τα περιγραφικά δεδομένα της επιφάνειας πηγής στο πίνακα δεδομένων των κελιών. H μέθοδος DasyMapR::joinMaxAreaSurfaceDataFrame μπορεί να κληθεί και να συνενώσει τα χαρακτηριστικά που δεν έχουν συμπεριληφθεί αρχικά ως ιδιότητες της νέας επιφάνειας (του κανάβου).
srf.grd.max.full<-joinMaxAreaSurfaceDataFrames(the.surface = NUTS3_OCMG, the.EtrsSurface = srf.grd.max)
kable(head(srf.grd.max.full@data[,c(7:12)]))
Κάτι που πρέπει να επισημανθεί εδώ είναι ότι ο χρήστης πρέπει πάντα να έχει υπόψη του τι δεδομένα έχει και τι αναγωγές πρέπει να κάνει για να έχει ορθά αποτελέσματα. Για παράδειγμα στον τελευταίο πίνακα στο κελί ως τιμές του πληθυσμού εμφανίζονται οι $POP91 ή $pl2001. Αυτό είναι λάθος διότι η τιμές αυτές αναφέρονται στις αρχικές επιφάνειες απαριθμήσεις δηλαδή στους νομούς. Ακόμη και η τιμή της πυκνότητας ανά κελί είναι εσφαλμένη αφού δεν έχουν γίνει κατάλληλες αναγωγές.
Ο χρόνος που χρειάζεται για να γίνουν οι υπολογισμοί μπορούν να μειωθεί με την χρήση της μεθόδου DasyMapR::etrsSurfacePar() όπου γίνεται χρήση περισσοτέρων πόρων του συστήματος (επεξεργαστών και πυρήνων) με την χρήση των πακέτων foreach και doParaller της R. Περισσότερες πληροφορίες για αυτό μπορούν να βρεθούν στο κυρίως κείμενο.
Αριθμητικά (ποσοτικά) Δεδομένα - Η χρήση της PropCal
Για την καλύτερη παρακολούθηση του παραδείγματος και κάνοντας χρίση των ίδιων δεδομένων που αντλήσαμε από το GeoData.open.gov θα κρατήσουμε μόνο τις στήλες που είναι χρήσιμες για το τρέχον παράδειγμά.
# Θα κρατήσουμε μόνο τον πλυθισμό του 2001
Π.ΠΕΛΟΠΟΝΝΗΣΟΥ.2001<-Π.ΠΕΛΟΠΟΝΝΗΣΟΥ[,c(3,7)]
kable(Π.ΠΕΛΟΠΟΝΝΗΣΟΥ.2001@data)
Η τιμή στόχος που θέλουμε να επιμερίσουμε στον κάναβο είναι αυτή του πληθυσμού. Αν επιμερίσουμε αυτή την τιμή του χαρακτηριστικού pl$2001 θα αποδώσουμε στο κελί λανθασμένη τιμή (όπως προαναφέρθηκε) . Θα πρέπει να μετατραπεί από απόλυτη τιμή σε πυκνότητα πληθυσμού. Στο πακέτο έχει αναπτυχθεί η μέθοδος ():DasyMapR::ActullVal2Density() που θα κάνει την μετατροπή
Π.ΠΕΛΟΠΟΝΝΗΣΟΥ.2001<-ActuallVal2Density(input.surface = Π.ΠΕΛΟΠΟΝΝΗΣΟΥ.2001,surface.value.col = 2,area.unit = 1e+06)
kable(Π.ΠΕΛΟΠΟΝΝΗΣΟΥ.2001@data)
Η τιμή $VALUE πλέον αναφέρεται σε κατοίκους/km^2^ και εκφράζει πυκνότητα. Το @data$VALUE (slot) πλέον περιέχει και την τιμή του εμβαδού σε km^2^.
Θα επιμερίσουμε την τιμή με την χρήση της μεθόδου etrsSourceSurface() για την τιμή (VALUE) της πυκνότητας Η μέθοδος που είναι καταλληλότερη για τον υπολογισμό της τιμής του κελιού μετά την προβολή της επιφάνειας θα είναι αυτή του αναλογικού υπολογισμού σε σχέση με την επιφάνεια $CELLVALUE = (V~i~ ∗ Share~i~)$ όπου V~i~, η τιμή της επιφάνειας και Share~i~, το μερίδιο της επιφάνειας i μέσα στο κελί
srf.grd.prop<-etrsSourceSurface(input.surface = Π.ΠΕΛΟΠΟΝΝΗΣΟΥ.2001,over.method.type = "PropCal",surface.value.col = 4,cell.size = 10000)
par("mar"=c(.1,.1,.1,.1))
plot(Π.ΠΕΛΟΠΟΝΝΗΣΟΥ)
plot(srf.grd.prop)
Το αντικείμενο που δημιουργήθηκε είναι αντικείμενο της κλάσης EtrsSourceSurface και ο πίνακας που είναι συσχετισμένος με την επιφάνεια περιεχέι τα χαρακτηριστικά
kable(head(srf.grd.prop@data))
Λεπτομέρειες για τον υπολογισμό της $CELLVALUE μπορεί κάποιος να βρει αν φορτώσει το προσωρινό αρχείο με όνομα ".surface.detailed.table.rds" που αποθηκεύετε στο working directory. Για το παράδειγμα της συγκεκριμένης εργασίας φορτώνουμε το εν λόγω αρχείο (.rds) και εμφανίζουμε ορισμένες γραμμές του.Σε αυτές εμφανίζονται οι τιμές για τμήμα του κελιού πριν την άθροιση και τον τελικό υπολογισμό (για αυτό και εμφανίζονται row.names(srf.grd.prop@data) με την προσθήκη αριθμού ώστε να είναι μοναδικός ο κωδικός αριθμός ID.
wd<-getwd()
setwd(wd)
surface.detailed.table<-readRDS(".surface.detailed.table.rds")
kable(head(surface.detailed.table,5))
Για πολλά δεδομένα προτείνεται να γίνει χρήση των προσφερόμενων parallel computing μεθόδων
Ανακεφαλαιώνοντας - Η επιφάνεια Εισόδου
Ίσως να μην έχει γίνει κατανοητό πως οι διαφορετικές μέθοδοι επηρεάζουν την γεωγραφική απεικόνιση της επιφάνειας πηγής για αυτό στο επόμενο σχήμα απεικονίζονται οι δύο περιπτώσεις εφαρμογής των δύο μεθόδων.
# Η κλήση της EtrsTrasnform για αλλαγη του CRS
Π.ΠΕΛΟΠΟΝΝΗΣΟΥ.ETRS <- EtrsTransform(Π.ΠΕΛΟΠΟΝΝΗΣΟΥ)
par("mar"=c(.1,.1,.1,.1))
plot(Π.ΠΕΛΟΠΟΝΝΗΣΟΥ.ETRS,border=2,lwd=2)
plot(srf.grd.max,col=rgb(0,1,0,0.1),add=TRUE)
plot(Π.ΠΕΛΟΠΟΝΝΗΣΟΥ.ETRS,border=2,lwd=2)
plot(srf.grd.prop,col=rgb(0,1,0,0.1),add=TRUE)
Στην πρώτη περίπτωση η επιφάνεια μπορεί να μην καλύπτεται ή υπερκαλύπτεται ενώ στην δεύτερη πάντα υπερκαλύπτεται. Είναι διακριτό ότι η μέθοδος που επιλέγεται είναι σημαντική για την απόκρυψη ή μη πληροφορίας και την παραγωγή σφαλμάτων (στην πρώτη εικόνα τμήματα ξηράς εξαιρούνται από την ανάλυση μας στην δεύτερη τμήματα της λογίζονται πλέον ως ξηρά). Περισσότερα σχετικά με τα σφάλματα της μεθόδου μπορεί να αναζητηθούν εδω
Για να είναι "ρεαλιστικότερο" το παράδειγμα μας θα αλλάξουμε επίπεδο αναφοράς και από το επίπεδο της περιφέρειας θα περάσουμε στο επίπεδο της περιφερειακής ενότητας (νομού ή NUTS3) και της ανάλυσης από τον κάναβο του 10km στον κάναβο 1km.
Επιλέγεται εδώ ο νόμος Αργολίδας
par(mar = c(0.1, 0.1, 0.1, 0.1))
ΑΡΓΟΛΙΔΑ <- NUTS3_OCMG[which(!is.na(match(NUTS3_OCMG[["NAME"]], "Ν. ΑΡΓΟΛΙΔΟΣ"))),]
plot(Π.ΠΕΛΟΠΟΝΝΗΣΟΥ)
plot(ΑΡΓΟΛΙΔΑ, add = TRUE, lwd = 2, border = 2)
# Εδώ καλείται η EtrsTransform απευθείας
ΑΡΓΟΛΙΔΑ.ETRS<-EtrsTransform(ΑΡΓΟΛΙΔΑ)
# Με την κλήση της etrsSourceSurface παράγεται η επιφανεια πηγή
source.surface <- etrsSourceSurface(input.surface = Π.ΠΕΛΟΠΟΝΝΗΣΟΥ.2001[3,],
over.method.type = "PropCal", surface.value.col = 4, cell.size = 1000)
plot(source.surface,col=rgb(0,1,0,0.01),lwd=.5,border="lightgrey")
plot(ΑΡΓΟΛΙΔΑ.ETRS, add = TRUE, lwd = 2, border = 2)
Η Βοηθιτική επιφάνεια
Θα γίνει χρήση των δεδομένων CORINE 2000 δηλαδή των δεδομένων κάλυψη γης για την Ελλάδα και το έτος 2000, όπως διατίθενται από Ευρωπαϊκή Υπηρεσία Περιβάλλοντος . Τα δεδομένα αυτά μπορούν να μεταμορφωθούν από τον ισότοπο της ΕΕΑ και επιλέχθηκαν ως κατάλληλα από την άποψη της χρονικής κατανομής σε σχέση με το φαινόμενο (πληθυσμό 2001) που θέλουμε να ανακατανείμουμε.
Για την προετοιμασία της επιφάνειας επίσης θα γίνει η χρήση της etrsAncillarySurface. Εδώ κρίνεται σκόπιμο να γίνουν ορισμένες επισημάνσεις στο χρήστη του πακέτου. Παρόλο που οι τρεις μέθοδοι δημιουργίας των επιφανειών etrsSurface, etsrSourceSurface, etrsAncillarySurface επί της ουσίας κάνουν ακριβώς τους ίδιους υπολογισμούς, και οι κλάσεις περιέχουν τα ίδια @slots και θα μπορούσε ο χρήστης να αποφασίζει για την "ορθή" χρήση τους η επιλογή αυτή έγινε για να καθοδηγήσει τον χρήστη για το πώς θα πρέπει να ενεργήσει. Θα γίνουμε πιο συγκεκριμένη στην πορεία ολοκλήρωσης του τρέχοντος παραδείγματος.
Προετοιμασία των Δεδομένων
Εδώ ο χρήστης θα πρέπει να κάνει τις επιλογές που στηρίζονται στην εμπιρεία του η/και την γνώση του για το θέμα (Ως εκ τούτου η μέθοδος συχνά αναφέρεται ως Intelligent Dasymetric Maping. Ας δούμε τις κατηγορίες καλύψεων που περιέχονται στην Corine DB:
data("CLC2000_CODES")
kable(head(CLC2000_CODES))
...
Είναι προφανές ότι για το παράδειγμα, ενδιαφέρουν μόνο οι δύο κατηγορίες 111 και 112 που σχετίζονται με την κατοικία. Αυτές θα μεταφορτώνουν από την ΕΕΑ αφού διατίθενται για μεταφόρτωση σε αρχεία .shp κατα τις προαναφερόμενες κατηγόρίες.
Με τον παρακάτω κώδικα θα μπορούσε κάποιος να δημιουργησει αν dataset με κάλυψη για όποια κατηγορία ή έτος απαιτείται για την ανάλυσή του. Παρατίθεται εδώ για να αποτελέσει οδηγό σε ενδεχομένη προετοιμασία άλλων δεδομένων από τον χρήστη.
# Τα δύο .shp files περιέχονται στο folder corine
setwd(system.file("data/corine",package="DasyMapR"))
# Με βάση τα όρια της περιοχής ...
bb<-bbox(ΑΡΓΟΛΙΔΑ.ETRS)
# διμιουργησε νέα .shp files που περιέχουν όσα δεδομένα χρειάζομαι
ogr2ogr(".","clc_cliped",spat = c(bb[,1],bb[,2]))
dsn<-setwd("clc_cliped/")
# φόρτωσε τα ως SpatialPolygonsDataFrames
CLC2000_POLY_ARGOLIDA111<-readOGR(".","clc00_v2_code_111")
CLC2000_POLY_ARGOLIDA112<-readOGR(".","clc00_v2_code_112")
# και ένωσε τα
CLC2000.ARGOLIDA.RES<-rbind.SpatialPolygonsDataFrame(CLC2000_POLY_ARGOLIDA111,CLC2000_POLY_ARGOLIDA112,makeUniqueIDs = T)
# τέλος σώστα στο δίσκο ως dataset
setwd(system.file("data",package = "DasyMapR"))
# Αφαιρέθηκαν 2 πολύγωνα με πρόβλημα στην γεωμετρία
CLC2000.ARGOLIDA.RES<-CLC2000.ARGOLIDA.RES[-which(row.names(CLC2000.ARGOLIDA.RES)==8),]
CLC2000.ARGOLIDA.RES<-CLC2000.ARGOLIDA.RES[-which(row.names(CLC2000.ARGOLIDA.RES)%in% "01"),]
devtools::use_data(CLC2000.ARGOLIDA.RES,overwrite = T)
Στο πακέτο DasyMapR υπάρχουν ήδη τα δεδομένα αυτά στο φάκελο .\data ως datasets που το συνοδεύουν.
#Φορτώνουμε τα δεδομένα
data("CLC2000.ARGOLIDA.RES")
# Τι Περιέχει το σετ;
kable(head(CLC2000.ARGOLIDA.RES@data,3))
Και σε πιο άνθρωπο-φιλική μορφή
ΠΕΡΙΟΧΕΣ.ΚΑΤΟΙΚΙΑΣ <- merge(x = CLC2000.ARGOLIDA.RES,y = CLC2000_CODES,by.x="code_00",by.y="Code_00")
kable(head(ΠΕΡΙΟΧΕΣ.ΚΑΤΟΙΚΙΑΣ@data,3))
Θα θεωρήσουμε ότι ο λόγος της πυκνότητας $111 : 112 = 3 : 1$. δηλαδή ότι στις περιοχές με πυκνό αστικό ιστό ζουν 3 φορές περισσότερου άνθρωποι ανά μονάδα επιφανείας από ότι στις αραιά δομημένες περιοχές.
par("mar"=c(.1,.1,.1,.1))
ReDens111<-round(3/4,2)
ReDens112<-round(1/4,2)
ΠΕΡΙΟΧΕΣ.ΚΑΤΟΙΚΙΑΣ@data[which(ΠΕΡΙΟΧΕΣ.ΚΑΤΟΙΚΙΑΣ@data[,"code_00"] == 111),"ReDens"]<-ReDens111
ΠΕΡΙΟΧΕΣ.ΚΑΤΟΙΚΙΑΣ@data[which(ΠΕΡΙΟΧΕΣ.ΚΑΤΟΙΚΙΑΣ@data[,"code_00"] == 112),"ReDens"]<-ReDens112
kable(head(ΠΕΡΙΟΧΕΣ.ΚΑΤΟΙΚΙΑΣ@data))
plot(ΑΡΓΟΛΙΔΑ.ETRS,lwd=2,border=2)
plot(ΠΕΡΙΟΧΕΣ.ΚΑΤΟΙΚΙΑΣ,col="purple",add=TRUE)
data("NUTSV9_LEAC")
plot(NUTSV9_LEAC,add=TRUE,border="lightgrey")
Καλώντας την etrsAncillarySurface θα υπολογίσουμε την "σχετική πυκνότητα" του κάθε κελιού πλέον με βάση τις παραδοχές που έγιναν: Όλος ο πλυθησμός της Αργολίδας το έτος 2001 κατοικουσε σις περιοχές με συνεχή και διακεκκομενη αστική δόμηση και με σχετική πυκνότητα μάλιστα 1:3. Θα μπορούσαμε να συμπεριλάβουμε και άλλες περιοχές και να του αποδώσουμε κάποιο ποσοστό. (π.χ αγροτικές περιοχές 5 %)
par("mar"=c(.1,.1,.1,.1))
data("NUTSV9_LEAC")
the.ancillary.surface.bf <- etrsAncillarySurface(input.surface = ΠΕΡΙΟΧΕΣ.ΚΑΤΟΙΚΙΑΣ,over.method.type = "PropCal",surface.value.col = 3,cell.size = 1000,binary = FALSE)
plot(the.ancillary.surface.bf,col= "purple")
plot(ΑΡΓΟΛΙΔΑ.ETRS, add = TRUE, lwd = 2, border = 2)
plot(NUTSV9_LEAC,add=TRUE,border="lightgrey")
kable(head(the.ancillary.surface.bf@data))
surface.detailed.table.bf<-readRDS(".surface.detailed.table.rds")
Θα επαναλάβουμε του ίδιους υπολογισμούς θέτοντας την παράμετρο binary = TRUE η οποία δίνει διαφοροποιημένα αποτελέσματα. Για να έχουμε αποτελέσματα θα πρέπει να θέσουμε την σχετική πυκνότητα σε 1. Στην ουσία λειτουργεί σαν την εφαρμογή της MaxArea στον υπολογισμό της βοηθητικής επιφάνειας. Πρακτικά σε φυσικούς όρους σημαίνει αν ένα κελί έχει την ιδιότητα να φιλοξενεί πληθυσμό με τιμή μεγαλύτερη από 50% Στην απεικόνιση που ακολουθεί μόνο τα κελιά με μοβ χρώμα έχουν τιμή για την $WCELLWEGHT
par("mar"=c(.1,.1,.1,.1))
ΠΕΡΙΟΧΕΣ.ΚΑΤΟΙΚΙΑΣ@data[which(ΠΕΡΙΟΧΕΣ.ΚΑΤΟΙΚΙΑΣ@data[,"code_00"] == c(111,112)),"ReDens"]<-1
the.ancillary.surface.bt <- etrsAncillarySurface(input.surface = ΠΕΡΙΟΧΕΣ.ΚΑΤΟΙΚΙΑΣ,over.method.type = "PropCal",surface.value.col = 3,cell.size = 1000,binary = TRUE)
plot(the.ancillary.surface.bt,col= "purple4")
plot(the.ancillary.surface.bt[which(the.ancillary.surface.bt[["WCELLWEIGHT"]]==0),],col= "lightgrey",add=T)
plot(ΑΡΓΟΛΙΔΑ.ETRS, add = TRUE, lwd = 2, border = 2)
plot(NUTSV9_LEAC,add=TRUE,border="lightgrey")
kable(head(the.ancillary.surface.bt@data))
surface.detailed.table.bt<-readRDS(".surface.detailed.table.rds")
Ας δούμε και τα αποτελέσματά του πίνακα καθώς και του προσωρινού αρχείου για τον έλεγχο τον υπολογισμών
wd<-getwd()
setwd(wd)
kable(head(surface.detailed.table.bf))
Ανακεφαλαιώνοντας
Είναι κατανοητό ότι και εδώ τα ζητήματα της κλίμακας του φαινομένου (κάλυψη) και της ανάλυσης (μέγεθος κελιού 1km) είναι σημαντικό να λαμβάνονται υπόψη από τον χρήστη για την παραγωγή αποτελεσμάτων αξιοποιήσιμων και "ρεαλιστικών" ώστε να διατηρείται η πληροφορία και να μην παράγεται νέα αλλά ψευδής.
Για την συνέχεια του παραδείγματος θα επιλέξουμε την κατανομή που προέκυψε από την εφαρμογή της etrsAncillarySurface(...,"PropCAl",...) δηλαδή την the.ancillary.surface.bf
Οι δασυμετρικοί υπολογισμοί
H μέθοδος για την εφαρμογή των δασυμετρικών υπολογισμών που θα κλιθεί είναι etsrDasymetricSurface και δέχεται ως ορίσματα μία επιφάνεια πηγή κλάσης EtrsSourceSurfaceτην βοηθητική επιφάνεια κλάσης EtrsAncillarySurface και μία μεταβλητή logical actuall.value όπου εάν δοθεί η τιμή TRUE θα εφαρμοστεί η ActuallVal2Density() που η χρήση της εξηγήθηκε παραπάνω
par("mar"=c(.1,.1,.1,.1))
dasymetric.surface<-EtrsDasymetricSurface(input.surface.grided = source.surface,ancillary.grided =the.ancillary.surface.bf ,actuall.value = FALSE)
plot(dasymetric.surface)
plot(ΑΡΓΟΛΙΔΑ.ETRS,add=T,border="red",lwd=2)
kable(head(dasymetric.surface@data))
Η χρήση του NUTS και LAEA κανάβου
Το παράδειγμα που χρησιμοποιήθηκε όπως υπονοήθηκε δεν είναι αρκετά ακριβές αφού δεν χορηγήσαμε επαρκή στοιχεία τόσο για την κατανομή του πληθυσμού όσο και λόγω της μικρής κλίμακας των δεδομένων της κάλυψης. Θα μπορούσαμε να χρησιμοποιήσουμε δεδομένα από τις UMZ200
Ειδικά για το πληθυσμό προσφέρονται δεδομένα για το πληθυσμό στον ETRS-LAEA κάναβο για τα έτη 2006 2011 που έχουν προκύψει από την εφαρμογή πιο περίπλοκων αλγορίθμων μπορούν να στο ESPOn.
data("GEOSTAT_grid_EU_POP_2006_1k_V1_1_1")
kable(head(GEOSTAT_grid_EU_POP_2006_1k_V1_1_1))
GR_POP_2006<-GEOSTAT_grid_EU_POP_2006_1k_V1_1_1[which(GEOSTAT_grid_EU_POP_2006_1k_V1_1_1[,'CNTR_CODE'] %in% "EL"),]
Μία πιο προσεκτική ματιά σε αυτά τα δεδομένα δείχνει δεν ακολουθούν την προτεινόμενη κωδικοποίηση των κελιών της INSPIRE Specification on Geographical Grid Systems Για να τέτοια προβλήματα αναπτύχθηκε η μέθοδος etrsReverseCellCode.
GR_POP_2006<-as.data.frame(GR_POP_2006)
GR_POP_2006<-etrsReverseCellCode(df = GR_POP_2006,cell.code.col = 1)
kable(head(GR_POP_2006))
Και για να απεικονισθούν τα δεδομένα μας στη μορφή του κανάβου μπορούμε να καλέσουμε την EtrsGrid και στην συνέχεια θα συγχωνεύσουμε τα πολύγωνα με τις εγγραφές στο data frame με του κωδικούς των κελλιών (row.names,"ID"(),ID).
par("mar"=c(.1,.1,.1,.1))
GR251<-NUTSV9_LEAC[which(!is.na(match(NUTSV9_LEAC[["N3CD"]],"GR251"))),]
GR251<-EtrsTransform(GR251)
GR251.grd<-etrsGrid(GR251,cell.size = 1000)
GR251.grd<-merge(GR251.grd,GR_POP_2006,by=0,all=F)
plot(GR251.grd)
plot(ΑΡΓΟΛΙΔΑ.ETRS,add=TRUE,border=2,lwd=3)
<<<<<<< HEAD
Η Δασυμετρική επιφάνεια σε Βοηθητική
=======
Η Δασυμετρική σε Βοηθητική
8e7ff89b1d44f91b64c476d81766c5eed543a545
Όπως αναφέρθηκε τα διοικητικά όρια συχνά αλλάζουν τέτοιο παράδειγμα είναι η περίπτωση του Ευρωπαϊκού χώρου όπου το ιστορικό των αλλαγων συχνά δυσκολεύει την χρήση και σύγκριση των διαθέσιμων δεδομένων. Τα γεωγραφικά όρια μπορούν να μεταμορφωθούν από το της EUROSTAT
Ας δούμε κατ' αρχήν τις αλλαγές των ορίων στην υπό μελέτη περιοχή την Αργολίδα :
Φορτώνουμε τα δεδομένα και τα :
par("mar"=c(.1,.1,.1,.1))
# Όρια NUTS 2006
GR25<-NUTSV9_LEAC[which(!is.na(match(NUTSV9_LEAC[["N2CD"]],"GR25"))),]
kable(head(GR25@data))
#Νέα όρια NUTS 2013
data("NUTS_2013_01M_EL")
kable(head(NUTS_2013_01M_EL@data))
NUTS_2013_01M_65<-EtrsTransform(NUTS_2013_01M_EL[grep('^EL65',NUTS_2013_01M_EL[["NUTS_ID"]]) ,])
plot(NUTS_2013_01M_65)
plot(GR25)
Για όποιον ενδιαφέρεται για δεδομένα κοινωνικό-οικονομικά, περιβαλλοντικά κ.λ.π. ήδη γνωρίζει ότι μια σημαντική πηγή στοιχείων είναι EurostatDatabase Με την χρήση του package eurostat μπορούμε να συνεχίσουμε την ανάλυση μας και να απεικονίσουμε κοινωνικό - οικονομικά δεδομένα στον κάναβο. Π.Χ. μπορούμε να κατεβάσουμε δεδομένα για το ακαθάριστο εθνικό προϊόν ανά κάτοικο και να τα προβάλουμε στο κελί επιμερίζοντας τα με βάση τον αριθμό κατοίκων ανά κελί από τα προηγούμενα βήματα.
Ας αναζητήσουμε κάποια στοιχεία από την Eurostat. Θα κάνουμε την αναζήτηση με βάση τo γεωγραφικό επίπεδο αναφοράς (NUTS3)
info<-search_eurostat("NUTS 3")
kable(info[c(23:30),c(1,2)])
επιλέγονται δεδομένα για το ακαθάριστο εθνικό προϊόν του 2001 (Gross domestic product (GDP) at current market prices by NUTS 3 regions) και εiδικότερα για την περιοχή ενδιαφέροντος. H Επόμενη γραμμή χρειάζεται σύνδεση και δεν θα τρέξει αν δεν υπάρχει οπότε στο παράδειγμά μας
#nama_10r_3gdp <- get_eurostat(id = "nama_10r_3gdp" ,filters = list(time=2006),time_format = "num")
Τα δεδομένα αυτά έχουν ήδη και περιέχονται στο πακέτο. Έτσι ι τα φορτώνονται από το /data
με την χρήση data()
data("nama_10r_3gdp")
dat<-nama_10r_3gdp
α περιορίσουμε την εργασία εδώ για λόγους οικονομίας στον νομό Αργολίδας (και "αναγκαστικά" και Αργολίδας)
GDP651_2006<-dat[grep('^EL651',dat$geo),]
kable(label_eurostat(GDP651_2006))
Για να γίνει γεωναφορά των δεδομένων αυτών θα χρησιμοποιήσουμε τα όρια όπως αυτά παρέχονται EUROSTAT.
par("mar"=c(.1,.1,.1,.1))
NUTS_2013_01M_651<-EtrsTransform(NUTS_2013_01M_EL[grep('^EL651',NUTS_2013_01M_EL[["NUTS_ID"]]) ,])
kable(NUTS_2013_01M_651@data)
Αρχικά γίνεται η γεωαναφορά των στοιχείων
NUTS_2013_01M_651<-merge(NUTS_2013_01M_651,GDP651_2006[1,],by.x="NUTS_ID",by.y="geo",all=FALSE)
H χρήση της merge ενδεχομένως να προκαλέσει αλλαγή στα rownames του dataframe οπότε έστω και προληπτικά τα διορθώνουμε
row.names(NUTS_2013_01M_651@data)<-sapply(slot(NUTS_2013_01M_651, "polygons"), function(x) slot(x, "ID"))
plot(NUTS_2013_01M_651)
kable(NUTS_2013_01M_651@data)
Στην συνέχεια θα χρησιμοποιήσουμε την etrsSourceSurface για να προβάλλουμε στον καναβο το "φαινόμενο" της κατανομής του κατά κεφαλήν ακαθάριστου εθνικού προϊόντος
par("mar"=c(.1,.1,.1,.1))
NUTS_2013_01M_651_GDP<-etrsSourceSurface(input.surface = NUTS_2013_01M_651,over.method.type = "PropCal",surface.value.col = 7,cell.size = 1000)
plot(NUTS_2013_01M_651_GDP)
kable(head(NUTS_2013_01M_651_GDP@data))
Η προηγούμενη επιφάνεια dasymetric.surface που υπολογίστηκε ως Δασυμετρική επιφάνεια τώρα θα μετατραπεί σε βοηθητική επιφάνεια και θα χρησιμοποιεί για την κατανομή της τιμής της GDP. Σε αυτήν απεικονιζόταν ο αριθμός των κατοίκων ανά κελί που είχε την ιδιότητα του συνεχούς ή αστικού ιστού
kable(head(dasymetric.surface@data))
POP_2001_ancillary<- etrsDasymetric2Ancillary(dasymetric.surface)
row.names(POP_2001_ancillary@data)<-sapply(slot(POP_2001_ancillary, "polygons"), function(x) slot(x, "ID"))
Τέλος με την κλίση της etrsPropWeightedValue θα υπολογιστεί η τιμή του ακαθάριστου προϊόντος σε κάθε κελί
par("mar"=c(0.1,0.1,0.1,0.1))
DASY_GPD<-etrsPropWeightedValue(input.surface.grided = NUTS_2013_01M_651_GDP, ancillary.grided = POP_2001_ancillary)
plot(DASY_GPD)
plot(ΑΡΓΟΛΙΔΑ.ETRS,add=T,border=2,lwd=2)
kable(head(DASY_GPD))
Τέλος η μετατροπή των αποτελέσματα σε raster θα έκανε τα ευκολότερα διαθέσιμα για χρήση. Θα χρησιμοποιήσουμε την etrsDasymetric2Rasterπου εξαρτάται απο το πακέτο raster για να επιτύχουμε το ζητούμενο αποτέλεσμα
par("mar"=c(0.1,0.1,0.1,0.1))
DASY_GPD_RASTER<-etrsDasymetric2Raster(dasymetric.surface = DASY_GPD)
rw.colors<-grey.colors
image(DASY_GPD_RASTER,col=rw.colors(5))
plot(ΑΡΓΟΛΙΔΑ.ETRS,add=T,border=2,lwd=2)
etsakl/DasyMapR documentation built on May 16, 2019, 9:07 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
library(knitr) library(DasyMapR) opts_chunk$set( tidy = TRUE, cache = TRUE, warning = FALSE, message = FALSE )
Το πακέτο DasyMapR
Το πακέτο της R που αναπτύχθηκε περιέχει μία σειρά εργαλείων με σκοπό την υλοποίηση δασυμετρικών χαρτών με την μέθοδο του επιμερισμού των δεδομένων στον κανονικό ETRS-LAEA κάναβο. Πληροφορίες για τον γεωγραφικό κανάβο μπορούν να βρεθούν INSPIRE Specification on Geographical Grid Systems
Τι περιέχεται σε αυτό το κείμενο
- Εγκατάσταση του πακέτου
- Τι υπάρχει στο πακέτο
- Φόρτωση Δεδομένων
- Απόδοση τιμής στο κελί - Κατηγορικά(ποιοτικά) Δεδομένα - MaxArea
- Απόδοση τιμής στο κελί - Αριθμητικά (ποσοτικά) Δεδομένα - PropCal
- Επιφάνεια Πηγή
- Βοηθητική Επιφάνεια
- Οι δασυμετρικοί υπολογισμοί
- Η Δασυμετρική επιφάνεια σε Βοηθητική
Εγκατάσταση
Το πακέτο δεν φιλοξενείται στο CRAN (θα χρειαστεί πολλή δουλειά ακόμη για αυτό) αλλά στο προσωπικό αποθετήριο του σπουδαστή στο github. Για να εγκατασταθεί το DasyMapR θα πρέπει καταρχήν να εγκατασταθεί το πακέτο devtools και στην συνέχεια με την χρήση της συνάρτησης install_github() να γίνει η εγκατάσταση του πακέτου
install.packages("devtools") library(devtools) install_github("etsakl/DasyMapR",build_vignettes=TRUE) library(DasyMapR)
DasyMapR τι περιέχεται
Το πακέτο αποτελείται αποτελείται από από S4 - Classes από τις μεθόδους τους methods, την τεκμηρίωσή τους files.Rd και τα δεδομένα data που θα χρησιμοποιηθούν σε αυτό το κείμενο με την τεκμηρίωσή τους και το ίδιο το κείμενο ως vignette που επίσης συνοδεύει το πακέτο.
library(DasyMapR) library(knitr) DasyMapR.contains<-sort(ls("package:DasyMapR")) kable(as.data.frame(DasyMapR.contains))
Περισσότερες πληροφορίες για το πακέτο θα μπορούσε να βρει κάποιος στα file του πακέτου help("DasyMapR") ή στο github
Φόρτωση Δεδομένων
Το πακέτο συνοδεύεται από κάποια δεδομένα κατά κανόνα ανοικτά δημόσια δεδομένα που η προσπάθεια για να συμβάλλουμε στην αξιοποίηση τους έγινε και η αφορμή για ανάπτυξη του συγκεκριμένου πακέτου που σκοπεύει να συμβάλει στην οπτικοποίηση και ανάλυσή τους. Προφανώς για λόγους οικονομίας χρόνου μεταφόρτωσης και αποθηκευτικού χώρου τα δεδομένα αυτά είναι υποσύνολα των δεδομένων που προσφέρονται από τους οργανισμούς που τα διαθέτουν. Τεκμηρίωση και πληροφορίες υπάρχουν στα help files του πακέτου help("DasyMapR")` ή στο github
# Αν τρέξει ο κώδικας παρουσιάζονται τα διαθέσιμα δεδομένα DasyMapR.data<-data(package ="DasyMapR") print(DasyMapR.data)
Θα φορτώσουμε τα δεδομένα που μεταμορφώθηκαν από το GeoData.open.gov και θα επιλέξουμε τους νομούς της περιφέρειας Πελοποννήσου για τους δασυμετρικούς μας υπολογισμούς.
# Φόρτωσε τα δεδομένα data("NUTS3_OCMG") # Δίάλεξε κάποια pp <- c("Ν. ΑΡΓΟΛΙΔΟΣ","Ν. ΚΟΡΙΝΘΙΑΣ","Ν. ΑΡΚΑΔΙΑΣ","Ν. ΜΕΣΣΗΝΙΑΣ","Ν. ΛΑΚΩΝΙΑΣ") Π.ΠΕΛΟΠΟΝΝΗΣΟΥ <- NUTS3_OCMG[which(!is.na(match(NUTS3_OCMG[["NAME"]],pp))),] par("mar"=c(.1,.1,.1,.1)) # και σχεδίασε τα plot(NUTS3_OCMG) plot(Π.ΠΕΛΟΠΟΝΝΗΣΟΥ,border=2,lwd=2,add=T) plot(Π.ΠΕΛΟΠΟΝΝΗΣΟΥ,border=2,lwd=2)
ας δούμε και τα περιγραφικά δεδομένα των νομών που περιέχει το αρχείο.
# Τα περιγράφικά σε πίνακα kable(head(Π.ΠΕΛΟΠΟΝΝΗΣΟΥ@data[,3:10]))
Κατηγορικά Δεδομένα - Η χρήση της MaxArea
Η Νομοί έχουν ένα ποιοτικό χαρακτηριστικό αυτό της "Ονομασίας" Π.ΠΕΛΟΠΟΝΝΗΣΟΥ[["ΝΑΜΕ"]] που τους προσδίδει την ιδιότητα του Νομού (περιοχή διοικητικής διαίρεσης). Θα ήταν χρήσιμο αυτή η ιδιότητα να αποδοθεί στα κελιά του κανάβου αφού με βάση αυτή την ιδιότητα στη συνέχεια θα μπορούσαν να του π.χ. αποδοθούν και άλλα χαρακτηριστικά που έχουν ως επιφάνεια απαρίθμησης, αυτή ακριβώς την ιδιότητα δηλαδή του νόμου. Το κελί σε αυτήν την περίπτωση των κατηγορικών δεδομένων παίρνει την τιμή της επιφάνειας που καταλαμβάνει το μεγαλύτερο μέρος του κελιού. Η μέθοδος etrsSurface που δέχεται ορίσματα την επιφάνεια input.surface την μέθοδο Max.Area και το μέγεθος του κελιού cell.size. Θα πρέπει να αναφερθεί εδώ ότι η επιφάνεια NUTS3_OCMG είναι σε σύστημα συντεταγμένων ΕΓΣΑ87(GGRS'87 (EPSG:2100)) όμως τα αποτελέσματα δηλαδή η τελικά υπολογιζόμενη etrsSurface θα έχει πλέον σύστημα αναφοράς το ETRS-LAEA καλώντας κατά τους υπολογισμούς την DasyMapR::etsrTransform() του πακέτου και στέλνει μήνυμα warning ενημερώνοντας για την μετατροπή.
library(DasyMapR) par("mar"=c(.1,.1,.1,.1)) # Καλεί την EtrsTransForm srf.grd.max <- etrsSurface(input.surface = Π.ΠΕΛΟΠΟΝΝΗΣΟΥ,over.method.type = "MaxArea",cell.size = 10000) plot(Π.ΠΕΛΟΠΟΝΝΗΣΟΥ) plot(srf.grd.max)
Το αντικείμενο που δημιουργήθηκε δηλαδή srf.grd.max είναι αντικείμενο κλάσης EtrsSurface και τα περιγραφικά δεδομένα στο slot (@data) είναι σε μορφή πίνακα :
kable(head(srf.grd.max@data))
Όπως φαίνεται στο πίνακα τα στοιχεία που περιέχονται είναι ο κωδικός $CELLCODEτου κελιού οι συντεταγμένες του κάτω άκρου του κελιού $EASTOFORIGIN,$NORTHOFORIGIN" και ο αριθμός αναγνώρισης της επιφάνειας $FEATURE από την οποία έχει πάρει την ιδιότητα το κελί. Ο κωδικός του κελιού που επαναλαμβάνεται στην αρχή είναι η ιδιότητα row.names() του data.frame @data και χρησιμεύει ως αναγνωριστικό σύνδεσης με το γεωμετρικό μέρος του αντικείμενου. (@polygons). Περισσότερες πληροφορίες για την δομή του αντικειμένου EtrsSurface αλλά και του SpatialPolygonsDataFrame παρέχονται στο κυρίως κείμενο και στα help files.
Θα ήταν χρήσιμο να συμπεριλάβει κάποιος ιδιότητες από τα περιγραφικά δεδομένα της επιφάνειας πηγής στο πίνακα δεδομένων των κελιών. H μέθοδος DasyMapR::joinMaxAreaSurfaceDataFrame μπορεί να κληθεί και να συνενώσει τα χαρακτηριστικά που δεν έχουν συμπεριληφθεί αρχικά ως ιδιότητες της νέας επιφάνειας (του κανάβου).
srf.grd.max.full<-joinMaxAreaSurfaceDataFrames(the.surface = NUTS3_OCMG, the.EtrsSurface = srf.grd.max) kable(head(srf.grd.max.full@data[,c(7:12)]))
Κάτι που πρέπει να επισημανθεί εδώ είναι ότι ο χρήστης πρέπει πάντα να έχει υπόψη του τι δεδομένα έχει και τι αναγωγές πρέπει να κάνει για να έχει ορθά αποτελέσματα. Για παράδειγμα στον τελευταίο πίνακα στο κελί ως τιμές του πληθυσμού εμφανίζονται οι $POP91 ή $pl2001. Αυτό είναι λάθος διότι η τιμές αυτές αναφέρονται στις αρχικές επιφάνειες απαριθμήσεις δηλαδή στους νομούς. Ακόμη και η τιμή της πυκνότητας ανά κελί είναι εσφαλμένη αφού δεν έχουν γίνει κατάλληλες αναγωγές.
Ο χρόνος που χρειάζεται για να γίνουν οι υπολογισμοί μπορούν να μειωθεί με την χρήση της μεθόδου DasyMapR::etrsSurfacePar() όπου γίνεται χρήση περισσοτέρων πόρων του συστήματος (επεξεργαστών και πυρήνων) με την χρήση των πακέτων foreach και doParaller της R. Περισσότερες πληροφορίες για αυτό μπορούν να βρεθούν στο κυρίως κείμενο.
Αριθμητικά (ποσοτικά) Δεδομένα - Η χρήση της PropCal
Για την καλύτερη παρακολούθηση του παραδείγματος και κάνοντας χρίση των ίδιων δεδομένων που αντλήσαμε από το GeoData.open.gov θα κρατήσουμε μόνο τις στήλες που είναι χρήσιμες για το τρέχον παράδειγμά.
# Θα κρατήσουμε μόνο τον πλυθισμό του 2001 Π.ΠΕΛΟΠΟΝΝΗΣΟΥ.2001<-Π.ΠΕΛΟΠΟΝΝΗΣΟΥ[,c(3,7)] kable(Π.ΠΕΛΟΠΟΝΝΗΣΟΥ.2001@data)
Η τιμή στόχος που θέλουμε να επιμερίσουμε στον κάναβο είναι αυτή του πληθυσμού. Αν επιμερίσουμε αυτή την τιμή του χαρακτηριστικού pl$2001 θα αποδώσουμε στο κελί λανθασμένη τιμή (όπως προαναφέρθηκε) . Θα πρέπει να μετατραπεί από απόλυτη τιμή σε πυκνότητα πληθυσμού. Στο πακέτο έχει αναπτυχθεί η μέθοδος ():DasyMapR::ActullVal2Density() που θα κάνει την μετατροπή
Π.ΠΕΛΟΠΟΝΝΗΣΟΥ.2001<-ActuallVal2Density(input.surface = Π.ΠΕΛΟΠΟΝΝΗΣΟΥ.2001,surface.value.col = 2,area.unit = 1e+06) kable(Π.ΠΕΛΟΠΟΝΝΗΣΟΥ.2001@data)
Η τιμή $VALUE πλέον αναφέρεται σε κατοίκους/km^2^ και εκφράζει πυκνότητα. Το @data$VALUE (slot) πλέον περιέχει και την τιμή του εμβαδού σε km^2^.
Θα επιμερίσουμε την τιμή με την χρήση της μεθόδου etrsSourceSurface() για την τιμή (VALUE) της πυκνότητας Η μέθοδος που είναι καταλληλότερη για τον υπολογισμό της τιμής του κελιού μετά την προβολή της επιφάνειας θα είναι αυτή του αναλογικού υπολογισμού σε σχέση με την επιφάνεια $CELLVALUE = (V~i~ ∗ Share~i~)$ όπου V~i~, η τιμή της επιφάνειας και Share~i~, το μερίδιο της επιφάνειας i μέσα στο κελί
srf.grd.prop<-etrsSourceSurface(input.surface = Π.ΠΕΛΟΠΟΝΝΗΣΟΥ.2001,over.method.type = "PropCal",surface.value.col = 4,cell.size = 10000) par("mar"=c(.1,.1,.1,.1)) plot(Π.ΠΕΛΟΠΟΝΝΗΣΟΥ) plot(srf.grd.prop)
Το αντικείμενο που δημιουργήθηκε είναι αντικείμενο της κλάσης EtrsSourceSurface και ο πίνακας που είναι συσχετισμένος με την επιφάνεια περιεχέι τα χαρακτηριστικά
kable(head(srf.grd.prop@data))
Λεπτομέρειες για τον υπολογισμό της $CELLVALUE μπορεί κάποιος να βρει αν φορτώσει το προσωρινό αρχείο με όνομα ".surface.detailed.table.rds" που αποθηκεύετε στο working directory. Για το παράδειγμα της συγκεκριμένης εργασίας φορτώνουμε το εν λόγω αρχείο (.rds) και εμφανίζουμε ορισμένες γραμμές του.Σε αυτές εμφανίζονται οι τιμές για τμήμα του κελιού πριν την άθροιση και τον τελικό υπολογισμό (για αυτό και εμφανίζονται row.names(srf.grd.prop@data) με την προσθήκη αριθμού ώστε να είναι μοναδικός ο κωδικός αριθμός ID.
wd<-getwd() setwd(wd) surface.detailed.table<-readRDS(".surface.detailed.table.rds") kable(head(surface.detailed.table,5))
Για πολλά δεδομένα προτείνεται να γίνει χρήση των προσφερόμενων parallel computing μεθόδων
Ανακεφαλαιώνοντας - Η επιφάνεια Εισόδου
Ίσως να μην έχει γίνει κατανοητό πως οι διαφορετικές μέθοδοι επηρεάζουν την γεωγραφική απεικόνιση της επιφάνειας πηγής για αυτό στο επόμενο σχήμα απεικονίζονται οι δύο περιπτώσεις εφαρμογής των δύο μεθόδων.
# Η κλήση της EtrsTrasnform για αλλαγη του CRS Π.ΠΕΛΟΠΟΝΝΗΣΟΥ.ETRS <- EtrsTransform(Π.ΠΕΛΟΠΟΝΝΗΣΟΥ) par("mar"=c(.1,.1,.1,.1)) plot(Π.ΠΕΛΟΠΟΝΝΗΣΟΥ.ETRS,border=2,lwd=2) plot(srf.grd.max,col=rgb(0,1,0,0.1),add=TRUE) plot(Π.ΠΕΛΟΠΟΝΝΗΣΟΥ.ETRS,border=2,lwd=2) plot(srf.grd.prop,col=rgb(0,1,0,0.1),add=TRUE)
Στην πρώτη περίπτωση η επιφάνεια μπορεί να μην καλύπτεται ή υπερκαλύπτεται ενώ στην δεύτερη πάντα υπερκαλύπτεται. Είναι διακριτό ότι η μέθοδος που επιλέγεται είναι σημαντική για την απόκρυψη ή μη πληροφορίας και την παραγωγή σφαλμάτων (στην πρώτη εικόνα τμήματα ξηράς εξαιρούνται από την ανάλυση μας στην δεύτερη τμήματα της λογίζονται πλέον ως ξηρά). Περισσότερα σχετικά με τα σφάλματα της μεθόδου μπορεί να αναζητηθούν εδω Για να είναι "ρεαλιστικότερο" το παράδειγμα μας θα αλλάξουμε επίπεδο αναφοράς και από το επίπεδο της περιφέρειας θα περάσουμε στο επίπεδο της περιφερειακής ενότητας (νομού ή NUTS3) και της ανάλυσης από τον κάναβο του 10km στον κάναβο 1km. Επιλέγεται εδώ ο νόμος Αργολίδας
par(mar = c(0.1, 0.1, 0.1, 0.1)) ΑΡΓΟΛΙΔΑ <- NUTS3_OCMG[which(!is.na(match(NUTS3_OCMG[["NAME"]], "Ν. ΑΡΓΟΛΙΔΟΣ"))),] plot(Π.ΠΕΛΟΠΟΝΝΗΣΟΥ) plot(ΑΡΓΟΛΙΔΑ, add = TRUE, lwd = 2, border = 2) # Εδώ καλείται η EtrsTransform απευθείας ΑΡΓΟΛΙΔΑ.ETRS<-EtrsTransform(ΑΡΓΟΛΙΔΑ) # Με την κλήση της etrsSourceSurface παράγεται η επιφανεια πηγή source.surface <- etrsSourceSurface(input.surface = Π.ΠΕΛΟΠΟΝΝΗΣΟΥ.2001[3,], over.method.type = "PropCal", surface.value.col = 4, cell.size = 1000) plot(source.surface,col=rgb(0,1,0,0.01),lwd=.5,border="lightgrey") plot(ΑΡΓΟΛΙΔΑ.ETRS, add = TRUE, lwd = 2, border = 2)
Η Βοηθιτική επιφάνεια
Θα γίνει χρήση των δεδομένων CORINE 2000 δηλαδή των δεδομένων κάλυψη γης για την Ελλάδα και το έτος 2000, όπως διατίθενται από Ευρωπαϊκή Υπηρεσία Περιβάλλοντος . Τα δεδομένα αυτά μπορούν να μεταμορφωθούν από τον ισότοπο της ΕΕΑ και επιλέχθηκαν ως κατάλληλα από την άποψη της χρονικής κατανομής σε σχέση με το φαινόμενο (πληθυσμό 2001) που θέλουμε να ανακατανείμουμε.
Για την προετοιμασία της επιφάνειας επίσης θα γίνει η χρήση της etrsAncillarySurface. Εδώ κρίνεται σκόπιμο να γίνουν ορισμένες επισημάνσεις στο χρήστη του πακέτου. Παρόλο που οι τρεις μέθοδοι δημιουργίας των επιφανειών etrsSurface, etsrSourceSurface, etrsAncillarySurface επί της ουσίας κάνουν ακριβώς τους ίδιους υπολογισμούς, και οι κλάσεις περιέχουν τα ίδια @slots και θα μπορούσε ο χρήστης να αποφασίζει για την "ορθή" χρήση τους η επιλογή αυτή έγινε για να καθοδηγήσει τον χρήστη για το πώς θα πρέπει να ενεργήσει. Θα γίνουμε πιο συγκεκριμένη στην πορεία ολοκλήρωσης του τρέχοντος παραδείγματος.
Προετοιμασία των Δεδομένων
Εδώ ο χρήστης θα πρέπει να κάνει τις επιλογές που στηρίζονται στην εμπιρεία του η/και την γνώση του για το θέμα (Ως εκ τούτου η μέθοδος συχνά αναφέρεται ως Intelligent Dasymetric Maping. Ας δούμε τις κατηγορίες καλύψεων που περιέχονται στην Corine DB:
data("CLC2000_CODES") kable(head(CLC2000_CODES))
...
Είναι προφανές ότι για το παράδειγμα, ενδιαφέρουν μόνο οι δύο κατηγορίες 111 και 112 που σχετίζονται με την κατοικία. Αυτές θα μεταφορτώνουν από την ΕΕΑ αφού διατίθενται για μεταφόρτωση σε αρχεία .shp κατα τις προαναφερόμενες κατηγόρίες.
Με τον παρακάτω κώδικα θα μπορούσε κάποιος να δημιουργησει αν dataset με κάλυψη για όποια κατηγορία ή έτος απαιτείται για την ανάλυσή του. Παρατίθεται εδώ για να αποτελέσει οδηγό σε ενδεχομένη προετοιμασία άλλων δεδομένων από τον χρήστη.
# Τα δύο .shp files περιέχονται στο folder corine setwd(system.file("data/corine",package="DasyMapR")) # Με βάση τα όρια της περιοχής ... bb<-bbox(ΑΡΓΟΛΙΔΑ.ETRS) # διμιουργησε νέα .shp files που περιέχουν όσα δεδομένα χρειάζομαι ogr2ogr(".","clc_cliped",spat = c(bb[,1],bb[,2])) dsn<-setwd("clc_cliped/") # φόρτωσε τα ως SpatialPolygonsDataFrames CLC2000_POLY_ARGOLIDA111<-readOGR(".","clc00_v2_code_111") CLC2000_POLY_ARGOLIDA112<-readOGR(".","clc00_v2_code_112") # και ένωσε τα CLC2000.ARGOLIDA.RES<-rbind.SpatialPolygonsDataFrame(CLC2000_POLY_ARGOLIDA111,CLC2000_POLY_ARGOLIDA112,makeUniqueIDs = T) # τέλος σώστα στο δίσκο ως dataset setwd(system.file("data",package = "DasyMapR")) # Αφαιρέθηκαν 2 πολύγωνα με πρόβλημα στην γεωμετρία CLC2000.ARGOLIDA.RES<-CLC2000.ARGOLIDA.RES[-which(row.names(CLC2000.ARGOLIDA.RES)==8),] CLC2000.ARGOLIDA.RES<-CLC2000.ARGOLIDA.RES[-which(row.names(CLC2000.ARGOLIDA.RES)%in% "01"),] devtools::use_data(CLC2000.ARGOLIDA.RES,overwrite = T)
Στο πακέτο DasyMapR υπάρχουν ήδη τα δεδομένα αυτά στο φάκελο .\data ως datasets που το συνοδεύουν.
#Φορτώνουμε τα δεδομένα data("CLC2000.ARGOLIDA.RES") # Τι Περιέχει το σετ; kable(head(CLC2000.ARGOLIDA.RES@data,3))
Και σε πιο άνθρωπο-φιλική μορφή
ΠΕΡΙΟΧΕΣ.ΚΑΤΟΙΚΙΑΣ <- merge(x = CLC2000.ARGOLIDA.RES,y = CLC2000_CODES,by.x="code_00",by.y="Code_00") kable(head(ΠΕΡΙΟΧΕΣ.ΚΑΤΟΙΚΙΑΣ@data,3))
Θα θεωρήσουμε ότι ο λόγος της πυκνότητας $111 : 112 = 3 : 1$. δηλαδή ότι στις περιοχές με πυκνό αστικό ιστό ζουν 3 φορές περισσότερου άνθρωποι ανά μονάδα επιφανείας από ότι στις αραιά δομημένες περιοχές.
par("mar"=c(.1,.1,.1,.1)) ReDens111<-round(3/4,2) ReDens112<-round(1/4,2) ΠΕΡΙΟΧΕΣ.ΚΑΤΟΙΚΙΑΣ@data[which(ΠΕΡΙΟΧΕΣ.ΚΑΤΟΙΚΙΑΣ@data[,"code_00"] == 111),"ReDens"]<-ReDens111 ΠΕΡΙΟΧΕΣ.ΚΑΤΟΙΚΙΑΣ@data[which(ΠΕΡΙΟΧΕΣ.ΚΑΤΟΙΚΙΑΣ@data[,"code_00"] == 112),"ReDens"]<-ReDens112 kable(head(ΠΕΡΙΟΧΕΣ.ΚΑΤΟΙΚΙΑΣ@data)) plot(ΑΡΓΟΛΙΔΑ.ETRS,lwd=2,border=2) plot(ΠΕΡΙΟΧΕΣ.ΚΑΤΟΙΚΙΑΣ,col="purple",add=TRUE) data("NUTSV9_LEAC") plot(NUTSV9_LEAC,add=TRUE,border="lightgrey")
Καλώντας την etrsAncillarySurface θα υπολογίσουμε την "σχετική πυκνότητα" του κάθε κελιού πλέον με βάση τις παραδοχές που έγιναν: Όλος ο πλυθησμός της Αργολίδας το έτος 2001 κατοικουσε σις περιοχές με συνεχή και διακεκκομενη αστική δόμηση και με σχετική πυκνότητα μάλιστα 1:3. Θα μπορούσαμε να συμπεριλάβουμε και άλλες περιοχές και να του αποδώσουμε κάποιο ποσοστό. (π.χ αγροτικές περιοχές 5 %)
par("mar"=c(.1,.1,.1,.1)) data("NUTSV9_LEAC") the.ancillary.surface.bf <- etrsAncillarySurface(input.surface = ΠΕΡΙΟΧΕΣ.ΚΑΤΟΙΚΙΑΣ,over.method.type = "PropCal",surface.value.col = 3,cell.size = 1000,binary = FALSE) plot(the.ancillary.surface.bf,col= "purple") plot(ΑΡΓΟΛΙΔΑ.ETRS, add = TRUE, lwd = 2, border = 2) plot(NUTSV9_LEAC,add=TRUE,border="lightgrey")
kable(head(the.ancillary.surface.bf@data)) surface.detailed.table.bf<-readRDS(".surface.detailed.table.rds")
Θα επαναλάβουμε του ίδιους υπολογισμούς θέτοντας την παράμετρο binary = TRUE η οποία δίνει διαφοροποιημένα αποτελέσματα. Για να έχουμε αποτελέσματα θα πρέπει να θέσουμε την σχετική πυκνότητα σε 1. Στην ουσία λειτουργεί σαν την εφαρμογή της MaxArea στον υπολογισμό της βοηθητικής επιφάνειας. Πρακτικά σε φυσικούς όρους σημαίνει αν ένα κελί έχει την ιδιότητα να φιλοξενεί πληθυσμό με τιμή μεγαλύτερη από 50% Στην απεικόνιση που ακολουθεί μόνο τα κελιά με μοβ χρώμα έχουν τιμή για την $WCELLWEGHT
par("mar"=c(.1,.1,.1,.1)) ΠΕΡΙΟΧΕΣ.ΚΑΤΟΙΚΙΑΣ@data[which(ΠΕΡΙΟΧΕΣ.ΚΑΤΟΙΚΙΑΣ@data[,"code_00"] == c(111,112)),"ReDens"]<-1 the.ancillary.surface.bt <- etrsAncillarySurface(input.surface = ΠΕΡΙΟΧΕΣ.ΚΑΤΟΙΚΙΑΣ,over.method.type = "PropCal",surface.value.col = 3,cell.size = 1000,binary = TRUE) plot(the.ancillary.surface.bt,col= "purple4") plot(the.ancillary.surface.bt[which(the.ancillary.surface.bt[["WCELLWEIGHT"]]==0),],col= "lightgrey",add=T) plot(ΑΡΓΟΛΙΔΑ.ETRS, add = TRUE, lwd = 2, border = 2) plot(NUTSV9_LEAC,add=TRUE,border="lightgrey")
kable(head(the.ancillary.surface.bt@data)) surface.detailed.table.bt<-readRDS(".surface.detailed.table.rds")
Ας δούμε και τα αποτελέσματά του πίνακα καθώς και του προσωρινού αρχείου για τον έλεγχο τον υπολογισμών
wd<-getwd() setwd(wd) kable(head(surface.detailed.table.bf))
Ανακεφαλαιώνοντας
Είναι κατανοητό ότι και εδώ τα ζητήματα της κλίμακας του φαινομένου (κάλυψη) και της ανάλυσης (μέγεθος κελιού 1km) είναι σημαντικό να λαμβάνονται υπόψη από τον χρήστη για την παραγωγή αποτελεσμάτων αξιοποιήσιμων και "ρεαλιστικών" ώστε να διατηρείται η πληροφορία και να μην παράγεται νέα αλλά ψευδής.
Για την συνέχεια του παραδείγματος θα επιλέξουμε την κατανομή που προέκυψε από την εφαρμογή της etrsAncillarySurface(...,"PropCAl",...) δηλαδή την the.ancillary.surface.bf
Οι δασυμετρικοί υπολογισμοί
H μέθοδος για την εφαρμογή των δασυμετρικών υπολογισμών που θα κλιθεί είναι etsrDasymetricSurface και δέχεται ως ορίσματα μία επιφάνεια πηγή κλάσης EtrsSourceSurfaceτην βοηθητική επιφάνεια κλάσης EtrsAncillarySurface και μία μεταβλητή logical actuall.value όπου εάν δοθεί η τιμή TRUE θα εφαρμοστεί η ActuallVal2Density() που η χρήση της εξηγήθηκε παραπάνω
par("mar"=c(.1,.1,.1,.1)) dasymetric.surface<-EtrsDasymetricSurface(input.surface.grided = source.surface,ancillary.grided =the.ancillary.surface.bf ,actuall.value = FALSE) plot(dasymetric.surface) plot(ΑΡΓΟΛΙΔΑ.ETRS,add=T,border="red",lwd=2) kable(head(dasymetric.surface@data))
Η χρήση του NUTS και LAEA κανάβου
Το παράδειγμα που χρησιμοποιήθηκε όπως υπονοήθηκε δεν είναι αρκετά ακριβές αφού δεν χορηγήσαμε επαρκή στοιχεία τόσο για την κατανομή του πληθυσμού όσο και λόγω της μικρής κλίμακας των δεδομένων της κάλυψης. Θα μπορούσαμε να χρησιμοποιήσουμε δεδομένα από τις UMZ200
Ειδικά για το πληθυσμό προσφέρονται δεδομένα για το πληθυσμό στον ETRS-LAEA κάναβο για τα έτη 2006 2011 που έχουν προκύψει από την εφαρμογή πιο περίπλοκων αλγορίθμων μπορούν να στο ESPOn.
data("GEOSTAT_grid_EU_POP_2006_1k_V1_1_1") kable(head(GEOSTAT_grid_EU_POP_2006_1k_V1_1_1)) GR_POP_2006<-GEOSTAT_grid_EU_POP_2006_1k_V1_1_1[which(GEOSTAT_grid_EU_POP_2006_1k_V1_1_1[,'CNTR_CODE'] %in% "EL"),]
Μία πιο προσεκτική ματιά σε αυτά τα δεδομένα δείχνει δεν ακολουθούν την προτεινόμενη κωδικοποίηση των κελιών της INSPIRE Specification on Geographical Grid Systems Για να τέτοια προβλήματα αναπτύχθηκε η μέθοδος etrsReverseCellCode.
GR_POP_2006<-as.data.frame(GR_POP_2006) GR_POP_2006<-etrsReverseCellCode(df = GR_POP_2006,cell.code.col = 1) kable(head(GR_POP_2006))
Και για να απεικονισθούν τα δεδομένα μας στη μορφή του κανάβου μπορούμε να καλέσουμε την EtrsGrid και στην συνέχεια θα συγχωνεύσουμε τα πολύγωνα με τις εγγραφές στο data frame με του κωδικούς των κελλιών (row.names,"ID"(),ID).
par("mar"=c(.1,.1,.1,.1)) GR251<-NUTSV9_LEAC[which(!is.na(match(NUTSV9_LEAC[["N3CD"]],"GR251"))),] GR251<-EtrsTransform(GR251) GR251.grd<-etrsGrid(GR251,cell.size = 1000) GR251.grd<-merge(GR251.grd,GR_POP_2006,by=0,all=F) plot(GR251.grd) plot(ΑΡΓΟΛΙΔΑ.ETRS,add=TRUE,border=2,lwd=3)
<<<<<<< HEAD
Η Δασυμετρική επιφάνεια σε Βοηθητική
=======
Η Δασυμετρική σε Βοηθητική
8e7ff89b1d44f91b64c476d81766c5eed543a545
Όπως αναφέρθηκε τα διοικητικά όρια συχνά αλλάζουν τέτοιο παράδειγμα είναι η περίπτωση του Ευρωπαϊκού χώρου όπου το ιστορικό των αλλαγων συχνά δυσκολεύει την χρήση και σύγκριση των διαθέσιμων δεδομένων. Τα γεωγραφικά όρια μπορούν να μεταμορφωθούν από το της EUROSTAT Ας δούμε κατ' αρχήν τις αλλαγές των ορίων στην υπό μελέτη περιοχή την Αργολίδα :
Φορτώνουμε τα δεδομένα και τα :
par("mar"=c(.1,.1,.1,.1)) # Όρια NUTS 2006 GR25<-NUTSV9_LEAC[which(!is.na(match(NUTSV9_LEAC[["N2CD"]],"GR25"))),] kable(head(GR25@data)) #Νέα όρια NUTS 2013 data("NUTS_2013_01M_EL") kable(head(NUTS_2013_01M_EL@data)) NUTS_2013_01M_65<-EtrsTransform(NUTS_2013_01M_EL[grep('^EL65',NUTS_2013_01M_EL[["NUTS_ID"]]) ,]) plot(NUTS_2013_01M_65) plot(GR25)
Για όποιον ενδιαφέρεται για δεδομένα κοινωνικό-οικονομικά, περιβαλλοντικά κ.λ.π. ήδη γνωρίζει ότι μια σημαντική πηγή στοιχείων είναι EurostatDatabase Με την χρήση του package eurostat μπορούμε να συνεχίσουμε την ανάλυση μας και να απεικονίσουμε κοινωνικό - οικονομικά δεδομένα στον κάναβο. Π.Χ. μπορούμε να κατεβάσουμε δεδομένα για το ακαθάριστο εθνικό προϊόν ανά κάτοικο και να τα προβάλουμε στο κελί επιμερίζοντας τα με βάση τον αριθμό κατοίκων ανά κελί από τα προηγούμενα βήματα.
Ας αναζητήσουμε κάποια στοιχεία από την Eurostat. Θα κάνουμε την αναζήτηση με βάση τo γεωγραφικό επίπεδο αναφοράς (NUTS3)
info<-search_eurostat("NUTS 3") kable(info[c(23:30),c(1,2)])
επιλέγονται δεδομένα για το ακαθάριστο εθνικό προϊόν του 2001 (Gross domestic product (GDP) at current market prices by NUTS 3 regions) και εiδικότερα για την περιοχή ενδιαφέροντος. H Επόμενη γραμμή χρειάζεται σύνδεση και δεν θα τρέξει αν δεν υπάρχει οπότε στο παράδειγμά μας
#nama_10r_3gdp <- get_eurostat(id = "nama_10r_3gdp" ,filters = list(time=2006),time_format = "num")
Τα δεδομένα αυτά έχουν ήδη και περιέχονται στο πακέτο. Έτσι ι τα φορτώνονται από το /data
με την χρήση data()
data("nama_10r_3gdp") dat<-nama_10r_3gdp
α περιορίσουμε την εργασία εδώ για λόγους οικονομίας στον νομό Αργολίδας (και "αναγκαστικά" και Αργολίδας)
GDP651_2006<-dat[grep('^EL651',dat$geo),] kable(label_eurostat(GDP651_2006))
Για να γίνει γεωναφορά των δεδομένων αυτών θα χρησιμοποιήσουμε τα όρια όπως αυτά παρέχονται EUROSTAT.
par("mar"=c(.1,.1,.1,.1)) NUTS_2013_01M_651<-EtrsTransform(NUTS_2013_01M_EL[grep('^EL651',NUTS_2013_01M_EL[["NUTS_ID"]]) ,]) kable(NUTS_2013_01M_651@data)
Αρχικά γίνεται η γεωαναφορά των στοιχείων
NUTS_2013_01M_651<-merge(NUTS_2013_01M_651,GDP651_2006[1,],by.x="NUTS_ID",by.y="geo",all=FALSE)
H χρήση της merge ενδεχομένως να προκαλέσει αλλαγή στα rownames του dataframe οπότε έστω και προληπτικά τα διορθώνουμε
row.names(NUTS_2013_01M_651@data)<-sapply(slot(NUTS_2013_01M_651, "polygons"), function(x) slot(x, "ID")) plot(NUTS_2013_01M_651) kable(NUTS_2013_01M_651@data)
Στην συνέχεια θα χρησιμοποιήσουμε την etrsSourceSurface για να προβάλλουμε στον καναβο το "φαινόμενο" της κατανομής του κατά κεφαλήν ακαθάριστου εθνικού προϊόντος
par("mar"=c(.1,.1,.1,.1)) NUTS_2013_01M_651_GDP<-etrsSourceSurface(input.surface = NUTS_2013_01M_651,over.method.type = "PropCal",surface.value.col = 7,cell.size = 1000) plot(NUTS_2013_01M_651_GDP)
kable(head(NUTS_2013_01M_651_GDP@data))
Η προηγούμενη επιφάνεια dasymetric.surface που υπολογίστηκε ως Δασυμετρική επιφάνεια τώρα θα μετατραπεί σε βοηθητική επιφάνεια και θα χρησιμοποιεί για την κατανομή της τιμής της GDP. Σε αυτήν απεικονιζόταν ο αριθμός των κατοίκων ανά κελί που είχε την ιδιότητα του συνεχούς ή αστικού ιστού
kable(head(dasymetric.surface@data)) POP_2001_ancillary<- etrsDasymetric2Ancillary(dasymetric.surface) row.names(POP_2001_ancillary@data)<-sapply(slot(POP_2001_ancillary, "polygons"), function(x) slot(x, "ID"))
Τέλος με την κλίση της etrsPropWeightedValue θα υπολογιστεί η τιμή του ακαθάριστου προϊόντος σε κάθε κελί
par("mar"=c(0.1,0.1,0.1,0.1)) DASY_GPD<-etrsPropWeightedValue(input.surface.grided = NUTS_2013_01M_651_GDP, ancillary.grided = POP_2001_ancillary) plot(DASY_GPD) plot(ΑΡΓΟΛΙΔΑ.ETRS,add=T,border=2,lwd=2) kable(head(DASY_GPD))
Τέλος η μετατροπή των αποτελέσματα σε raster θα έκανε τα ευκολότερα διαθέσιμα για χρήση. Θα χρησιμοποιήσουμε την etrsDasymetric2Rasterπου εξαρτάται απο το πακέτο raster για να επιτύχουμε το ζητούμενο αποτέλεσμα
par("mar"=c(0.1,0.1,0.1,0.1)) DASY_GPD_RASTER<-etrsDasymetric2Raster(dasymetric.surface = DASY_GPD) rw.colors<-grey.colors image(DASY_GPD_RASTER,col=rw.colors(5)) plot(ΑΡΓΟΛΙΔΑ.ETRS,add=T,border=2,lwd=2)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.