radial KDE in R

radial-KDE.r

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

# get the full result set from the table
d <- dbGetQuery(connection,"SELECT t1, t2, ST_Azimuth(ST_StartPoint(line)::geography, ST_EndPoint(line)::geography) AS azimuth, ST_Distance(ST_StartPoint(line)::geography,ST_EndPoint(line)::geography)/1000 AS distance, white AS commuters FROM wbcommute_38300 WHERE t1 != t2 AND t1 = 684;")

dbDisconnect(connection)

# this is the number of commuter-kilometers
# density estimation weight field
d$weight = d$commuters * d$distance

# reverse direction
d$azimuth = abs(d$azimuth - 2*pi)
# move from 12oclock = 0 to 3oclock = 0
d$azimuth =  (d$azimuth+(1/2*pi)) %% (2*pi)

#calculate the density function
dens = density(
	d$azimuth,
	bw=0.1,
	weights=d$weight / sum(d$weight),
	n=1000
)
# subset back to the original 360 degrees
# and make a simpler data structure
a = dens$x		# x becomes angle
r = dens$y		# y becomes radius
table = data.frame(a,r)

# restrict the table to 360 degrees
range = table$a >= 0 & table$a <= 2*pi
table = table[range,]

# find points on the circle
table$x = 0.5 + table$r * cos(table$a)
table$y = 0.5 + table$r * sin(table$a)

# plot the circle
png('test.png',width=600,height=600)
	plot(
		1,
		type='n',
		xlim=c(0,1),
		ylim=c(0,1)
	)
	points(0.5,0.5)
	lines(table$x,table$y)
	# radians = 0
	lines(c(0.5,(0.5+1*cos(0))),c(0.5,(0.5+1*sin(0))),col='red')
	# radians = 1/2*pi
	lines(c(0.5,(0.5+1*cos(1/2*pi))),c(0.5,(0.5+1*sin(1/2*pi))),col='green')
dev.off()