projecting GIS coordinates for analysis with spatstat package

An embedded and charset-unspecified text was scrubbed...
Name: not available
URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20090301/0bdc1e45/attachment-0002.pl>
https://answers.google.com/answers/threadview?id=577262

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.
With a conversion strategy in hand, you could relabel the results of a  
density map:

library(MASS)
mapfit <- kde2d(datat$xcoord, data$ycoord)
contour(mapfit)

The aspect ratio is "off" with this plot and the spatialstats package  
give you a less distorted map. On my Mac I can resize the quartz  
window and the lebeling remains in the proper aspect, but I don't know  
if that works for all other platforms.

So you could then multiply your ycoord by 69.10537434873548 and your x  
coord by 45.23874340706528 and get miles (from Greenwich, UK to the  
East and the Equator to the South) as your axis labels

mapdat$xmiles <- mapdat$xcoord*45.23874340706528
  mapdat$ymiles <- mapdat$ycoord*69.10537434873548
  mapfit.mi <- kde2d(mapdat$xmiles, mapdat$ymiles)
  contour(mapfit.mi)

If the conversion worked then your coordinates have about a 12 mile E- 
W range and a 6-7 mile N-S range, correct? When I do the same sort of  
on a renamed dataset, data2:

 > data2$ymiles <- data2$ycoord*69.10537434873548
 > data2$xmiles <- data2$xcoord*45.23874340706528
 > coordinates(data2) <- c("xmiles","ymiles")
 > ppp_data = as(data2["itype"], "ppp")
Warning message:
In ppp(cc[, 1], cc[, 2], window = W, marks = marks) :
   data contain duplicated points
 > density_data = density.ppp(ppp_data)
 > plot(density_data, col=brewer.pal(9, "Reds"))

I get a maximal density of 6 per <something> .  6 events per square  
mile certainly looks like a sensible estimate for that datatset. All  
you are missing now are the axis labels or perhaps an overlay of a  
local map. Sarkar's book may be helpful for that purpose.
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Event.per.sq.mile.pdf
Type: application/pdf
Size: 110759 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20090301/0254181a/attachment-0003.pdf>
-------------- next part --------------

--  
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.