Advertisement
Guest User

Untitled

a guest
Aug 26th, 2016
80
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 1.53 KB | None | 0 0
  1. library("RPostgreSQL")
  2. driver <- dbDriver("PostgreSQL")
  3. connection <- dbConnect(driver, dbname="", password="", user="", host="localhost")
  4.  
  5. # get the full result set from the table
  6. 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;")
  7.  
  8. dbDisconnect(connection)
  9.  
  10. # this is the number of commuter-kilometers
  11. # density estimation weight field
  12. d$weight = d$commuters * d$distance
  13.  
  14. # reverse direction
  15. d$azimuth = abs(d$azimuth - 2*pi)
  16. # move from 12oclock = 0 to 3oclock = 0
  17. d$azimuth = (d$azimuth+(1/2*pi)) %% (2*pi)
  18.  
  19. #calculate the density function
  20. dens = density(
  21. d$azimuth,
  22. bw=0.1,
  23. weights=d$weight / sum(d$weight),
  24. n=1000
  25. )
  26. # subset back to the original 360 degrees
  27. # and make a simpler data structure
  28. a = dens$x # x becomes angle
  29. r = dens$y # y becomes radius
  30. table = data.frame(a,r)
  31.  
  32. # restrict the table to 360 degrees
  33. range = table$a >= 0 & table$a <= 2*pi
  34. table = table[range,]
  35.  
  36. # find points on the circle
  37. table$x = 0.5 + table$r * cos(table$a)
  38. table$y = 0.5 + table$r * sin(table$a)
  39.  
  40. # plot the circle
  41. png('test.png',width=600,height=600)
  42. plot(
  43. 1,
  44. type='n',
  45. xlim=c(0,1),
  46. ylim=c(0,1)
  47. )
  48. points(0.5,0.5)
  49. lines(table$x,table$y)
  50. # radians = 0
  51. lines(c(0.5,(0.5+1*cos(0))),c(0.5,(0.5+1*sin(0))),col='red')
  52. # radians = 1/2*pi
  53. lines(c(0.5,(0.5+1*cos(1/2*pi))),c(0.5,(0.5+1*sin(1/2*pi))),col='green')
  54. dev.off()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement