Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- 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()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement