ignore extdata folder

.gitignore

@@ -2,3 +2,4 @@
 .Rhistory
 .RData
 .Ruserdata
+inst/extdata

inst/COPYRIGHTS

+The R package and documentation itself is copyrighted by BNOSAC and the R package maintainer Jan Wijffels
+
+  (C) BNOSAC (www.bnosac.be/contact) 
+  
+Contains data of administrative identifiers and shapefiles provided by the Directie Statistiek - Statistics Belgium. 
+See how enrichment is done in inst/dev/01_import.R
+
+  (C) Algemene Directie Statistiek - Statistics Belgium, released under Licentie open data Algemene Directie Statistiek (see file LICENCE)
+
+

inst/LICENCE

+Licentie open data Algemene Directie Statistiek - Statistics Belgium
+Het hergebruik
+Je kan de “informatie” of “overheidsgegevens” die door de “producent” binnen de vrijheden en voorwaarden voorzien in onderhavige licentie hergebruiken.
+De “producent” garandeert de “hergebruiker” het persoonlijke, niet-exclusieve en kosteloze recht op hergebruik van informatie onderworpen aan onderhavige licentie, wereldwijd en voor onbeperkte duur, binnen de hieronder vermelde vrijheden en voorwaarden.
+- de “informatie” reproduceren, kopiëren, publiceren en doorgeven;
+- de “informatie” verspreiden en herverdelen;
+- de “informatie” gebruiken om info aan te passen, te wijzigen, te extraheren en om te vormen, in het bijzonder om “gegevens af te leiden;”
+- de “informatie” commercieel benutten, bijvoorbeeld door ze met andere “gegevens” te combineren, of door ze in je eigen product of applicatie te gebruiken.
+Onder de volgende voorwaarden:
+- De verwijzingen naar de “informatie” vermelden: de bron (ten minste de naam van de “producent”) en de datum van de laatste bijwerking.
+- De “hergebruiker” kan onder andere aan deze voorwaarde voldoen door een of meerdere hyperlinks (URL) aan te geven die verwijzen naar de “informatie” en die een effectieve vermelding van het auteurschap garanderen.
+- Het vermelden van de verwijzingen dient geen officieel karakter aan het hergebruik van de “informatie” te geven, noch enige erkenning of instemming van de “producent” of van enige ander publieke entiteit, van de “hergebruiker” of van het hergebruik ervan. Algemene Directie Statistiek - Statistics Belgium 22 oktober 2015
+Koning Albert II laan 16 - 1000 Brussel
+E-mail: mailto:statbel.opendata@economie.fgov.be
+Verantwoordelijkheid
+De “informatie” wordt ter beschikking gesteld zoals ze door de “producent” wordt aangemaakt of ontvangen, zonder enige andere niet door onderhavige licentie voorziene uitdrukkelijke of stilzwijgende garantie.
+De “producent”:
+- garandeert de “informatie” kosteloos ter beschikking te stellen binnen de vrijheden en voorwaarden vastgelegd in onderhavige licentie.
+- garandeert niet dat de “informatie” geen gebreken of onregelmatigheden bevat.
+- garandeert niet dat de “informatie” non-stop kan worden verschaft.
+- kan niet verantwoordelijk worden gesteld voor enig verlies, nadeel of schade aan derden door het hergebruik.
+De “hergebruiker” is als enige verantwoordelijk voor het hergebruik van de “informatie”, dat een legitiem doel moet hebben. Het hergebruik mag derden niet misleiden wat de inhoud van de “informatie,” de bron en de datum van bijwerking betreft.
+Intellectuele eigendomsrechten
+De “producent” garandeert dat de “informatie” geen intellectuele eigendomsrechten bevat die aan derden toekomen.
+De eventuele “Intellectuele eigendomsrechten” van de “producent” op documenten met “informatie” staan het vrije hergebruik van de “informatie” niet in de weg.
+Wanneer de “producent” de “Intellectuele eigendomsrechten” bezit op documenten die “informatie” bevatten, staat hij het hergebruik van deze “informatie” niet-exclusief, kosteloos en wereldwijd toe aan de “hergebruiker”, die er volledig van gebruik kan maken overeenkomstig de vrijheden en voorwaarden vastgelegd in onderhavige licentie.
+De “hergebruiker” behoudt alle intellectuele eigendomsrechten op het product dat hij maakt op basis van de geproduceerde of door de producent ter beschikking gestelde informatie, overeenkomstig de vrijheden en voorwaarden voorzien in onderhavige licentie.
+De licentie geldt voor een onbepaalde duur en neemt automatisch een einde ingeval van niet-naleving van de in onderhavige licentie vastgelegde voorwaarden. Algemene Directie Statistiek - Statistics Belgium 22 oktober 2015
+Koning Albert II laan 16 - 1000 Brussel
+E-mail: mailto:statbel.opendata@economie.fgov.be
+De compatibiliteit van de Licentie
+Om het hergebruik van de “gegevens” te vergemakkelijken, werd deze licentie ontwikkeld om compatibel te zijn met elke licentie die ten minste het vermelden van het auteurschap vereist. Ze is onder andere compatibel met de licenties:
+- “Open Government Licence” (OGL) van het Verenigd Koninkrijk
+- “Creative Commons Attribution 2.0” (CC-BY 2.0) van Creative Commons
+- “Open Data Commons Attribution” (ODC-BY) van de Open Knowledge Foundation
+Het toepasbare recht
+Onderhavige licentie wordt door het Belgisch recht geregeld.
+Er kan een beroep worden ingesteld bij de Commissie voor de toegang tot en het hergebruik van bestuursdocumenten bij niet-naleving van een of meerdere clausules van de licentie. Algemene Directie Statistiek - Statistics Belgium 22 oktober 2015
+Koning Albert II laan 16 - 1000 Brussel
+E-mail: mailto:statbel.opendata@economie.fgov.be
\ No newline at end of file

