I am working on creating an R package for doing fire department analysis and am trying to create a function that can display emergency incident densities. The following code sort of does the trick, but I need a display that shows the number of incidents per square mile. I believe the code below shows incidents per square unit (in this case, degrees lat/long).
To solve this problem, I believe that I need to convert the coordinates (currently WGS84) to some projection that is based on miles rather than degrees lat/long. Does anybody know the code for projecting coordinates so that my density plot will show incidents per sq-mile? If there is a simpler way of displaying incident densities than using the spatstat package, please let me know. Thanks, Markus #create data data = data.frame(xcoord=c(-123.1231, -123.0245, -123.1042, -123.1555, -123.1243, -123.0984, -123.1050, -123.0909, -123.1292, -123.0973, -123.0987, -123.1016, -123.2355, -123.1005, -123.1130, -123.1308, -123.1281, -123.1281, -123.1275, -123.1269, -123.1595, -123.1202, -123.1756, -123.0791, -123.0791, -123.0969, -123.0969, -123.0905, -123.0718, -123.0969, -123.1337, -123.1531, -123.1362, -123.1550, -123.0725, -123.1249, -123.1249, -123.1249, -123.1249, -123.1249, -123.1777, -123.1237, -123.1912, -123.0256, -123.1347, -123.1246, -123.1931, -123.0971, -123.0281, -123.0928), ycoord=c(49.27919, 49.23780, 49.24881, 49.27259, 49.26057, 49.25654, 49.25000, 49.28119, 49.27908, 49.28442, 49.28318, 49.27293, 49.25805, 49.28137, 49.22528, 49.26066, 49.27841, 49.27841, 49.28019, 49.27414, 49.24220, 49.27744, 49.23474, 49.28229, 49.28229, 49.27671, 49.27671, 49.25974, 49.26510, 49.27671, 49.29036, 49.26100, 49.27989, 49.26103, 49.27216, 49.27548, 49.27548, 49.27548, 49.27548, 49.27548, 49.23475, 49.27759, 49.24524, 49.26271, 49.20531, 49.26337, 49.23862, 49.28447, 49.20871, 49.28306), itype=c("Emergency Medical Service", "Rescue", "Service Call", "Alarm Activation", "Hazardous Condition", "Motor Vehicle Accident", "Emergency Medical Service", "Emergency Medical Service", "Fire", "Alarm Activation", "Emergency Medical Service", "Motor Vehicle Accident", "Emergency Medical Service", "Emergency Medical Service", "Emergency Medical Service", "Alarm Activation", "Alarm Activation", "Alarm Activation", "Emergency Medical Service", "Emergency Medical Service", "Emergency Medical Service", "Alarm Activation", "Emergency Medical Service", "Hazardous Condition", "Hazardous Condition", "Motor Vehicle Accident", "Motor Vehicle Accident", "Motor Vehicle Accident", "Alarm Activation", "Motor Vehicle Accident", "Emergency Medical Service", "Motor Vehicle Accident", "Alarm Activation", "Emergency Medical Service", "Emergency Medical Service", "Fire", "Fire", "Fire", "Fire", "Fire", "Motor Vehicle Accident", "Emergency Medical Service", "Emergency Medical Service", "Motor Vehicle Accident", "Alarm Activation", "Emergency Medical Service", "Alarm Activation", "Fire", "Emergency Medical Service", "Emergency Medical Service")) #add necessary libraries library(sp) library(maptools) library(spatstat) library(RColorBrewer) #add coordinates to data coordinates(data) = c("xcoord", "ycoord") #convert coordinates to spatstat point pattern dataset ppp_data = as(data["itype"], "ppp") #determine density of point pattern density_data = density.ppp(ppp_data) #plot density plot(density_data, col=brewer.pal(9, "Reds")) -- Markus Weisner, Firefighter Charlottesville Fire Department 203 Ridge Street Charlottesville, Virginia 22901 (434) 970-3240 [[alternative HTML version deleted]] ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.