R smooth ratio raster map

kde-ratio-map.r

# make dead-ends ratio map

library('ks')
library('raster')

# connect to database
library("RPostgreSQL")
driver <- dbDriver("PostgreSQL")
connection <- dbConnect(driver, dbname="", password="", user="", host="localhost")

cities = read.csv('cities.csv',stringsAsFactors=FALSE)

# iterate over cities
for( i in 1:length(cities$name)){
	city = cities[i,'name']
	epsg = cities[i,'epsg']
	print(paste('working on',city))

	d = dbGetQuery(connection,paste(
		"SELECT",
			"km,flags,",
			"car_comp AS cc,",
			"bike_comp AS bc,",
			"foot_comp AS fc,",
			"car_directness AS cd,",
			"bike_directness AS bd,",
			"foot_directness AS fd,",
			"ST_X(ST_Transform(ST_Centroid(geom_way),",epsg,")) AS x,",	# meters
			"ST_Y(ST_Transform(ST_Centroid(geom_way),",epsg,")) AS y",	# meters
		"FROM",city
		))

	# calculate grid resolution
	cell_size = 250 # square cell size in METERS
	bandwidth = 2000^2 # meters^2
	xres = round((max(d$x)-min(d$x)) / cell_size)
	yres = round((max(d$y)-min(d$y)) / cell_size)

	# subset for each mode
	car = d[!is.na(d$cc),c('x','y','km','cc','cd')]
	#bike = d[!is.na(d$bc),c('x','y','km','bc','bd')]
	#foot = d[!is.na(d$fc),c('x','y','km','fc','fd')]

	# calculate the KDE
	total_kde = kde(
		x=matrix(c(car$x,car$y),length(car$x),2),
		w=car$km / sum(car$km) * length(car$km),
		H=matrix(c(bandwidth,0,0,bandwidth),2),
		xmin=c(min(d$x),min(d$y)),
		xmax=c(max(d$x),max(d$y)),
		gridsize=c(xres,yres),
		approx.cont=FALSE,
		verbose=TRUE
	)
	# rasterize
	total = raster(total_kde)
	projection(total) = CRS(paste0('+init=epsg:',epsg))
	# scale to the length of road
	total = (total / sum(total_kde$estimate)) * sum(car$km)

	# deadends
	deadcar = car[car$cc==-1,]
	dead_kde = kde(
		x=matrix(c(deadcar$x,deadcar$y),length(deadcar$x),2),
		w=deadcar$km / sum(deadcar$km) * length(deadcar$km),
		H=total_kde$H, # reuse H from the total version
		xmin=c(min(d$x),min(d$y)),
		xmax=c(max(d$x),max(d$y)),
		gridsize=c(xres,yres),
		verbose=TRUE
	)
	# rasterize
	dead = raster(dead_kde)
	projection(dead) = CRS(paste0('+init=epsg:',epsg))
	# scale to the length of dead-end road
	dead = (dead / sum(dead_kde$estimate)) * sum(deadcar$km)

	# Percent Dead-ends
	percent = dead / total
	percent[is.na(percent)] = 0
	percent[percent>1] = 1

	banded_raster = addLayer(total,dead,percent)
	names(banded_raster) = c('car_total','car_dead','car_percent')

	# output the file as geotiff
	writeRaster(
		banded_raster,
		filename=paste0('data/ratio-rasters/',city,'_car.tiff'),
		format='GTiff',
		overwrite=TRUE
	)
}

dbDisconnect(connection)