inst/NEWS

+                CHANGES IN BelgiumMaps.StatBel VERSION 1.0
+
+    o    Initial version of the package
+

inst/dev/01_import.R

+library(rgdal)
+library(sp)
+library(rgeos)
+library(tmap)
+library(utils)
+library(data.table)
+
+settings <- list()
+settings$hierarchy <- list()
+settings$hierarchy$municipalities <- file.path(getwd(), "inst", "extdata", "TF_PSNL_INC_TAX_MUNTY_tcm325-272513.zip")
+settings$hierarchy$nis.sectors <- file.path(getwd(), "inst", "extdata", "TF_PSNL_INC_TAX_SECTOR_tcm325-278417.zip")
+settings$nis.sectors <- file.path(getwd(), "inst", "extdata", "scbel01012011_gen13_tcm325-275679", "scbel01012011_gen13.shp")
+settings$agglomeraties <- file.path(getwd(), "inst", "extdata", "SH_CENSUS_2011_AGGLOMERATIONS_200M.shp_tcm325-276234")
+
+
+
+##
+## Import NIS hierarchy
+##
+unzip(settings$hierarchy$municipalities, list = TRUE)
+unzip(settings$hierarchy$nis.sectors, list = TRUE)
+hierarchy <- list()
+hierarchy$municipalities <- read.table(unz(settings$hierarchy$municipalities, filename = "TF_PSNL_INC_TAX_MUNTY.txt"), sep="|", header = TRUE, 
+                                       encoding = "UTF-8", stringsAsFactors = FALSE, quote = "")
+hierarchy$nis.sectors <- read.table(unz(settings$hierarchy$nis.sectors, filename = "TF_PSNL_INC_TAX_SECTOR.txt"), sep="|", header = TRUE, 
+                                       encoding = "UTF-8", stringsAsFactors = FALSE, quote = "")
+hierarchy$municipalities <- unique(hierarchy$municipalities[, c("CD_MUNTY_REFNIS", "TX_MUNTY_DESCR_NL", "TX_MUNTY_DESCR_FR", 
+                                                         "CD_DSTR_REFNIS", "TX_ADM_DSTR_DESCR_NL", "TX_ADM_DSTR_DESCR_FR", 
+                                                         "CD_PROV_REFNIS", "TX_PROV_DESCR_NL", "TX_PROV_DESCR_FR", 
+                                                         "CD_RGN_REFNIS", "TX_RGN_DESCR_NL", "TX_RGN_DESCR_FR")])
+hierarchy$nis.sectors <- unique(hierarchy$nis.sectors[, c("CD_MUNTY_REFNIS", "CD_SECTOR", "CD_REFNIS_SECTOR", "TX_SECTOR_DESCR_NL", "TX_SECTOR_DESCR_FR")])
+sum(duplicated(hierarchy$municipalities$CD_MUNTY_REFNIS)) ## doublecheck to make sure names did not change over the years
+sum(duplicated(hierarchy$nis.sectors$CD_REFNIS_SECTOR)) ## doublecheck to make sure names did not change over the years
+hierarchy <- merge(hierarchy$nis.sectors, hierarchy$municipalities, by = "CD_MUNTY_REFNIS")
+BE_ADMIN_HIERARCHY <- hierarchy[, c("CD_REFNIS_SECTOR", "CD_SECTOR", "TX_SECTOR_DESCR_NL", "TX_SECTOR_DESCR_FR", 
+                                    "CD_MUNTY_REFNIS", "TX_MUNTY_DESCR_NL", "TX_MUNTY_DESCR_FR", 
+                                    "CD_DSTR_REFNIS", "TX_ADM_DSTR_DESCR_NL", "TX_ADM_DSTR_DESCR_FR", 
+                                    "CD_PROV_REFNIS", "TX_PROV_DESCR_NL", "TX_PROV_DESCR_FR", 
+                                    "CD_RGN_REFNIS", "TX_RGN_DESCR_NL", "TX_RGN_DESCR_FR")]
+
+##
+## Import shape of agglomerations
+##
+ogrListLayers(settings$agglomeraties)
+ogrInfo(settings$agglomeraties, layer = "SH_CENSUS_2011_AGGLOMERATIONS_200M")
+
+BE_ADMIN_AGGLOMERATIONS <- readOGR(settings$agglomeraties, layer = "SH_CENSUS_2011_AGGLOMERATIONS_200M", stringsAsFactors = FALSE)
+colnames(BE_ADMIN_AGGLOMERATIONS@data) <- tolower(make.names(colnames(BE_ADMIN_AGGLOMERATIONS@data), allow_=FALSE))
+BE_ADMIN_AGGLOMERATIONS <- spTransform(x=BE_ADMIN_AGGLOMERATIONS, CRSobj = CRS("+proj=longlat +datum=WGS84"))
+
+##
+## Import shape of NIS sectors
+##
+ogrListLayers(settings$nis.sectors)
+ogrInfo(settings$nis.sectors, layer = "scbel01012011_gen13")
+
+
+BE_ADMIN_SECTORS <- readOGR(settings$nis.sectors, layer = "scbel01012011_gen13", stringsAsFactors = FALSE)
+colnames(BE_ADMIN_SECTORS@data) <- tolower(make.names(colnames(BE_ADMIN_SECTORS@data), allow_=FALSE))
+
+isutf8 <- function(x){
+  for(column in colnames(x)){
+    if(is.character(x[[column]])){
+      Encoding(x[[column]]) <- "UTF-8"
+    }
+  }
+  x
+}
+BE_ADMIN_SECTORS@data <- isutf8(BE_ADMIN_SECTORS@data)
+BE_ADMIN_SECTORS$nuts0 <- "BE"
+BE_ADMIN_SECTORS@data <- setnames(BE_ADMIN_SECTORS@data, old = "nuts3.new", new = "nuts3")
+BE_ADMIN_SECTORS@data <- BE_ADMIN_SECTORS@data[, c(setdiff(names(BE_ADMIN_SECTORS), c("nuts0", "nuts1", "nuts2", "nuts3")), c("nuts0", "nuts1", "nuts2", "nuts3"))]
+## make brussels into a province
+BE_ADMIN_SECTORS$prov.nl <- ifelse(is.na(BE_ADMIN_SECTORS$prov.nl), BE_ADMIN_SECTORS$reg.nl, BE_ADMIN_SECTORS$prov.nl)
+BE_ADMIN_SECTORS$prov.fr <- ifelse(is.na(BE_ADMIN_SECTORS$prov.fr), BE_ADMIN_SECTORS$reg.fr, BE_ADMIN_SECTORS$prov.fr)
+## add area in km2
+BE_ADMIN_SECTORS$gis.area.km2 <- BE_ADMIN_SECTORS$gis.area.h / 100
+BE_ADMIN_SECTORS$cad.area.km2 <- BE_ADMIN_SECTORS$cad.area.h / 100
+## convert to WGS84 coordinates instead of Lambert
+BE_ADMIN_SECTORS <- spTransform(x=BE_ADMIN_SECTORS, CRSobj = CRS("+proj=longlat +datum=WGS84"))
+
+## remove redundant data (not needed as in sector)
+BE_ADMIN_SECTORS$nis.012011 <- NULL
+BE_ADMIN_SECTORS$sec012011 <- NULL
+
+## add information of hierarchy based on NIS definition in 20/01/2011
+## mark that BE_ADMIN_HIERARCHY contains sectors not in BE_ADMIN_SECTORS$cs012011 (TAXES not linked to sector or to sector in other period namely cs102001)
+##   this can be checked with x <- subset(BE_ADMIN_HIERARCHY, !CD_REFNIS_SECTOR %in%  BE_ADMIN_SECTORS$cs012011)
+## also sectors seem to exist in the BE_ADMIN_SECTORS which are not in the taxes data first fix these
+x <- subset(BE_ADMIN_SECTORS@data, !cs012011 %in% BE_ADMIN_HIERARCHY$CD_REFNIS_SECTOR)
+x$CD_REFNIS_SECTOR <- x$cs012011
+x$CD_MUNTY_REFNIS <- substr(x$CD_REFNIS_SECTOR, 1, 5)
+x$CD_SECTOR <- substr(x$CD_REFNIS_SECTOR, 6, nchar(x$CD_REFNIS_SECTOR))
+x <- setnames(x, old = c("sector.nl", "sector.fr"), new = c("TX_SECTOR_DESCR_NL", "TX_SECTOR_DESCR_FR"))
+x <- x[, c("CD_REFNIS_SECTOR", "CD_SECTOR", "TX_SECTOR_DESCR_NL", "TX_SECTOR_DESCR_FR", "CD_MUNTY_REFNIS")]
+x <- merge(x, unique(BE_ADMIN_HIERARCHY[, c(#"CD_REFNIS_SECTOR", "CD_SECTOR", "TX_SECTOR_DESCR_NL", "TX_SECTOR_DESCR_FR", 
+  "CD_MUNTY_REFNIS", "TX_MUNTY_DESCR_NL", "TX_MUNTY_DESCR_FR", 
+  "CD_DSTR_REFNIS", "TX_ADM_DSTR_DESCR_NL", "TX_ADM_DSTR_DESCR_FR", 
+  "CD_PROV_REFNIS", "TX_PROV_DESCR_NL", "TX_PROV_DESCR_FR", 
+  "CD_RGN_REFNIS", "TX_RGN_DESCR_NL", "TX_RGN_DESCR_FR")]), by = "CD_MUNTY_REFNIS")
+x <- x[, names(BE_ADMIN_HIERARCHY)]
+BE_ADMIN_HIERARCHY <- rbind(BE_ADMIN_HIERARCHY, x) 
+## mark that sectors can be split or aggregated from year to year (see subset(BE_ADMIN_SECTORS, cs102001 == "52063A372")@data) and Hmisc::describe(BE_ADMIN_SECTORS@data) (look at uniques in sector ids)
+## there seems to be one element per cs012011 but doubles in older NIS versions. will use version cs012011
+BE_ADMIN_SECTORS$CD_REFNIS_SECTOR <- BE_ADMIN_SECTORS$cs012011
+BE_ADMIN_SECTORS <- merge(BE_ADMIN_SECTORS, BE_ADMIN_HIERARCHY, by.x = "CD_REFNIS_SECTOR", by.y = "CD_REFNIS_SECTOR", all.x=FALSE, all.y=FALSE)
+BE_ADMIN_SECTORS@data <- BE_ADMIN_SECTORS@data[, c(setdiff(names(BE_ADMIN_SECTORS), colnames(BE_ADMIN_HIERARCHY)), colnames(BE_ADMIN_HIERARCHY))]
+
+## add nuts information to hierarchy
+BE_ADMIN_HIERARCHY <- merge(BE_ADMIN_HIERARCHY, 
+                            unique(na.exclude(BE_ADMIN_SECTORS@data[, c("CD_MUNTY_REFNIS", "nuts0", "nuts1", "nuts2", "nuts3")])), 
+                            by = "CD_MUNTY_REFNIS", all.x=TRUE, all.y=FALSE)
+BE_ADMIN_HIERARCHY <- BE_ADMIN_HIERARCHY[, c("CD_REFNIS_SECTOR", "CD_SECTOR", "TX_SECTOR_DESCR_NL", "TX_SECTOR_DESCR_FR", 
+                                             "CD_MUNTY_REFNIS", "TX_MUNTY_DESCR_NL", "TX_MUNTY_DESCR_FR", 
+                                             "CD_DSTR_REFNIS", "TX_ADM_DSTR_DESCR_NL", "TX_ADM_DSTR_DESCR_FR", 
+                                             "CD_PROV_REFNIS", "TX_PROV_DESCR_NL", "TX_PROV_DESCR_FR", 
+                                             "CD_RGN_REFNIS", "TX_RGN_DESCR_NL", "TX_RGN_DESCR_FR", "nuts0", "nuts1", "nuts2", "nuts3")]
+
+x <- split(BE_ADMIN_SECTORS, BE_ADMIN_SECTORS$CD_MUNTY_REFNIS)
+x <- lapply(x, FUN=function(data){
+  SpatialPolygonsDataFrame(gUnaryUnion(spgeom = data),
+                           data = data.frame(CD_MUNTY_REFNIS = unique(data$CD_MUNTY_REFNIS),
+                                             TX_MUNTY_DESCR_NL = unique(data$TX_MUNTY_DESCR_NL),
+                                             TX_MUNTY_DESCR_FR = unique(data$TX_MUNTY_DESCR_FR),
+                                             CD_DSTR_REFNIS = unique(data$CD_DSTR_REFNIS),
+                                             TX_ADM_DSTR_DESCR_NL = unique(data$TX_ADM_DSTR_DESCR_NL),
+                                             TX_ADM_DSTR_DESCR_FR = unique(data$TX_ADM_DSTR_DESCR_FR),
+                                             CD_PROV_REFNIS = unique(data$CD_PROV_REFNIS),
+                                             TX_PROV_DESCR_NL = unique(data$TX_PROV_DESCR_NL),
+                                             TX_PROV_DESCR_FR = unique(data$TX_PROV_DESCR_FR),
+                                             CD_RGN_REFNIS = unique(data$CD_RGN_REFNIS),
+                                             TX_RGN_DESCR_NL = unique(data$TX_RGN_DESCR_NL),
+                                             TX_RGN_DESCR_FR = unique(data$TX_RGN_DESCR_FR),
+                                             nuts0 = paste(unique(data$nuts0), collapse=","),
+                                             nuts1 = paste(unique(data$nuts1), collapse=","),
+                                             nuts2 = paste(unique(data$nuts2), collapse=","),
+                                             nuts3 = paste(unique(data$nuts3), collapse=","),
+                                             SURFACE.GIS.h = sum(data@data$gis.area.h),
+                                             SURFACE.CAD.h = sum(data@data$cad.area.h),
+                                             SURFACE.GIS.km2 = sum(data@data$gis.area.km2),
+                                             SURFACE.CAD.km2 = sum(data@data$cad.area.km2),
+                                             stringsAsFactors=FALSE))
+})
+BE_ADMIN_MUNTY <- do.call(sbind, x)
+plot(BE_ADMIN_MUNTY)
+
+x <- split(BE_ADMIN_SECTORS, BE_ADMIN_SECTORS$CD_DSTR_REFNIS)
+x <- lapply(x, FUN=function(data){
+  SpatialPolygonsDataFrame(gUnaryUnion(spgeom = data),
+                           data = data.frame(CD_DSTR_REFNIS = unique(data$CD_DSTR_REFNIS),
+                                             TX_ADM_DSTR_DESCR_NL = unique(data$TX_ADM_DSTR_DESCR_NL),
+                                             TX_ADM_DSTR_DESCR_FR = unique(data$TX_ADM_DSTR_DESCR_FR),
+                                             CD_PROV_REFNIS = unique(data$CD_PROV_REFNIS),
+                                             TX_PROV_DESCR_NL = unique(data$TX_PROV_DESCR_NL),
+                                             TX_PROV_DESCR_FR = unique(data$TX_PROV_DESCR_FR),
+                                             CD_RGN_REFNIS = unique(data$CD_RGN_REFNIS),
+                                             TX_RGN_DESCR_NL = unique(data$TX_RGN_DESCR_NL),
+                                             TX_RGN_DESCR_FR = unique(data$TX_RGN_DESCR_FR),
+                                             nuts0 = paste(unique(data$nuts0), collapse=","),
+                                             nuts1 = paste(unique(data$nuts1), collapse=","),
+                                             nuts2 = paste(unique(data$nuts2), collapse=","),
+                                             nuts3 = paste(unique(data$nuts3), collapse=","),
+                                             SURFACE.GIS.h = sum(data@data$gis.area.h),
+                                             SURFACE.CAD.h = sum(data@data$cad.area.h),
+                                             SURFACE.GIS.km2 = sum(data@data$gis.area.km2),
+                                             SURFACE.CAD.km2 = sum(data@data$cad.area.km2),
+                                             stringsAsFactors=FALSE))
+})
+BE_ADMIN_DISTRICT <- do.call(sbind, x)
+plot(BE_ADMIN_DISTRICT)
+
+x <- split(BE_ADMIN_SECTORS, BE_ADMIN_SECTORS$CD_PROV_REFNIS)
+x <- lapply(x, FUN=function(data){
+  SpatialPolygonsDataFrame(gUnaryUnion(spgeom = data),
+                           data = data.frame(CD_PROV_REFNIS = unique(data$CD_PROV_REFNIS),
+                                             TX_PROV_DESCR_NL = unique(data$TX_PROV_DESCR_NL),
+                                             TX_PROV_DESCR_FR = unique(data$TX_PROV_DESCR_FR),
+                                             CD_RGN_REFNIS = unique(data$CD_RGN_REFNIS),
+                                             TX_RGN_DESCR_NL = unique(data$TX_RGN_DESCR_NL),
+                                             TX_RGN_DESCR_FR = unique(data$TX_RGN_DESCR_FR),
+                                             nuts0 = paste(unique(data$nuts0), collapse=","),
+                                             nuts1 = paste(unique(data$nuts1), collapse=","),
+                                             nuts2 = paste(unique(data$nuts2), collapse=","),
+                                             SURFACE.GIS.h = sum(data@data$gis.area.h),
+                                             SURFACE.CAD.h = sum(data@data$cad.area.h),
+                                             SURFACE.GIS.km2 = sum(data@data$gis.area.km2),
+                                             SURFACE.CAD.km2 = sum(data@data$cad.area.km2),
+                                             stringsAsFactors=FALSE))
+})
+BE_ADMIN_PROVINCE <- do.call(sbind, x)
+plot(BE_ADMIN_PROVINCE)
+
+
+x <- split(BE_ADMIN_SECTORS, BE_ADMIN_SECTORS$CD_RGN_REFNIS)
+x <- lapply(x, FUN=function(data){
+  SpatialPolygonsDataFrame(gUnaryUnion(spgeom = data),
+                           data = data.frame(CD_RGN_REFNIS = unique(data$CD_RGN_REFNIS),
+                                             TX_RGN_DESCR_NL = unique(data$TX_RGN_DESCR_NL),
+                                             TX_RGN_DESCR_FR = unique(data$TX_RGN_DESCR_FR),
+                                             nuts0 = paste(unique(data$nuts0), collapse=","),
+                                             nuts1 = paste(unique(data$nuts1), collapse=","),
+                                             SURFACE.GIS.h = sum(data@data$gis.area.h),
+                                             SURFACE.CAD.h = sum(data@data$cad.area.h),
+                                             SURFACE.GIS.km2 = sum(data@data$gis.area.km2),
+                                             SURFACE.CAD.km2 = sum(data@data$cad.area.km2),
+                                             stringsAsFactors=FALSE))
+})
+BE_ADMIN_REGION <- do.call(sbind, x)
+plot(BE_ADMIN_REGION)
+
+BE_ADMIN_BELGIUM <- SpatialPolygonsDataFrame(gUnaryUnion(spgeom = BE_ADMIN_SECTORS),
+                                             data = data.frame(nuts0 = paste(unique(BE_ADMIN_SECTORS$nuts0), collapse=","),
+                                                               SURFACE.GIS.h = sum(BE_ADMIN_SECTORS@data$gis.area.h),
+                                                               SURFACE.CAD.h = sum(BE_ADMIN_SECTORS@data$cad.area.h),
+                                                               SURFACE.GIS.km2 = sum(BE_ADMIN_SECTORS@data$gis.area.km2),
+                                                               SURFACE.CAD.km2 = sum(BE_ADMIN_SECTORS@data$cad.area.km2)))
+
+
+save(BE_ADMIN_HIERARCHY, file = file.path(getwd(), "data", "BE_ADMIN_HIERARCHY.RData"), compress = 'xz')
+save(BE_ADMIN_SECTORS, file = file.path(getwd(), "data", "BE_ADMIN_SECTORS.RData"), compress = 'xz')
+save(BE_ADMIN_MUNTY, file = file.path(getwd(), "data", "BE_ADMIN_MUNTY.RData"), compress = 'xz')
+save(BE_ADMIN_DISTRICT, file = file.path(getwd(), "data", "BE_ADMIN_DISTRICT.RData"), compress = 'xz')
+save(BE_ADMIN_PROVINCE, file = file.path(getwd(), "data", "BE_ADMIN_PROVINCE.RData"), compress = 'xz')
+save(BE_ADMIN_REGION, file = file.path(getwd(), "data", "BE_ADMIN_REGION.RData"), compress = 'xz')
+save(BE_ADMIN_BELGIUM, file = file.path(getwd(), "data", "BE_ADMIN_BELGIUM.RData"), compress = 'xz')
+save(BE_ADMIN_AGGLOMERATIONS, file = file.path(getwd(), "data", "BE_ADMIN_AGGLOMERATIONS.RData"), compress = 'xz')
+
+
+# library(rmapshaper)
+# x <- ms_simplify(BE_ADMIN_REGION, keep = 0.8)
+# plot(BE_ADMIN_BELGIUM)
+# save(BE_ADMIN_SECTORS, file = "BE_ADMIN_SECTORS.RData", compress = 'xz')
+# save(BE_ADMIN_MUNTY, file = "BE_ADMIN_MUNTY.RData", compress = 'xz')
+# save(BE_ADMIN_DISTRICT, file = "BE_ADMIN_DISTRICT.RData", compress = 'xz')
+# save(BE_ADMIN_PROVINCE, file = "BE_ADMIN_PROVINCE.RData", compress = 'xz')
+# save(BE_ADMIN_REGION, file = "BE_ADMIN_REGION.RData", compress = 'xz')
+# save(BE_ADMIN_BELGIUM, file = "BE_ADMIN_BELGIUM.RData", compress = 'xz')
+# plot(x)
+# x <- gSimplify(BE_ADMIN_BELGIUM, tol = 10)
+# 
+# Hmisc::describe(BE_ADMIN_SECTORS@data)

