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.
