Markus Weisner
2009-Mar-01 20:22 UTC
[R] projecting GIS coordinates for analysis with spatstat package
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]]
David Winsemius
2009-Mar-01 21:20 UTC
[R] projecting GIS coordinates for analysis with spatstat package
https://answers.google.com/answers/threadview?id=577262 On Mar 1, 2009, at 3:22 PM, Markus Weisner wrote:> 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 at 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.