inst/dev/09_example_visualisation.R

+library(sp)
+library(leaflet)
+library(BelgiumMaps.StatBel)
+
+## Get taxes / statistical sector
+tempfile <- tempfile()
+download.file("http://statbel.fgov.be/nl/binaries/TF_PSNL_INC_TAX_SECTOR_tcm325-278417.zip", tempfile)
+unzip(tempfile, list = TRUE)
+taxes <- read.table(unz(tempfile, filename = "TF_PSNL_INC_TAX_SECTOR.txt"), sep="|", header = TRUE, 
+                    encoding = "UTF-8", stringsAsFactors = FALSE, quote = "", na.strings = c("", "C"))
+colnames(taxes)[1] <- "CD_YEAR"
+
+## Get taxes in last year
+taxes <- subset(taxes, CD_YEAR == max(taxes$CD_YEAR))
+taxes <- taxes[, c("CD_YEAR", "CD_REFNIS_SECTOR", 
+                   "MS_NBR_NON_ZERO_INC", "MS_TOT_NET_TAXABLE_INC", "MS_AVG_TOT_NET_TAXABLE_INC", 
+                   "MS_MEDIAN_NET_TAXABLE_INC", "MS_INT_QUART_DIFF", "MS_INT_QUART_COEFF", 
+                   "MS_INT_QUART_ASSYM")]
+
+## Join taxes with the map
+data(BE_ADMIN_SECTORS, package = "BelgiumMaps.StatBel")
+data(BE_ADMIN_DISTRICT, package = "BelgiumMaps.StatBel")
+data(BE_ADMIN_MUNTY, package = "BelgiumMaps.StatBel")
+
+mymap <- merge(BE_ADMIN_SECTORS, taxes, by = "CD_REFNIS_SECTOR", all.x=TRUE, all.y=FALSE)
+
+## Map the data
+mymap <- subset(mymap, TX_RGN_DESCR_NL %in% "Brussels Hoofdstedelijk Gewest")
+pal <- colorBin(palette = rev(heat.colors(11)), domain = mymap$MS_AVG_TOT_NET_TAXABLE_INC, 
+                bins = c(0, round(quantile(mymap$MS_AVG_TOT_NET_TAXABLE_INC, na.rm=TRUE, probs = seq(0.1, 0.9, by = 0.1)), 0), +Inf),
+                na.color = "#cecece")
+m <- leaflet(mymap) %>%
+  addTiles() %>%
+  addLegend(title = "Net Taxable Income (EURO)", pal = pal, values = ~MS_AVG_TOT_NET_TAXABLE_INC, position = "bottomleft", na.label = "Missing") %>%
+  addPolygons(color = ~pal(MS_AVG_TOT_NET_TAXABLE_INC), stroke = FALSE, smoothFactor = 0.2, fillOpacity = 0.85, 
+              popup = sprintf("%s: %s<br>%s: %s<br><br>%s €: Average net taxable income<br>%s €: Median net taxable income<br>%s declarations", 
+                              mymap$TX_SECTOR_DESCR_NL, mymap$TX_MUNTY_DESCR_NL, 
+                              mymap$TX_SECTOR_DESCR_FR, mymap$TX_MUNTY_DESCR_FR,  
+                              mymap$MS_AVG_TOT_NET_TAXABLE_INC, mymap$MS_MEDIAN_NET_TAXABLE_INC, 
+                              mymap$MS_NBR_NON_ZERO_INC))
+m <- addPolylines(m, data = BE_ADMIN_DISTRICT, weight = 1.5, color = "black")
+m
+
+library(htmlwidgets)
+m$sizingPolicy$defaultHeight <- "100%"
+m$sizingPolicy$defaultWidth <- "100%"
+saveWidget(m, file="be_net_taxable_income_sector.html")

