[R] projecting GIS coordinates for analysis with spatstat package

[David Winsemius]

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.