DasyMapR - Ένας οδηγός για την χρήση του πακέτου

Το πακέτο DasyMapR

Το πακέτο της R που αναπτύχθηκε περιέχει μία σειρά εργαλείων με σκοπό την υλοποίηση δασυμετρικών χαρτών με την μέθοδο του επιμερισμού των δεδομένων στον κανονικό ETRS-LAEA κάναβο. Πληροφορίες για τον γεωγραφικό κανάβο μπορούν να βρεθούν INSPIRE Specification on Geographical Grid Systems

Τι περιέχεται σε αυτό το κείμενο

• Εγκατάσταση του πακέτου
• Τι υπάρχει στο πακέτο
• Φόρτωση Δεδομένων
• Απόδοση τιμής στο κελί - Κατηγορικά(ποιοτικά) Δεδομένα - MaxArea
• Απόδοση τιμής στο κελί - Αριθμητικά (ποσοτικά) Δεδομένα - PropCal
• Επιφάνεια Πηγή
• Βοηθητική Επιφάνεια
• Οι δασυμετρικοί υπολογισμοί
• Η Δασυμετρική επιφάνεια σε Βοηθητική
• H Δασυμετρική Επιφάνεια σε Raster

Εγκατάσταση

Το πακέτο δεν φιλοξενείται στο 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))

DasyMapR.contains

ActuallVal2Density
candidates_addresses
candidates_addresses_el_cor
CheckEtrsResolution
CheckEtrsValidity
CLC2000.ARGOLIDA.RES
CLC2000_CODES
CLC2000_GEO_CODES
clc_v2_code_11x_el
CORINE2000_CODES
CreateEtrsVectorGrid
DASY_GPD
dasymapPlot
EAP.application.loc
EAP.application.NUTS
EAP.application.pre_err
EAP_NUTS_V9
ETRS
etrsAncillarySurface
EtrsAncillarySurface
etrsAncillarySurface.default
etrsCellCenter
EtrsCellCodes
etrsCells
EtrsCells
EtrsCheckCodeColumns
etrsDasymetric2Ancillary
etrsDasymetric2Raster
etrsDasymetric2Source
EtrsDasymetricSurface
etrsGrid
etrsGrid2Spdf
etrsGrid.default
etrsMaxArea
etrsPoint2Grid
etrsPoints
EtrsPoints
etrsPropValue
etrsPropWeightedValue
etrsReverseCellCode
etrsSourceSurface
EtrsSourceSurface
etrsSurface
EtrsSurface
etrsSurface2Spdf
etrsSurface.default
etrsSurfacePar
EtrsTableCodes
EtrsTransform
GeoPc_GR_Places_etrs_k
GeoPc_GR_Places_etrs_k_el
GeoPc_GR_Places_etrs_k_en
GEOSTAT_grid_EU_POP_2006_1k_V1_1_1 GEOSTAT_grid_POP_1K_2011_V2_0_el
joinMaxAreaSurfaceDataFrames
mosaic_100m_sealing_el
nama_10r_3gdp
nama_10r_3popgdp
NUTS_2013_01M_EL
NUTS3_OCMG
NUTSV9_LEAC
NUTSV9_LEAC_GR_1km_grided
oikismoi
oikismoi_ETRS
oria_ota_etrs
pc_el_NUTS2010
pc_el_NUTS2013
pc_regions
pntsattr2surface
raster2Ancillary
wgsTransform

Περισσότερες πληροφορίες για το πακέτο θα μπορούσε να βρει κάποιος στα 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(0.1, 0.1, 0.1, 0.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)
# Καλεί την 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² και εκφράζει πυκνότητα. Το @data$VALUE (slot) πλέον περιέχει και την τιμή του εμβαδού σε km².

Θα επιμερίσουμε την τιμή με την χρήση της μεθόδου 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(0.1, 0.1, 0.1, 0.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(0.1, 0.1, 0.1, 0.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 = 0.5, border = "lightgrey")
plot(ΑΡΓΟΛΙΔΑ.ETRS, add = TRUE, lwd = 2, border = 2)

Η Βοηθιτική επιφάνεια

Θα γίνει χρήση των δεδομένων CORINE 2000 δηλαδή των δεδομένων κάλυψη γης για την Ελλάδα και το έτος 2000, όπως διατίθενται από Ευρωπαϊκή Υπηρεσία Περιβάλλοντος . Τα δεδομένα αυτά μπορούν να μεταμορφωθούν από τον ισότοπο της ΕΕΑ και επιλέχθηκαν ως κατάλληλα από την άποψη της χρονικής κατανομής σε σχέση με το φαινόμενο (πληθυσμό 2001) που θέλουμε να ανακατανείμουμε. Για την προετοιμασία της επιφάνειας επίσης θα γίνει η χρήση της etrsAncillarySurface. Εδώ κρίνεται σκόπιμο να γίνουν ορισμένες επισημάνσεις στο χρήστη του πακέτου. Παρόλο που οι τρεις μέθοδοι δημιουργίας των επιφανειών etrsSurface, etsrSourceSurface, etrsAncillarySurface επί της ουσίας κάνουν ακριβώς τους ίδιους υπολογισμούς, και οι κλάσεις περιέχουν τα ίδια @slots και θα μπορούσε ο χρήστης να αποφασίζει για την “ορθή” χρήση τους η επιλογή αυτή έγινε για να καθοδηγήσει τον χρήστη για το πώς θα πρέπει να ενεργήσει. Θα γίνουμε πιο συγκεκριμένη στην πορεία ολοκλήρωσης του τρέχοντος παραδείγματος.

data("CLC2000_CODES")
kable(head(CLC2000_CODES))

# Τα δύο .shp files περιέχονται στο folder corine
system.file("extdata", "clc00_v2_code_111", package = "DasyMapR")
system.file("extdata", "clc00_v2_code_112", 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)

# Φορτώνουμε τα δεδομένα
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))

par(mar = c(0.1, 0.1, 0.1, 0.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")

par(mar = c(0.1, 0.1, 0.1, 0.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")

par(mar = c(0.1, 0.1, 0.1, 0.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(0.1, 0.1, 0.1, 0.1))
dasymetric.surface <- EtrsDasymetricSurface(input.surface.grided = source.surface, 
    ancillary.grided = the.ancillary.surface.bf, actuall.value = FALSE)

## ~~~ ETRS validity ~~~

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(0.1, 0.1, 0.1, 0.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)

Η Δασυμετρική επιφάνεια σε Βοηθητική

Όπως αναφέρθηκε τα διοικητικά όρια συχνά αλλάζουν τέτοιο παράδειγμα είναι η περίπτωση του Ευρωπαϊκού χώρου όπου το ιστορικό των αλλαγων συχνά δυσκολεύει την χρήση και σύγκριση των διαθέσιμων δεδομένων. Τα γεωγραφικά όρια μπορούν να μεταμορφωθούν από το της EUROSTAT Ας δούμε κατ’ αρχήν τις αλλαγές των ορίων στην υπό μελέτη περιοχή την Αργολίδα :

Φορτώνουμε τα δεδομένα και τα :

par(mar = c(0.1, 0.1, 0.1, 0.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(0.1, 0.1, 0.1, 0.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)

par(mar = c(0.1, 0.1, 0.1, 0.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)

## ~~~ ETRS validity ~~~

plot(DASY_GPD)
plot(ΑΡΓΟΛΙΔΑ.ETRS, add = T, border = 2, lwd = 2)

kable(head(DASY_GPD))

H Δασυμετρική Επιφάνεια σε Raster

Τέλος η μετατροπή των αποτελέσματα σε 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)