inst/extdata/SH_CENSUS_2011_AGGLOMERATIONS_200M.shp_tcm325-276234/SH_CENSUS_2011_AGGLOMERATIONS_200M.cpg

+UTF-8
\ No newline at end of file

inst/extdata/SH_CENSUS_2011_AGGLOMERATIONS_200M.shp_tcm325-276234/SH_CENSUS_2011_AGGLOMERATIONS_200M.prj

+PROJCS["Belge_1972_Belgian_Lambert_72",GEOGCS["GCS_Belge 1972",DATUM["D_Belge_1972",SPHEROID["International_1924",6378388,297]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["standard_parallel_1",51.16666723333333],PARAMETER["standard_parallel_2",49.8333339],PARAMETER["latitude_of_origin",90],PARAMETER["central_meridian",4.367486666666666],PARAMETER["false_easting",150000.013],PARAMETER["false_northing",5400088.438],UNIT["Meter",1]]
\ No newline at end of file

inst/extdata/SH_CENSUS_2011_AGGLOMERATIONS_200M.shp_tcm325-276234/SH_CENSUS_2011_AGGLOMERATIONS_200M.qpj

+PROJCS["Belge 1972 / Belgian Lambert 72",GEOGCS["Belge 1972",DATUM["Reseau_National_Belge_1972",SPHEROID["International 1924",6378388,297,AUTHORITY["EPSG","7022"]],TOWGS84[-106.869,52.2978,-103.724,0.3366,-0.457,1.8422,-1.2747],AUTHORITY["EPSG","6313"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4313"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["standard_parallel_1",51.16666723333333],PARAMETER["standard_parallel_2",49.8333339],PARAMETER["latitude_of_origin",90],PARAMETER["central_meridian",4.367486666666666],PARAMETER["false_easting",150000.013],PARAMETER["false_northing",5400088.438],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],AUTHORITY["EPSG","31370"]]

inst/extdata/scbel01012011_gen13_tcm325-275679/scbel01012011_gen13.CPG

+UTF-8
\ No newline at end of file

inst/extdata/scbel01012011_gen13_tcm325-275679/scbel01012011_gen13.prj

+PROJCS["Belge_Lambert_1972",GEOGCS["GCS_Belge_1972",DATUM["D_Belge_1972",SPHEROID["International_1924",6378388.0,297.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",150000.01256],PARAMETER["False_Northing",5400088.4378],PARAMETER["Central_Meridian",4.367486666666666],PARAMETER["Standard_Parallel_1",49.8333339],PARAMETER["Standard_Parallel_2",51.16666733333333],PARAMETER["Latitude_Of_Origin",90.0],UNIT["Meter",1.0]]
\ No newline at end of file

inst/extdata/scbel01012011_gen13_tcm325-275679/scbel01012011_gen13.shp.xml

+<?xml version="1.0" encoding="UTF-8"?>
+<metadata xml:lang="fr"><Esri><CreaDate>20151207</CreaDate><CreaTime>11395700</CreaTime><ArcGISFormat>1.0</ArcGISFormat><SyncOnce>TRUE</SyncOnce><DataProperties><lineage><Process ToolSource="c:\program files (x86)\arcgis\desktop10.2\ArcToolbox\Toolboxes\Data Management Tools.tbx\CalculateField" Date="20150306" Time="100642">CalculateField scbel01012011_gen10 Gis_area_h [Area_ha] VB #</Process><Process ToolSource="c:\program files (x86)\arcgis\desktop10.2\ArcToolbox\Toolboxes\Data Management Tools.tbx\CalculateField" Date="20150306" Time="101024">CalculateField scbel01012011_gen10 Cad_area_h [Cad_area_a]/100 VB #</Process></lineage><itemProps><itemLocation><linkage Sync="TRUE">file://\\mineco2000\dfs-eco\E8\0536-849-GIS\GIS\BASE\diffusisec2011\diffusionsectgeneralise\scbel01012011_gen13</linkage><protocol Sync="TRUE">Local Area Network</protocol></itemLocation></itemProps><copyHistory><copy source="Z:\E8\0536-849-GIS\GIS\BASE\e4.sec2011\scbel01012011_gen13" dest="\\mineco2000\dfs-eco\E8\0536-849-GIS\GIS\BASE\diffusisec2011\diffusionsectgeneralise\scbel01012011_gen13" date="20151207" time="11395700"></copy></copyHistory></DataProperties></Esri></metadata>