population.R

## Synopsis
## ========
##
## Standalone R-Script to convert geostatistical population raster
## data (in CSV format) from Statistics Autria into MAT and CSV files.
##
## License
## =======
##
## Copyright (c) 2015 Markus Strauss
##
## Permission is hereby granted, free of charge, to any person obtaining a copy
## of this software and associated documentation files (the "Software"), to deal
## in the Software without restriction, including without limitation the rights
## to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
## copies of the Software, and to permit persons to whom the Software is
## furnished to do so, subject to the following conditions:
##
## The above copyright notice and this permission notice shall be included in all
## copies or substantial portions of the Software.
##
## THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
## IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
## FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
## AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
## LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
## OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
## SOFTWARE.
##
## Revisions
## =========
##
## 2015-10-06: Initial version.
## 2018-08-09: Standalone version with CSV export.
##

require(sp)
require(spatstat)
require(data.table)
require(mapdata)
require(maptools)
require(rmatio)

## define coordinate systems
epsg3035 <- CRS("+init=epsg:3035 +units=km")
wgs84    <- CRS("+proj=longlat +ellps=WGS84")

## load map of Austria
at <- map('worldHires','austria', fill=T, plot=F)
at_sp <- map2SpatialPolygons(at, IDs=at$names, proj4string = wgs84)
rm(at)

## convert map of Austria into epsg3035
at_sp_epsg3035 = spTransform( at_sp, epsg3035 )
at_window_epsg3035 = as.owin( at_sp_epsg3035 )

## load population data
pop <- fread( file.path( '.', 'GEOSTAT_grid_1K_pop_AT_2006.csv' ))
pat <- "1kmN([0-9]+)E([0-9]+)"

## the population data are coded such that the coordinates refer to
## the south-western corner of the corr. grid cell; pop_delta ensures
## the centering: set it to the grid width halved (for the 1km grid:
## 0.5)
pop_delta <- 0.5

pop[,y:=as.numeric(sub(pat, "\\1", pop$GRD_ID))]
pop[,x:=as.numeric(sub(pat, "\\2", pop$GRD_ID))]

pop[,GRD_ID:=NULL]
pop[,METHD_CL:=NULL]
pop[,YEAR:=NULL]

pop_sp <- pop

## change coordinates to refer to cell centers
pop_sp$x <- pop_sp$x + pop_delta
pop_sp$y <- pop_sp$y + pop_delta

coordinates(pop_sp) <- c('x','y')
proj4string(pop_sp) <- epsg3035

pop <- as.data.frame(pop_sp)[c('x','y','POP_TOT')]

pop_ppp <- ppp(pop_sp$x, pop_sp$y, window=at_window_epsg3035, marks=pop_sp$POP_TOT)
pop_ppp <- as.ppp( pop_ppp )  # remove rejected points

data <- list(pop=list(x=pop$x, y=pop$y, tot=pop$POP_TOT))

## export to CSV
write.csv(data, 'population-epsg3035.csv', row.names=F)

## export to MATLAB
write.mat(data,
          filename=file.path('.', 'population-epsg3035.mat'))