[R] projecting GIS coordinates for analysis with spatstat package

Sun Mar 1 23:17:32 CET 2009

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

On Mar 1, 2009, at 4:20 PM, David Winsemius wrote:

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